Comprehensive circular RNA expression profile of lung adenocarcinoma with bone metastasis: Identification of potential biomarkers.

Background: Lung adenocarcinoma (LUAD) has a significant tendency to metastasize to the bone, with severe comorbidities. Recent studies have reported that circular RNAs (circRNAs) are involved in various cancer metastasis-related physiological cellular processes. However, their role in LUAD with bone metastasis (LUAD-BM) remains unknown.
Methods: Bone metastasis (BM) circRNAs were identified using high-throughput sequencing and validated by quantitative reverse transcription-PCR (qRT-PCR). Bioinformatic analyses were used to predict the potential functions of the differentially expressed circRNAs. The effects of circ_0096442 on the growth and metastasis of A549 cells were detected in a co-culture system of A549 and bone marrow-derived cells.
Results: There were 598 (238 upregulated and 360 downregulated) 390 (187 upregulated and 203 downregulated) and 644 (336 upregulated and 308 downregulated) differentially expressed circRNAs between LUAD-BM and LUAD, LUAD-BM and healthy individuals, and LUAD and healthy individuals, respectively. These differentially expressed circRNAs play important roles in cellular components, biological processes, and molecular functions. Moreover, they map several pathways related to BM, including DNA repair, DNA damage, and osteoclast differentiation. The results validated by qRT-PCR for the five most dysregulated circRNAs are consistent with the sequencing data. Additionally, circ_0096442 was found to promote the growth and metastasis of LUAD in a bone microenvironment.
Conclusion: Our findings provide a novel and important circRNA expression profile of LUAD-BM and suggest that circ_0096442 may be a biomarker for LUAD-BM.

Introduction

Bone is the third most common target organ for advanced cancer metastasis (Coleman, 2006). In advanced lung cancer, 30%–40% of patients develop bone metastasis (BM) with severe comorbidities (Cho et al., 2019; Segaliny et al., 2019). The presence of BM has serious consequences on the quality of life of patients, and dramatically reduces the overall survival rate of affected individuals (Borghaei et al., 2015). Despite the advancements in diagnostic and therapeutic strategies, BM detection at an early stage and its poor prognosis for lung cancer remain challenging. Thus, there is an urgent need to better understand BM pathogenesis to identify specific biomarkers and develop new treatment strategies for preventing and treating BM.

Genomic studies have indicated that each step of BM is related to a series of molecular events (Liu et al., 2014). It is currently accepted that microRNAs (miRNAs) play important roles in regulating BM physiological processes and serve as potential biomarkers (Hesse and Taipaleenmaki, 2019). With the advancement of high-throughput sequencing technology, circular RNAs (circRNAs) have been discovered to be pervasively expressed in human genes (Li and Han, 2019). Compared to miRNAs, circRNAs are more stable and not easily degraded by exonucleases. CircRNAs are involved in various cancer-related physiological cellular processes and are recognized as potential biomarkers in cancer metastasis (Li and Yang, 2021; Wei et al., 2021; Yarmishyn et al., 2022).

As the leading cause of cancer morbidity and mortality, lung cancer is usually diagnosed after developing locally or systematically (Bray et al., 2018; da Silva et al., 2019), with lung adenocarcinoma (LUAD) as the most common subtype (Torre et al., 2015; Yang et al., 2018). Among BM patients, cases without a history of cancer account for 25%–30%, and the diagnostic assessment may be delayed due to difficulties in differentiating BM from orthopedic degenerative diseases (Kitagawa et al., 2019). Therefore, this study simultaneously presents high-throughput sequencing analyses showing differential circRNA expression profiles and their regulatory interaction networks from three pairs of groups, which are LUAD with bone metastasis (LUAD-BM) vs. LUAD, LUAD-BM vs. healthy individuals, and LUAD vs. healthy individuals. Moreover, the effects of circ_0096442 (a potential identified biomarker) on the growth and metastasis of A549 cells were detected in a co-culture system of A549 and bone marrow-derived cells.

Materials and methods

Clinical samples

The study was approved by the ethics committee of our hospital, and all participants provided consent for blood donation. All LUAD patients were diagnosed based on histopathologic tests and had no history of other tumors. The emission-computed tomography of bone scintigraphy showed multiple lesions in patients with BM. Control patients were matched in age, sex, body mass index, and region (characteristics of the patients are shown in Supplementary Table S1).

RNA extraction and high-throughput sequencing

