Integrated analysis of circRNA-miRNA-mRNA network reveals potential prognostic biomarkers for radiotherapies with X-rays and carbon ions in non-small cell lung cancer.

BACKGROUND: This work was aimed at exploring the regulatory network of non-coding RNA (ncRNA) especially circular RNA (circRNA) and microRNA (miRNA), in the sensitivity of non-small cell lung cancer (NSCLC) cells to low linear energy transfer (LET) X-ray and high-LET carbon ion irradiations.
METHODS: The radioresistant NSCLC cell line A549-R11 was obtained from its parental cell line A549 through irradiation with X-rays of 2.0 Gy per fraction for 30 times. The sensitivities of A549, A549-R11 and H1299 cells exposed to X-rays and carbon ions were verified using the colony formation assay. A comprehensive circRNA-miRNA-mRNA network was constructed through the sequencing data in parental A549, acquired radioresistant A549-R11 and intrinsic radioresistant H1299 cells, and the network was further optimized according to the prognostic results from the TCGA and GEO databases.
RESULTS: Based on high-throughput sequencing of circRNAs, we found that 40 circRNAs were up-regulated while 184 circRNAs were down-regulated in the intersection of the sets of A549-R11 and H1299 cells. Subsequently, a circRNA- miRNA-mRNA network, including 14 interactive pairs and 8 circRNAs, 4 overall survival-associated miRNAs, and 4 mRNAs, was constructed through the high-throughput data screening and bioinformatics methods.
CONCLUSIONS: Our results provide a complete understanding to the regulatory mechanism of the sensitivities to low-LET X-ray and high-LET carbon ion irradiations, and might be helpful to screen potential biomarkers for predicting the Carbon-ion radiotherapy (CIRT) and X-ray radiotherapy responses in NSCLC. 2020 Annals of Translational Medicine. All rights reserved.

Introduction

Lung cancer is the most frequent primary malignance and the leading cause of cancer-related death worldwide. It has been reported that non-small cell lung cancer (NSCLC) accounts for approximately 80% of lung cancers (1). Radiotherapy as a universal treatment modality, however, produces a 5-year survival rate of merely 10–30% because of poor control of the primary tumor (2,3). Radioresistance has been considered as one of the most important reasons for treatment failure or local tumor recurrence (4). Carbon-ion radiotherapy (CIRT) is expected to be more effective even for photon-resistant tumors because of its physical and biological advantages (5,6). For early NSCLC patients, CIRT has demonstrated that local control and overall survival rates at 5 years are 80–90% and 40–50%, respectively (7). However, only a small proportion of the patients can benefit from CIRT because of high building and running costs of the accelerator facility. Therefore, it is absolutely necessary to screen potential biomarkers for predicting the CIRT and X-ray radiotherapy responses in locally-advanced NSCLC cases.

Over the past decade, with a great advance in the high-throughput sequencing, many non-coding RNAs (ncRNAs) such as microRNAs (miRNAs), long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs) have been identified and found to play essentially regulatory role in the carcinogenesis (8-11). CircRNA, unlike the well-known linear RNA, forms a ring structure through the linkage of the 3′ and 5′ ends with a covalent bond (12). Recently, increasing studies have revealed that circRNAs play regulatory roles mainly at the transcriptional and post-transcriptional levels, including acting as miRNA sponges (13), interaction with RNA-binding proteins (RBPs) (14), and the modulation of gene transcription (15). The recently uncovered circRNA has protein-coding potential in eukaryotes (16). A few circRNAs, such as circRNA 100876 (17), hsa_circ_0043256 (18), circRNA HIPK3 (19), have been reported to be involved in the regulation of NSCLC progression.