Blood samples were collected into anticoagulant tubes, mixed with exactly three volumes of RNASafer Reagent (Magen) and stored under strict standard operating procedure conditions temporarily at −80°C before batch analysis. Hipure PX Blood RNA Kits (Magen) were used to extract total RNA. Concentrations of isolated RNA were detected using a Qubit 3.0 fluorometer (Invitrogen, Carlsbad, California). Isolated RNA integrity was detected with an Agilent 2100 bioanalyzer (Applied Biosystems, Carlsbad, CA). KAPA RNA HyperPrep Kits with RiboErase for Illumina (Kapa Biosystems, Woburn, MA) were used to prepare a high-throughput sequencing library with 2 μg of total RNA. Illumina Hiseq X ten was used to perform paired-end (PE150) sequencing.

Quantitative reverse transcription-PCR validation

Twenty-five pairs of clinical samples from patients with LUAD-BM or LUAD were used to validate five candidate circRNAs by quantitative reverse transcription-PCR (qRT-PCR). Primer sequences and back-splicing sites of the candidate circRNAs are shown in Supplementary Table S2 and Supplementary Figure S1, respectively.

Cell culture and co-culture system

The human lung adenocarcinoma (HLA) A549 and bone marrow stromal HS-5 cells were purchased from the American Type Culture Collection (ATCC; Manassas, VA, United States) and cultured in Dulbecco’s Modified Eagle’s medium (DMEM; Hyclone, United States) containing 10% fetal bovine serum (FBS; Gibco, United States) and 1% penicillin/streptomycin at 37°C in a 5% CO2 incubator under humidified conditions. The co-culture system was set up using 6-well Transwell inserts (Corning, United States) with a 0.4-µm pore size.

Dual-luciferase reporter assay

Circ96442-wt and circ96442-mut were inserted into the psiCHECK2 dual-luciferase vector (Promega, Madison, WI, United States). Hsa-miR-326 mimic and negative control miRNA (miR-NC) were synthesized by GenePharma (GenePharma, Shanghai, China). After co-transfection of the reporter vector and miR-326 mimic or negative control in A549 cells for 48 h, Renilla luciferase activity was measured using a dual-luciferase assay kit (Promega, Madison, United States) against that of Firefly luciferase.

Fluorescence in situ hybridization

A549 cells were seeded on cell slides at the bottom of a 24-well plate and fixed with 4% paraformaldehyde. The CY3-labeled probe targeted the hsa_circ_0096442 junction site, and the FITC-labeled hsa-miR-326 probes were designed and synthesized by Geneseed (Guangzhou, China). Nuclei were stained with 4,6-diamidino-2-phenylindole (DAPI). The signals of the probes were detected using a fluorescence in situ hybridization (FISH) kit (GenePharma, Shanghai, China) according to the manufacturer’s instructions.

Statistical analysis

All values are expressed as the mean ± standard error of the mean (SEM). Statistical significance was determined using the Student’s t-test. Statistical analyses were performed using GraphPad Prism 8 (GraphPad Software, Inc., San Diego, CA, United States) and SPSS v.20.0 (SPSS Inc., Chicago, IL, United States). Differences with p < 0.05 were considered as statistically significant and were noted by asterisks (*, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001).

More method details (have been provided in Additional file 1).

Results

Differentially expressed circular RNAs among the different groups

There were 598 (238 upregulated and 360 downregulated) 390 (187 upregulated and 203 downregulated) and 644 (336 upregulated and 308 downregulated) differentially expressed circRNAs between LUAD-BM and LUAD, LUAD-BM and healthy individuals, LUAD and healthy individuals, respectively. Notably, 54 consistently dysregulated circRNAs were identified in the LUAD-BM group than in both control groups, including 19 previously unknown circRNAs (Supplementary Table S3). Differentially expressed circRNAs among the different groups were illustrated by hierarchical clustering analyses (Figure 1) and volcano plots (Figure 2).

FIGURE 1

Hierarchical clustering analysis of differentially expressed circRNAs among different groups. (A) The differentially expressed circRNAs of LUAD-BM, LUAD, and control. (B) The differentially expressed circRNAs between LUAD-BM and LUAD. (C) The differentially expressed circRNAs between LUAD-BM and healthy individuals. (D) The differentially expressed circRNAs between LUAD and healthy individuals. LUAD-BM: 201, 202, 203; LUAD: 101, 102, 103; Control: 001, 002, 003. LUAD-BM, lung adenocarcinoma with bone metastasis; LUAD, lung adenocarcinoma; circRNAs, circular RNA.