In present study, we developed a lung carcinoma cell line (A549-R11) by irradiating human lung carcinoma A549 cells with fractionated irradiation (2 Gy × 30 fractions), which is resistant to low- linear energy transfer (LET) X-rays. Moreover, these two cell lines (A549 and A549-R11) showed a similar sensitivity to high-LET carbon ions. To discover the underlying molecular regulation mechanisms of circRNAs, miRNAs, and mRNAs in the radiosensitivity, the next-generation sequencing (NGS) technology was used to detect the differentially expressed circRNAs among the radioresistant NSCLC cell line A549-R11, its parental cell line A549 and an intrinsic radioresistant NSCLC cell line H1299. The differential expressions of circRNAs might be useful predictors for the differential responses of NSCLC to low-LET X-rays and high-LET carbon ions.

Methods

Establishment and identification of radioresistant cell line

Human NSCLC cell lines A549 and H1299 were purchased from the Type Culture Collection of the Chinese Academy of Sciences (Shanghai, China) and cultured in RPMI 1640 Medium (Thermo Fisher Scientific Inc.) containing 10% heat-inactivated fetal bovine serum (FBS; Bailing Bio, Lanzhou, China) at 37 °C in a humidified 5% CO2 atmosphere. The radioresistant NSCLC cell line A549-R11 was obtained from its parental cell line A549 through irradiation with X-rays of 2.0 Gy per fraction for 30 times (the process was depicted in ).

Figure 1

The flow charts of the establishment of the radioresistant NSCLC cell line and the construction of circRNA-miRNA-mRNA network. (A) The screening procedure of the radioresistant cell line A549-R11. (B) The flowchart for constructing the circRNA-miRNA-mRNA network. DE, differentially expressed; MREs, miRNA response elements; OS, overall survival.

The sensitivities of A549, A549-R11 and H1299 cells exposed to X-rays and carbon ions were verified using the colony formation assay (20).

RNA extraction and identification of radioresistance-associated circRNAs, miRNAs and mRNAs

Three samples were obtained from each of the three cell cultural groups (A549, A549-R11 and H1299 cells). Total RNA was isolated using the TRIzol (Invitrogen, USA) according to the manufacturer’s protocol. All the sequencing procedures and analyses were performed in Ribobio (Guangzhou, China). Differentially expressed circRNAs, miRNAs and mRNAs between two groups (A549 cells vs. A549-R11 cells and A549 cells vs. H1299 cells) with statistical significance were defined as |fold change| ≥2 and P<0.05 estimated with the Student t-test.

qRT-PCR validation

qRT-PCR was used to confirm the expression of differentially expressed circRNAs obtained from the sequencing data. Total RNA was extracted and assays were done in triplicate using the primers indicated in Table S1.

Kaplan-Meier analysis