FIGURE 2

Volcano plots of differentially expressed circRNAs in peripheral blood among different groups. (A) The differentially expressed circRNAs between LUAD-BM and LUAD. (B) The differentially expressed circRNAs between LUAD-BM and healthy individuals. (C) The differentially expressed circRNAs between LUAD and healthy individuals. Red and green dots indicate upregulated and downregulated circRNAs, respectively. LUAD-BM, lung adenocarcinoma with bone metastasis; LUAD, lung adenocarcinoma; circRNAs, circular RNA.

Functional prediction of differentially expressed circular RNAs

The most significant enrichment pathways of Gene Ontology (GO) Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome analyses of differentially expressed circRNAs between LUAD-BM and LUAD, LUAD-BM and healthy individuals, and LUAD and healthy individuals are shown in Figures 3–5, respectively. Several pathways, such as DNA repair, DNA damage, and osteoclast differentiation, were related to BM and may be significantly dysregulated in LUAD-BM. In addition, a regulatory network of circRNA-miRNA-mRNA was constructed based on 35 consistently dysregulated circRNAs (excluding 19 previously unknown circRNAs) in the LUAD-BM group compared with both control groups (Supplementary Figure S2).

FIGURE 3

Functional predictions of differentially expressed circRNAs between LUAD-BM and LUAD. (A) Top 10 GO analysis to demonstrate the molecular functions of dysregulated circRNAs. (B) Top 15 significant enrichment KEGG pathways of dysregulated circRNAs. (C) Top 15 significant enrichment Reactome pathways of dysregulated circRNAs. LUAD-BM, lung adenocarcinoma with bone metastasis; LUAD, lung adenocarcinoma; circRNAs, circular RNA.

FIGURE 5

Functional predictions of differentially expressed circRNAs between LUAD and healthy individuals. (A) Top 10 GO analysis to demonstrate the molecular functions of dysregulated circRNAs. (B) Top 15 significant enrichment KEGG pathways of dysregulated circRNAs. (C) Top 15 significant enrichment Reactome pathways of dysregulated circRNAs. LUAD, lung adenocarcinoma; circRNAs, circular RNA.

Functional predictions of differentially expressed circRNAs between LUAD-BM and healthy individuals. (A) Top 10 GO analysis to demonstrate the molecular functions of dysregulated circRNAs. (B) Top 15 significant enrichment KEGG pathways of dysregulated circRNAs. (C) Top 15 significant enrichment Reactome pathways of dysregulated circRNAs. LUAD-BM, lung adenocarcinoma with bone metastasis.

Validation by quantitative reverse transcription-PCR

Five candidate circRNAs were validated by qRT-PCR to verify the RNA-seq data. Expression patterns consistent with the sequencing data were found (Figure 6). Among these five circRNAs, circ_0096442 had the most significant differential expression.

FIGURE 6

Validation of the five-candidate circular RNAs. The relative levels of the five circRNAs is confirmed by qRT–PCR. LUAD-BM, lung adenocarcinoma with bone metastasis; LUAD, lung adenocarcinoma; circRNAs, circular RNA.

Circ_0096442 promotes proliferation and migration/invasion of A549 cells in the bone microenvironment

To further confirm the function of circ_0096442 in lung cancer, we constructed stable A549 cells overexpressing circ_0096442 and A549 cells expressing short hairpin RNA (shRNA); the efficiency of circ_0096442 overexpression and knockdown was verified by qPCR (Figure 7A). A549 cells were co-cultured with HS-5 cells, and the proliferation, migration/invasion, and apoptosis of A549 cells in the co-culture system were examined. Cell Counting Kit-8 (CCK8) and colony formation assays showed that circ_0096442 overexpression promoted cell proliferation, whereas its knockdown significantly reduced the proliferative capabilities of A549 cells (Figures 7B,C). We evaluated the migration and invasion using a Transwell assay and found that circ_0096442 overexpression promoted migration and invasion of A549 cells, whereas its knockdown inhibited this (Figure 7D). In addition, compared to the control groups, circ_0096442 overexpression suppressed the apoptosis of A549 cells, whereas its knockdown promoted it. (Figure 7E). Circ_0096442 overexpression significantly inhibited the protein expression levels of pro-apoptosis-related proteins [Caspase3 and BCL2-associated X protein (BAX)] and significantly promoted the protein expression levels of anti-apoptosis-related gene BCL2. Contrasting results were obtained from studies on circ_0096442 knockdown cells (Figure 7F). Furthermore, wound healing assays found that circ_0096442 overexpression significantly promoted the wound healing ability, whereas its knockdown significantly decreased this healingability (Figure 7G). In summary, circ_0096442 overexpression promotes the proliferation, invasion, and migration of A549 cells in the bone microenvironment.

FIGURE 7

Circ_0096442 promotes proliferation and migration/invasion of A549 cells in bone microenvironment. (A) RT-qPCR detection showing the expression of circ_0096442 in overexpression and knockdown stable A549 cells. (B) Cell proliferation was detected by CCK8 assays. (C) Cell proliferation was detected by colony formation assays. (D) Transwell assays were applied to discover the migration/invasion ability of A549 cells in bone microenvironment. (E) Apoptosis assay by flow cytometric analysis. (F) Western blot analysis for the detection of pro-apoptotic and anti-apoptotic proteins. (G) Impact of circ_0096442 on A549 cellular migration assessed through wound healing assays.

circ_0096442 can bind to hsa-miR-326 and act as its sponge

The binding sites between circ_0096442 and hsa-miR-326 were predicted by the miRanda algorithm. Next, we designed and generated psiCHECK2 vectors carrying wild-type and mutant miRNA-target sites in the circ_0096442 sequence, enabling validation of the direct binding of miR-326 and circ_0096442 (Figure 8A). The dual-luciferase assay demonstrated that miR-326 inhibited luciferase activity with wt-circ_0096442 co-transfection compared to mimic NC control but did not influence luciferase activity with mut-circ-0096442 (Figure 8B). Additionally, FISH results showed that circ_0096442 was co-localized with miR-326 in the cytoplasm of A549 cell lines (Figure 8C). These results revealed that circ_0096442, through its expected binding site, functions as a miR-326 sponge.

FIGURE 8

Circ_0096442 acted as a sponge of miR-326 in A549 cells. (A) Schematic representation of potential binding sites between circ_0096442 and miR-326. (B) The luciferase activity of the WT psicheck2-circ_0096442 or Mut psicheck2-circ_0096442 after transfection with miR-326 mimics in A549 cells. (C) FISH showed the subcellular co-localization between circ_0096442 (red) and miR-326 (green). DAPI (blue) was used to stain the nuclear.

Discussion

As microarray analyses can only identify hundreds of circRNAs, this study adopted high-throughput sequencing to analyze all circRNAs attained from blood samples of patients with LUAD-BM, LUAD alone, or healthy individuals. After strict filtering of approximately 270 million reads, 40,207 circRNAs were detected, including 598 differentially expressed circRNAs between LUAD-BM and LUAD, 390 differentially expressed circRNAs between LUAD-BM and healthy individuals, and 644 differentially expressed circRNAs between LUAD and healthy individuals. These circRNAs may play important roles in regulating the pathological processes of BM.

Recently, miRNAs have become a subject of interest in lung cancer with BM, mainly presenting as osteolytic bone destruction due to enhanced osteoclast activity. For example, Valencia et al.. (2013) found serum levels of miR-326 strongly linked to bone turnover markers in lung cancer with BM. Moreover, another study found that exosome-bound miR-192 significantly aggravated BM in immunocompromised mice, suggesting that secreted factors derived from miR-192 overexpressed lung cancer cells promote osteolytic lesions and bone colonization (Valencia et al., 2014). Using human lung cancer cell lines in vivo and in vitro, Gong et al. (2014) showed that miR-355 overexpression significantly slows down osteoclast induction and tumor cell invasion. Further, Xu et al. (2018) found that miR-139-5p expression in serum from patients with LUAD-BM was significantly lower than that of patients with metastases at other sites. In a study conducted by Thomas et al. (2012), miR-33a downregulation in A549 lung cancer cells reduced osteolytic BM. Increasing evidence has demonstrated that circRNAs absorb miRNAs as competing endogenous RNAs to inhibit their function, thus indirectly promoting the downstream function of the target gene (Zheng et al., 2016; Cheng et al., 2018). Moreover, circRNA expression changes can directly influence the expression of their source gene (Memczak et al., 2013).