To investigate the impact of the differentially expressed miRNAs and mRNAs on prognostic overall survival (OS) of patients with lung cancer, Kaplan-Meier survival analysis was performed using OncoLnc (http://www.oncolnc.org/), which is a database for interactively exploring survival correlations and contains survival data for more than 8,000 patients from 21 cancer studies by The Cancer Genome Atlas (TCGA) (21). The OSs of patients based on the expression level of the miRNAs and mRNAs in the A549-R11 set were compared. Briefly, patients were classified into the high or low expression group based on the expression level of each differentially expressed miRNA or mRNA, and then the OS was analyzed.

To further confirm the above-mentioned results, another tool, Kaplan-Meier plotter was employed (22). The Kaplan-Meier plotter (http://kmplot.com/analysis/) is capable to assess the effect of 54,000 genes on survival in 21 cancer types. In particular, we focused on a cohort of patients with NSCLC (GSE31210), which includes 246 patients, to evaluate the prognostic value of the differentially expressed mRNAs.

Prediction of circRNA-miRNA-mRNA networks

CircRNA, serving as a miRNA sponge to bind miRNA competitively, indirectly regulates target gene, which is known as competing endogenous RNA (ceRNA) (23). Based on this theory, circRNA-miRNA-mRNA networks were constructed. The interactions between miRNAs and circRNAs were predicted through miRNA response elements (MREs) in the miRanda, RNA hybrid, and TargetScan databases. Moreover, these circRNA-miRNA pairs were further screened using the OncoLnc toolkit according to the correlation between the differentially expressed miRNAs and OS of patients with lung adenocarcinoma (LUAD). Next, we used miRDB, an online database for miRNA target prediction and functional annotations (http://www.mirdb.org/), to predict target mRNAs. These interactions were similarly screened using the OncoLnc and Kaplan-Meier plotter toolkits. The flowchart of the circRNA-miRNA-mRNA interaction construction is depicted in and the network was visually displayed using the Cystoscope software V3.7.1 (San Diego, CA, USA).

Statistical analysis

Data are represented as the mean ± standard deviation (SD). Statistical analysis was conducted using the unpaired Student’s t-test. A difference was considered significant when P<0.05.

Results

Radiosensitivities of A549, A549-R11, and H1299 cells to low- and high-LET radiations

Radioresistance was identified by comparing the survival fraction of A549-R11 or H1299 with A549 cells. Shown in are the survival curves of the three cell lines, after exposure to X-rays. A549 cells represented the most significant sensitivity among the three cell lines, whereas A549-R11 and H1299 cells showed obvious radioresistance. For carbon ion irradiation (), H1299 exhibited the most significant radioresistance, while A549 and A549-R11 cells displayed a similar cell survival curve. The survival fractions at 2Gy of the three cell lines exposed to X-rays and carbon ions are summarized in . A549-R11 and H1299 cells were more resistant to low-LET X-rays than A549 cells; however, A549-R11 cells possessed a similar sensitivity to high-LET carbon ions like their paternal A549 cells.

Figure 2

Survival fractions of A549, A549-R11 and H1299 cells exposed to low-LET X-rays (A) and high-LET carbon ions (B) measured with the clonogenic survival assay. *, P<0.05; **, P<0.01 vs. A549 cells.

Table 1

Summary of survival fraction and RR at 2 Gy

Treatment	X-rays	RR	Carbon ions	RR
SF2 (A549)	0.393	–	0.117
SF2 (A549-R11)	0.568	1.445	0.122	1.043
SF2 (H1299)	0.523	1.331	0.358	3.060

The data calculated from the survival curves. RR, radiosensitivity ratios.

Overview of circRNAs profiles

High-throughput next-generation RNA sequencing data were analyzed to explore the circRNA expression profile in radiosensitive A549, acquired radioresistant A549-R11 and intrinsic radioresistant H1299 cell lines. Using the CIRI2 and CIRCexplorer2 software, 28,082 distinct circRNA candidates containing at least one unique back- spliced read were found in the three cell lines. These circRNAs and their host genes are located in various genomic regions (Figure S1). The expression analysis of these transcripts revealed that a series of circRNAs were differentially expressed in the acquired radioresistant A549-R11 cells and intrinsic radioresistant H1299 cells compared with the radiosensitive A549 cells (). In the A549-R11 set, a total of 557 circRNAs were differentially expressed (|fold change| ≥2 and P<0.05), consisting of 164 up-regulated circRNAs and 393 down-regulated circRNAs in A549-R11 cells compared with A549 cells. In the H1299 set, the expression profiles of 2,056 circRNAs were different between H1299 and A549 cells, where 879 circRNAs were up-regulated while 1,177 circRNAs were down-regulated in H1299 cells. Then the different expressions were taken at the intersection of the sets of A549-R11 and H1299 (). We found that 40 were up-regulated and 184 were down-regulated in both sets, whereas 17 were inconsistent in two sets. The top twenty up- and down-regulated circRNAs ranked by fold change in an intersected set are shown in a hierarchical clustering map (). The distribution of the differentially expressed circRNAs on the chromosomes and the fold change are exhibited in .

Figure 3

CircRNA expression profile. Volcano plots exhibit significantly dysregulated miRNAs in A549-R11 (A) and H1299 (B) sets, where the horizontal lines represent 2-fold up and down expressed circRNAs, and the vertical lines indicate P=0.05. (C) The intersection was taken in the different expressions of A549-R11 and H1299 sets. (D) Hierarchical clustering shows the top twenty up-regulated and down-regulated circRNAs. (E) The distribution of the differentially expressed circRNAs on the chromosomes and the fold change. The green and red dots represent the down- and up-regulated circRNAs, respectively. Center distance represents the fold change.

Validation of circRNA expressions

To verify the sequencing data, twenty differentially expressed circRNAs were selected, including the top 10 up-regulated circRNAs and the top 10 down-regulated circRNAs in the paired radioresistant A549-R11 and radiosensitive A549 cell lines. The results demonstrated that the expression levels of the circRNAs were consistent with the sequencing data, except hsa_circ_0084606, which was down-regulated in the sequencing data but up-regulated in the qPCR experiment (). These results suggest the high reliability of the sequencing data.

Figure 4

Validation of the sequencing data. Twenty differentially expressed circRNAs including 10 up-regulated (A) and 10 down-regulated (B) were validated by qRT-PCR.

Expression profile of miRNAs and construction of OS-related circRNA-miRNA network

In our sequencing data, 399 miRNAs were identified to be differentially modulated. Among them, 314 miRNAs were down-regulated whereas 85 miRNAs were up-regulated (). The MiRanda, RNA hybrid, and TargetScan databases were used to predict the circRNA-miRNA pairs. We found that 81 of the above-mentioned total 224 intersected circRNAs had target relationships with 277 of the 399 differentially expressed miRNAs and formed 934 circRNA–miRNA pairs (Figure S2). However, these data were too large to further analysis. So, OS as a prognostic factor about differentially expressed miRNAs was introduced to screen these pairs. We found 17 miRNAs among the 399 differentially expressed miRNAs were significantly correlated with the survival of patients (Table S2). shows an intersection of the two miRNA groups (targeted miRNAs and OS related miRNAs). There were 7 miRNAs consisting of 6 down-regulated miRNAs and 1 up-regulated miRNA in this intersection. Their Kaplan-Meier survival analysis is shown in Figure S3. Finally, a circRNA-miRNA related OS network was constructed, including 11 circRNAs, 7 miRNAs, and 12 pairs ().

Figure 5

Construction of the OS-related circRNA-miRNA networks. (A) A scatter plot shows the sequencing data distribution of miRNAs between A549-R11 and A549 cells. The red and green points show that the miRNAs were up-regulated and down-regulated, respectively. (B) The intersection was taken between the OS-related miRNAs and targeted-miRNAs in NSCLC. (C) The circRNA-miRNA interaction network associated with OS. Blue symbols represent up-regulated expressions, whereas yellow symbols represent down-regulated expression. Triangular and quadrilateral nodes represent miRNAs and circRNAs, respectively.

Differentially expressed mRNA profiles and construction of circRNA-miRNA- mRNA network correlated with the prognosis of NSCLC

A total of 661 differentially expressed mRNAs were identified between A549 and A549-R11 cells, including 406 up-regulated and 255 down-regulated mRNAs (). The GO analysis revealed that, for these altered mRNAs, the top 3 enriched terms were cell part, cell and binging between A549 and A549-R11 cells. The KEGG analysis showed the top 10 pathways associated with the different mRNAs in the A549-R11 group (). Among the pathways, we identified the PI3K-Akt signaling pathway, cell adhesion molecules (CAMs) and the Ras signaling pathway, which have been reported to be related to tumorigenesis, progression, invasion and metastasis and drug resistance.

Figure 6

A regulatory network of circRNA-miRNA-mRNA associated with prognosis. (A) A volcano plot exhibits significantly dysregulated mRNAs between A549-R11 and A549 groups. The red points show that the miRNAs were up-regulated more than two-fold. The green points show that the miRNAs were down-regulated more than two-fold. (B) Enrichment of the top 10 GO terms (left) and KEGG pathways (right) of differentially expressed mRNAs in A549-R11 and A549 groups. (C) The Kaplan-Meier analysis of NSCLC-specific overall survival of patients with tumors expressing different levels of mRNAs. (D) A circRNA-miRNA-mRNA network correlated with the prognosis of NSCLC. Quadrilateral, triangular and rectangle nodes represent circRNAs, miRNAs and mRNAs, respectively.

As the above-mentioned 7 prognostic signatures and circRNA-targeted miRNAs were employed to predict miRNA-mRNA interactions, we found that 32 of the 661 differentially expressed mRNAs could form 58 pairs with the miRNAs, 27 of which were up-regulated and 5 of which were down-regulated as shown in Table S3 and Figure S4. Likewise, prognostic survival analysis was conducted for the 32 mRNAs to examine whether the mRNAs expression significantly correlate to the survival of patients with NSCLC using the Oncolnc and Kaplan-Meier plotter datasets. However, only 4 mRNAs (COL5A2, FBN2, FCDH7 and PXDN) were found to be significantly associated with OS, both indicating worse prognosis in the high expression group (). Subsequently, a circRNA-miRNA-mRNA visualized regulatory network was constructed using the Cytoscape 3.7.1 (), containing 14 circRNA-miRNA-mRNA interactive pairs and 8 circRNAs, 4 OS-associated miRNAs, and 4 mRNAs.

Discussion

Radioresistance remains a significant clinical problem for radiotherapy with conventional radiation such as X-rays, because various tumors display widely varying radiosensitivities. To explore the radiation response in NSCLC, we made an isogenic model of radioresistance via chronic exposure to fractionated radiations, which has been previously reported (24-28). Compared with cell lines of the same origin but with obviously different radiosensitivities, the advantage of this type of model is that it avoids the influence of confusing factors such as genetic and inherent variation, allowing for the clarification of molecules that are involved in radioresistance. Our results also verified this point. The number of different circRNA expressions in the A549-R11 set accounted for merely 27% of those in the H1299 set. Additionally, this model mimics what occurs in the clinical setting, because standard radiotherapy consists of about 2 Gy once a day, 5 days a week, a period for 5–8 weeks (26). This produced a subclone cell line A549-R11, which exhibited enhanced clonogenic survival after irradiation compared to its parental cell line A549. Although A549-R11 cells were specifically resistant to low-LET X-rays, they possessed a similar sensitivity to high-LET carbon ions like their parental A549 cells. Besides, we also introduced another NSCLC cell line H1299, which has intrinsic resistance to X-rays and carbon ions compared with A549 cells, into our study to verify our sequencing data from the A549-R11 set.

It has been reported that radiation could induce a change of circRNA expressions in tumor and normal cells (29-31). Additionally, circRNAs could regulate cellular radiosensitivity. CircRNA_014511 over-expression has been linked to radioresistance in bone marrow mesenchymal stem cells (32). Shuai et al. (33) demonstrated that circRNA_0000285 located in the HIPK3 gene locus involved in the radiosensitivity of nasopharyngeal carcinoma. Different from the studies mentioned above, we believe that some circRNAs could be taken as the biomarkers for NSCLC diagnosis and prognosis in CIRT or radiotherapy with conventional radiation (X- or γ-rays). In this study, we detected the circRNA expressions among the different radiosensitive NSCLC cell lines (including acquired and intrinsic radioresistance) rather than compared the expression changes before and after irradiation. This may produce more realistic environment and improve diagnostic accuracy. We identified a large number of differentially expressed circRNAs among the cell lines from various genomic locations. The circRNA sequencing data revealed that 40 circRNAs were significantly up-regulated and 184 circRNAs were down-regulated in the radioresistant cells compared with their parental cells. Moreover, there was another characteristic for the radioresistant A549-R11 cell line, that is it had a similar sensitivity to carbon ions compared with its parental A549 cell line. This suggests that the circRNAs with different expressions between the two cell lines might be potential biomarkers to distinguish the radiation response after irradiation with X-rays or carbon ions. In fact, we have already obtained a circRNA hsa_circ_0103301, which is highly expressed in A549-R11 and H1299 cells. It could regulate the sensitivity to X-rays but has no effect on NSCLC cells after carbon ion irradiation (unpublished data).

Generally, interruption of a single ncRNA may only regulate a suppression effect, whereas joint destruction of multiple ncRNAs may be more effective. Based on this viewpoint, an interfering circRNA targeting multiple miRNAs was proposed and subjected to trials for NSCLC radiotherapy. Firstly, 934 circRNA-miRNA pairs and 58 miRNA-mRNA pairs were individually predicted using the bioinformatics method. Next, given the significant OS about differentially expressed miRNAs and mRNAs, only 12 circRNA-miRNA pairs and 5 miRNA- mRNA pairs, were selected to generate a global regulatory network. Finally, a circRNA-miRNA-mRNA network including 8 circRNAs, 4 miRNAs, and 4 mRNAs was constructed based on negative correlations between miRNA and mRNA as well as miRNA and circRNA.

In the network we constructed, some ncRNAs and mRNAs have been reported to probably be involved in the initiation and progression of cancer. For example, Zheng et al. found that PXDH was highly expressed in ovarian cancer, which promoted proliferation, invasion and migration and was associated with poor prognosis (34). The similar results were observed in metastatic melanoma tumor (35). PCDH7 upregulation was observed in castrate-resistant prostate cancer (CRPC) and CRPC cells. The knockdown of PCDH7 decreased the activity of the PI3K/Akt pathway and inhibited cell proliferation (36). Recently, several studies showed that there was clinical significance of COL5A2 in patients with adenomas (37), breast cancer (38), colorectal cancer (39), and osteosarcoma (40), especial bladder cancer (41-43). Furthermore, some studies reported that the miRNAs in the constructed network might be differentially expressed in some cancers. Ma et al. uncovered 15 miRNAs, including hsa-miR-16-1-3p, were down-regulated in the plasma of EGFR-TKI resistant patients with NSCLC compared with that of EGFR-TKIs sensitive patients (44). Hsa-miR-16-1 was also mostly down-regulated in lung adenocarcinoma cell lines, and induced cell cycle arrest (45). It has been reported that over-expression of miR-34b could decrease proliferation, migration, and invasion and promote apoptosis in lung cancer (46) and colon cancer cells (47) by targeting different pathways. Lower miR-34b was significantly associated with an aggressive phenotype in urothelial carcinoma of the bladder patients compared with the nonaggressive subject (48). Moreover, down-regulated hsa-miR-146b-3p was significantly correlated with recurrence in osteosarcoma (49), and OS in hepatocellular carcinoma (22). MiR-449 mimics strongly inhibited proliferation, promotes apoptosis, and leads to cell cycle arrest in prostate cancer (50,51) and hepatocellular carcinoma cells (52). Also, we searched only one paper about 4 circRNAs in our network, namely up-regulated has_circ_0001944 may be involved in breast cancer brain metastasis (53). Other circRNAs detected in our ceRNA network have not been reported previously.

In addition, another cohort, CAARRY including 68 patients who received radiotherapy in the Kaplan-Meier plotter, was used to verify the 4 mRNAs mentioned above. CLO5A2 and PXDN were observed to have a significant correlation with the first progression of the lung cancer patients after radiotherapy (Figure S5), suggesting the accuracy of our prediction clearly.

Conclusions

We united high-throughput data screening and bioinformatics to construct a ncRNA network describing the possible regulatory mechanisms of NSCLC sensitivity to low-LET X-rays and high-LET carbon ions. Importantly, our results provide promising potential molecular markers for distinguishing the patients with NSCLC who will gain from radiotherapy with conventional X-rays or CIRT.

The article’s supplementary files as