The molecular mechanisms underlying the interactions between circRNAs, miRNAs and mRNAs in BM are still unclear. Therefore, a circRNA-miRNA-mRNA network of BM was constructed based on 35 consistently dysregulated circRNAs in patients with LUAD-BM compared with both control groups. The results indicate that the identified circRNAs could potentially interact with the miRNAs in the aforementioned studies. For example, according to our coexpression analysis, miR-326 can potentially bind to hsa_circ_0111889,hsa_circ_0096438,hsa_circ_0000657 and hsa_circ_0058988.

Although differentially expressed circRNAs in patients with BM were identified, the underlying mechanisms remain poorly understood. Each molecular pathway contributing to BM is regulated by various factors, through closely-controlled genes expressed by cancer cells and interactions between cells within the bone microenvironment (Ell and Kang, 2012). Functional annotations of target genes can predict the biological functions of the differentially expressed circRNAs. This study found that several pathways related to BM, such as DNA repair, DNA damage, and osteoclast differentiation, may be significantly dysregulated during LUAD-BM. Recent studies have demonstrated that DNA damage and repair are closely related to BM (Iglesias-Gato et al., 2018; Isaacsson Velho et al., 2019). Moreover, tumor cells colonizing the bone can promote osteoclast differentiation and osteoclastic bone resorption (Nakai et al., 2019). These findings are in accordance with our enrichment analyses, suggesting that these circRNAs might be involved in BM development and progression.

Bone is the most common metastatic site for primary breast and lung tumors (Suva et al., 2011). However, information about circRNAs during BM is highly limited to date. In a recent study, Xu et al. (2021) screened breast cancer bone-metastatic circRNAs using deep sequencing, validated the results using in situ hybridization, and demonstrated that the circular inhibitor of nuclear factor-kappa B kinase subunit beta (circIKBKB) was upregulated significantly. They further revealed a plausible mechanism for circIKBKB-mediated nuclear factor-kappa B (NF-κB) hyperactivation in bone-metastatic breast cancer. In a transcriptome sequencing study, Han et al. (2021) explored expression profiles of the non-coding RNAs in primary LUAD and LUAD-BM and identified 706 differentially expressed circRNAs. In non-small-cell lung cancer (NSCLC) bone metastasis cell line, circ_0060937 knockdown inhibited cell proliferation and invasion (Zhang et al., 2020). To the best of our knowledge, this is the first study to investigate the circRNA expression profile of LUAD-BM using circRNA sequencing, and above results have displayed that circ_0096442 was identified as the most significant biomarker in LUAD-BM. Basing on the literature review, we found that circ_0096442 was a novel biomolecular in LUAD-BM, and the molecular mechanisms of it in regulating tumor metastasis and growth were also not reported ever before. In other words, miR-326 was firstly reported was the down-stream target of circ_0096442 in regulating tumor metastasis and growth. In fact, previous studies have reported that miR-326 was a tumor inhibitor in lung cancer, which inhibited lung cancer metastasis, growth and chemotherapy resistance by HOTAIR/miR-326/phoox2a pathway, miR-326/ZEB1 pathway and RHOT1/miR-326/FOXM1 pathway (Wang et al., 2016; Liu et al., 2021; Zhang et al., 2021).

BM is a complex, multistep process requiring tumor cells to detach from their original sites, intravasate into blood vessels, and subsequently survive and proliferate in the bone (Cheung and Ewald, 2016). Therefore, peripheral blood samples, which may better reflect the BM processes, were selected for investigation. Furthermore, for biomarker discovery, this procedure is more easily accepted by patients as it is relatively economical and minimally invasive. Circulating non-coding RNAs in the peripheral blood system were reported to be useful biomarkers for the detection of diseases (Zhang et al., 2022a; Zhang et al., 2022b; Liu et al., 2022). However, it is possible that the changes in the expression profile of the disease reflect shifts in cell populations (Keller et al., 2011). Although patients in this study were matched for age, sex, body mass index, and region, it should be considered as a potential confounding factor that peripheral blood can be affected by various factors. Moreover, this study has a limited sample size; thus, a larger clinical sample size is needed to confirm our results in future studies. Furthermore, the mechanisms and functions of the circRNAs predicted by bioinformatic analyses should be further elucidated by more rigorous molecular biology experiments.

Conclusion

This study is the first to systematically characterize and annotate circRNA expression in patients with LUAD-BM. Our findings provide a novel and important circRNA expression profile of LUAD-BM and suggest circ_0096442 as a possible biomarker for LUAD-BM.