Transcriptome Sequencing reveals the expressed profiles of mRNA and ncRNAs and regulate network via ceRNA mediated molecular mechanism of lung adenocarcinoma bone metastasis in Xuanwei.

BACKGROUND: The most ordinary subtype of lung cancer is lung adenocarcinoma (LuAC), which is characterized by strong metastatic ability. And LuAC rates in Xuanwei leads to the poor prognosis and high death rate. In this study, we systematically explored the molecular mechanism of LuAC bone metastasis in Xuanwei by transcriptome sequencing.
METHODS: RNA Sequencing was conducted to explore the noncoding RNAs (ncRNAs) expression profiles in primary LuAC and LuAC bone metastasis. We identified differentially expressed mRNAs (DEmRNAs), miRNAs (DEmiRNAs), lncRNAs (DElncRNAs) and circRNAs (DEcircRNAs). Bioinformatics analyses the possible relationships and functions of the LuAC bone metastasis-related competing endogenous RNA (ceRNA). And qRT-PCR was performed to evaluate the expression of these differently expressed genes in serum.
RESULTS: A total of 2,141 DEmRNAs, 43 DEmiRNAs, 136 DElncRNAs and 706 DEcircRNAs were identified in the Xuanwei patients with primary LuAC vs. LuAC bone metastasis, respectively. The circRNA/miRNA/mRNA and lncRNA/miRNA/mRNA networks of LuAC in Xuanwei with bone metastasis were built, and the gene expression mechanisms regulated by ncRNAs were unveiled via the ceRNA regulatory networks. We observe that lncRNA (ADAMTS9-AS2, TEX41, DLEU2, LINC00152)-miR-223-3p-SCARB1 and hsa_circ_0000053-miR-196a-5p/miR-196b-5p-HOXA5 ceRNA networks might play an important role in bone metastasis of Xuanwei LuAC.
CONCLUSIONS: We comprehensively identified ceRNA regulatory networks of LuAC in Xuanwei with bone metastasis as well as revealed the contribution of different ncRNAs expression profiles. Our data demonstrate the association between mRNAs and ncRNAs in the metastasis mechanism of LuAC in Xuanwei with bone metastasis. 2021 Translational Cancer Research. All rights reserved.

Introduction

Lung cancer is one of the most malignant tumors that poses a threat to the health and life of the population (1). Its morbidity and mortality are the fastest growing, which is the main cause of the global cancer death rate (2). Among them, non-small cell lung cancer (NSCLC) account for 85% of lung cancer (3). For patients with early-stage NSCLC who have no evidence of mediastinal invasion or distant metastasis before surgery, surgical resection of the lung mass and the collection of mediastinal lymph nodes is the best treatment. However, the recurrence rate after surgery is very high (4,5). NSCLC can be divided into three subtypes of lung adenocarcinoma (LuAC), squamous cell carcinoma and large cell carcinoma. The most ordinary subtype of lung cancer is LuAC, which is characterized by strong metastatic ability and high mortality (6). The LuAC rates in Xuanwei, China, are among the highest around the world, especially in nonsmoking women (7). Various of risk factors for LuAC rates are reported, such as smoking, diet, atmospheric pollution, genetic factors, etc. Among them, smoky coal for cooking and heating is considered as a major risk factors in its pathogenesis (7). But because of all kinds of reason, the death rate hasn’t to be controlled effectively (8). And the patients are already at an advanced stage when they are diagnosed due to the scarcity of effective biomarkers for early detection (9). Several newer targeted therapies may shed new light on potential therapeutic breakthroughs due to surgical limitations, chemotherapy resistance, radiotherapy side effects and other factors, such as PD-1 noncoding RNAs (ncRNAs) and so on (10-13). Simultaneously, bone metastases are common in patients with (LuAC and are commonly osteolytic in nature (14). The survival time of lung cancer patients after bone metastasis is 6–10 months, while the 1-year survival rate is only 40–50% after the clinical treatment (15,16). The study of associated driven genes in lung cancer bone metastasis will provide new insight into the mechanism, potential therapeutic targets, promote the prognosis and survival in clinic (17,18). Although bone metastasis has important clinical significance in lung cancer, its pathological mechanism is still indistinct. Therefore, taking Xuanwei LuAC as the core research object, in-depth study of the pathogenesis of bone metastasis and systematic screening of early diagnosis and potential treatment targets will surely provide vital basis for clinical and scientific research.

Noncoding RNAs (ncRNAs) are central components of the human transcriptome that functionally regulate protein expression (19). The ncRNAs include long ncRNAs (lncRNAs, >200 bp, up to 100 kb), microRNAs (miRNAs <200 bp) and circular RNAs (circRNA, closed loop structure) (20,21). Competing endogenous RNAs (ceRNAs) hypothesis reveals a new mechanism for interaction between RNAs. MiRNAs can cause gene silencing by binding to mRNA. However, ceRNA can regulate gene expression by competitively binding to miRNAs. LncRNAs or circRNAs competitively binds to miRNA and acts as a sponge to regulate gene expression. CeRNA can bind miRNAs through miRNA response elements (MREs), affecting miRNA-induced gene silencing, thus leading to elevated levels of miRNA target genes (22,23). It has been reported that they play key roles in varieties of biological processes such as angiogenesis, cell proliferation, differentiation, apoptosis, autophagy, and immune response (24-28). The long non-coding RNA activated by TGF-β, lncRNA-ATB, acts as a ceRNA and competes with ZEB1 and ZEB2 to bind to miR-200 to drive the invasion and metastasis cascade of hepatocellular carcinoma (29). CircRNA_100290 can play a part in oral cancer by functioning as a sponge of miR-29 (30). ZEB1 initiates miR-181b-regulated ceRNA network to drive LuAC metastasis (31). LncRNA metastasis-associated LuAC transcript1MALAT1, which is a ceRNA that regulates autophagy-associated 7 (ATG7) gene expression via the sponge miR142-3p (32). As far as we all know, the ceRNA regulatory networks and its basic mechanism in LuAC bone metastasis have not been completely inquired.

High-through put transcriptome sequencing has been widely used to analyse particularly gene expression modes at different developmental phases. In present study, we investigated the mRNA, miRNA, lncRNA and circRNA expression in primary LuAC and LuAC bone metastasis in Xuanwei, China using RNA-sequencing. We structured the ceRNA network to detail the interactions or potential crosstalk among the key miRNA, mRNA, lncRNA and circRNA in primary LuAC and LuAC bone metastasis. Furthermore, through the integrated prediction and research of the ceRNA mechanism, we hope to screen a group of ncRNAs as biomarkers for the diagnosis of LuAC bone metastasis and come out a novel diagnostic and therapeutic strategy against LuAC bone metastasis in Xuanwei.

We present the following article in accordance with the MDAR checklist (available at http://dx.doi.org/10.21037/tcr-20-2376).

Methods

Study population

We used a retrospective cohort study design. In our study, the Xuanwei patients come from the Third Affiliated Hospital of Kunming Medical University. The detailed clinical data of patients are shown in . Fresh tumor tissue samples were procured from 3 cases of primary LuAC (P1-P3) and 3 cases of LuAC bone metastasis (M1-M3) after surgical resection. Tissue samples were frozen in liquid nitrogen prior to RNA extraction and stored at −80 °C. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The study was approved by ethics board of the Third Affiliated Hospital of Kunming Medical University (Yunnan Cancer Hospital) (No. KMUYXLL0046) and informed consent was taken from all the patients.

Table 1

Patient clinical information

Index	P1	P2	P3	M1	M2	M3
Age, years	50	73	64	43	74	65
Gender	Male	Female	Female	Female	Male	Female
Site of metastasis	–	–	–	Bone	Bone	Bone

RNA separation extraction and quality monitoring

RNA extraction of all tissue samples was performed using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manual instruction. The RNA quality monitoring was performed using an Agilent 2100 Bioanalyzer (Agilent Technologies) and NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA).

MiRNA library construction, sequencing, and data processing

The library was constructed that used the Truseq RNA sample Prep Kit (Illumina, Inc., San Diego, CA, USA). The adapter is connected to the 3’ end of the RNA and then to the 5’ adapters. Subsequently, the RNA was fragmented into short base pairs when added the lysis buffer. Synthesis of single-stranded cDNA was performed by RNA reverse transcription. Amplification of miRNA high-throughput sequencing libraries was performed by PCR. Sequencing was performed by a lumina HiSeq4000 (Illumina, San Diego, CA, USA) platform.

The internal pipeline processes raw data (original reads), including adapter pruning, read alignment, and read counts. By using the popular alignment tool Bowtie (http://bowtie-bio.sourceforge.net/index.shtml) to map the trimmed reads (also known as clean reads) to the human reference genome GRCh38 (33). Then, using the software miRDeep2 (https://www.mdcberlin.de/8551903/en/) to quantify the expression of miRNA (34). miRBase downloaded mature miRNA and miRNA precursors. Moreover, DEGseq package in R was used to identify DEmiRNAs between samples. The cut-off criteria are P value <0.05 and |log2 (Fold_change)| >2.

MiRNA and lncRNA sequencing and data processing

The library construction standards were as follows: (I) concentration of RNA ≥200 ng/mL; (II) RNA with amount >5 µg; (III) RIN >7; (IV) 1.8< OD260/280 <2.2. Fragmentation of rRNA-deleted RNA after removal of ribosomal RNA, using the Truseq RNA sample Prep Kit (Illumina, Inc., San Diego, CA, USA) was got to structure the cDNA library. Sequencing of libraries were performed by an Illumina HiSeq x-ten (Illumina Inc., San Diego, CA, USA) platform.

All the raw data obtained from high-throughput RNA-sequencing was translated into raw FASTQ sequence data by using Base Calling. To acquire the clean reads from the RNA sequencing results of lncRNA and mRNA, we used cutadapt (http://cufflinks.cbcb.umd.edu/) and Fastx-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/) to remove lower quality sequences, including adapter sequences, sequences with mass fraction <20 and sequences with N base rates of raw reads >10% and sequence less than 18 bp. The mapped reads were quantified with cuffquant, via using the Cufflinks (http://cufflinks.cbcb.umd.edu/) package of Cuffdiff program identify the DEmRNAs and DElncRNAs between samples. The P value <0.05, with |log2 (Fold-change)| >2 as the cutoff criteria.

CircRNA library construction, sequencing, and data processing

Removing the ribosomal RNA from total RNA, Digestion of linear RNA with RNaseR was carried out to isolate circRNAs. Then, the cDNA was constructed by further adding the reverse transcriptase to the RNA sample with RNase R. The 3’end of the cDNA is ligated with the base “A” and ligated to the linker, then the fragment size is selected, and finally PCR amplification is performed. Sequencing of libraries were performed by an IlluminaHiSeq x-ten (Illumina Inc., San Diego, CA, USA) platform.

Raw reads of FASTQ format were then processed through in-house Perl scripts and removed low quality reads using software SOAPnuke. Align clean RNA sequencing reads with reference genome, using CIRI (v2.0.5) and find circ (v1.2) predict circRNA. After combining the results of the two softwares, quantitative and differential expression analysis of circRNA was performed. We analyzed the circRNA abundance by using the normalized reads count. DEcircRNAs were evaluated by DESeq2 (R software package) (2). In the differential expression analysis, P value <0.05, with |log2 (Fold_-change)| >2 was as the cutoff criteria.

Functional analysis of DEmRNAs and DEncRNAs

To explore the functions of DEmRNAs and DEncRNAs, using the annotation of gene ontology (GO) and analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathways. However, for the differential expression genes, we use GO analysis (http://www.geneontology.org) to create hierarchical categories based on molecular function (MF), cellular composition (CC), and biological process (BP) characteristics to elucidate genes. In addition, KEGG (http://www.genome.jp/kegg/) was used to analyse the important pathway in DEmRNAs and DEncRNAs. Using Metascape (http://metascape.org/gp/index.html) bioinformatics resources to perform GO functional annotation and KEGG pathway analysis of differentially expressed genes. P<0.001 was selected as the cutoff standards of significance.

Analysis of the ceRNA regulatory network

The ceRNA analysis was performed on mRNAs, lncRNAs and circRNAs with significant correlations in expression levels, searching for potential MREs in the lncRNA, circRNA and mRNA sequences. The lncRNA/circRNA-miRNA-mRNA interaction was predicted by considering the overlap of lncRNA/circRNA with the same miRNA seed sequence binding site in the mRNA sequence. The miRNA binding sites were predicted by miRcode (http://www.mircode.org/) and TargetScan (http://www.targetscan.org/) predicts the miRNA-mRNA interactions.

QRT-PCR verify differential gene expression

Total RNA was extracted using TRIzol reagent (Invitrogen, USA) according to the instructions. Reverse transcription into cDNA using reverse Transcriptase Kit (TaKaRa, Dalian, China) using mRNA as a template. The miRNA was reverse transcribed using the miRcute miRNA First-strand cDNA Synthesis Kit (TaKaRa, Dalian, China). Designed the forward and reverse primers, qRT-PCR detection was performed on a BIO-RAD IQ5 RT-PCR detection system (Bio-Rad Laboratories Inc., Germany). Primer sequences of qRT-PCR reactions were shown in Table S1. Calculating the relative expression of DEmRNAs and DEncRNAs used the 2−ΔΔCt method. GAPDH was used as control. The lncRNAs, circRNAs and miRNAs expression levels were normalized against the snU6.

Statistical analysis

All statistical analyses were performed by GraphPad Prism 7.0 (GraphPad Software, La Jolla, CA, USA). Two independent sample t test was used to compare the mean of measurement data between groups. And the data was represented as means with standard deviation (SD). α=0.05 was used as the statistical test level in this study.

Results

Identification of DEmRNAs and DEncRNAs

Significant analysis of RNA sequencing was performed using Cuffdiff software and analyse differential expressed mRNAs and ncRNAs the standards of P<0.05. The clustering heatmap of DEmRNAs, DEmiRNAs, DElncRNAs and DEcircRNAs were showed in separately. In all of 2,141 DEmRNAs (1046 up-regulated and 1095 down-regulated DEmRNAs) were confirmed in the Xuanwei patients with primary LuAC vs. LuAC bone metastasis, respectively (). In addition, there are 43 DEmiRNAs, including 23 upregulated and 20 downregulated miRNAs and 136 DElncRNAs, including 76 upregulated and 60 downregulated lncRNAs; 706 DEcircRNAs, including 365 upregulated and 351 downregulated circRNAs in the Xuanwei patients with primary LuAC vs. LuAC bone metastasis, respectively ().

Figure 1

RNA-seq reveals distinct expression pattern of mRNAs, miRNAs, lncRNAs and circRNAs in primary lung adenocarcinoma patients and lung adenocarcinoma bone metastasis patients. (A,B,C,D) Unsupervised clustering analysis showing expression profiles of DEmRNAs, DEmiRNAs, DElncRNAs and DEcircRNAs in the patients with primary lung adenocarcinoma vs. lung adenocarcinoma bone metastasis.

Functional annotation of DEmRNAs and DEncRNAs

To investigate the function of ncRNAs in the Xuanwei patients with LuAC bone metastasis, we used GO annotation and KEGG pathway analysis to identify the functions of DEmRNAs and DEncRNAs between LuAC bone metastasis and primary LuAC patients in Xuanwei. LncRNAs might regulate the nearby protein-coding genes expression. In the three major GO categories of biological process, molecular function and cellular component, there are overly representative terms in GO analysis. The top 15 outcomes of gene-annotated DEmRNAs and DEncRNAs were analyzed for biological processes (BP), cellular components (CC) and molecular function (MF). Sort by P value and select the prediction term with P value <0.05. Based on the prediction of the dissemination of the target gene in GO, the gene functions of DEncRNAs were clarified. Each GO item enrichment statistically analyzes the number of genes and is represented in the form of histograms.

The most enriched BPs were cell division and cell adhesion that were showed by the GO analysis of DEmRNAs. The proteinaceous extracellular matrix, plasma membrane, chromosome and spindle microtubule were the meaningfully enriched in the CCs. The single-stranded DNA-dependent ATPase activity, microtubule motor activity and extracellular matrix structural constituent were the importantly enriched in the MF ().

Figure 2

GO analysis of the DEmRNAs and DEncRNA between primary lung adenocarcinoma patients and lung adenocarcinoma bone metastasis patients. (A) GO analysis of DEmRNAs in the lung adenocarcinoma bone metastasis patients; (B) GO analysis of DEmiRNAs in the lung adenocarcinoma bone metastasis patients; (C) GO analysis for nearby protein-coding genes of DElncRNAs in the lung adenocarcinoma bone metastasis patients; (D) GO analysis of DEcircRNAs in the lung adenocarcinoma bone metastasis patients. Gene numbers in the GO term were presented in histography. The −log10 (P value) yields an enrichment score representing the significance of GO term enrichment among differentially expressed mRNAs and ncRNA.

GO analysis for the targeted gene of DEmiRNAs revealed that the most expressively enriched BPs were extracellular structure organization, negative regulation of thymocyte apoptotic process, skeletal system development and ceramide biosynthetic process. The cell leading edge, receptor complex, apical part of cell and axon were the most noteworthy enriched in the CCs. The ubiquitin-protein transferase activity, metalloendopeptidase activity, single-stranded DNA binding and lipoprotein particle receptor activity were the most meaningfully enriched in the MF ().

The GO analysis for nearby protein-coding genes of DElncRNAs revealed that the most significantly enriched BPs were positive regulation of inflammatory response, positive regulation of GTPase activity and osteoblast differentiation. The microtubule cytoskeleton, primary cilium and ciliary basal body were the most significantly enriched in the CCs. The chemokine activity3’,5’-cyclic-nucleotide phosphodiesterase activity and transforming growth factor beta receptor binding were the importantly enriched in the MF ().

The GO analysis for the source gene of DEcircRNAs revealed that the most noteworthy enriched BPs were microtubule-based process, regulation of GTPase activity, regulation of mitotic cell cycle phase transition and cell morphogenesis involved in differentiation revealed that the most expressively enriched in the BPs. The microtubule organizing center, centrosome, spindle, microtubule, cell leading edge and transferase complex were the importantly enriched in the CCs. The phosphotransferase activity, alcohol group as acceptor, tubulin binding, nucleoside-triphosphatase regulator activity and kinase binding were the most noticeable enriched in the MF ().

To acquire further comprehensions of the different biological functions of DEmRNAs and DEncRNAs, we performed KEGG pathway analysis. The top 20 pathways in the DEmRNAs are shown (, Table S2). Intriguingly, some of identified KEGG pathways comprised cell adhesion molecules (CAMs), pathways in cancer, ECM-receptor interaction and p53 signaling pathway. The total 16 pathways in the target mRNAs were based on the DEmiRNAs were shown (, Table S3). The most meaningfully associated pathways were PI3K-Akt signaling pathway, cell cycle, TNF signaling pathway and AMPK signaling pathway. The total 11 pathways in the nearby mRNAs of DE lncRNAs were shown (, Table S4). The most notably pathways were FoxO signaling pathway, endocytosis, TGF-beta signaling pathway and MAPK signaling pathway. The total 12 pathways in the source gene of DEcircRNAs were shown (, Table S5). The importantly pathways were adherens junction, Rap1 signaling pathway and regulation of actin cytoskeleton. These pathways are suggested to be involved in the metastatic mechanisms of LuAC bone metastasis.

Figure 3

KEGG pathways for DEmRNAs and DEncRNA between primary lung adenocarcinoma patients and lung adenocarcinoma bone metastasis patients. (A) KEGG of DEmRNAs in the lung adenocarcinoma bone metastasis patients; (B) KEGG of DEmiRNAs in the lung adenocarcinoma bone metastasis patients; (C) KEGG of DElncRNAs in the lung adenocarcinoma bone metastasis patients; (D) KEGG of DEcircRNAs in the lung adenocarcinoma bone metastasis patients. Identification of linear transcripts in different biological pathways by KEGG pathway analysis. The −log10 (P value) produces an enrichment score indicating the importance of pathway correlation.

Construction of ceRNA regulatory network

The study found that RNA regulates MREs through a mechanism called the ceRNA hypothesis. MREs involved in ceRNA networks can regulate mRNA expression and have also been discovered (35). The role of lncRNAs and circRNAs as miRNA sponges can inhibit the regulation of target genes on miRNAs, thereby indirectly regulating gene expression and functioning. Accordingly, lncRNA/circRNA-miRNA-gene and lncRNA/circRNA (as a sponge), miRNA (as the center) and mRNA (as the target) relationship pairs were established. From the high-throughput sequencing data of tissue samples, we integrated the expression details and regulatory relationships of the mRNAs, miRNAs, lncRNAs and circRNAs to construct a ceRNA network (). Our results showed that the 67 mRNAs, 27 lncRNAs and 12 miRNAs were involved in the lncRNA-miRNA-mRNA ceRNA network (). However, the ceRNA networks were found on circRNA-miRNA-mRNA interactions include 40 circRNAs, 17 miRNAs and 71 mRNAs (). These ceRNA regulatory networks contain circRNAs/ lncRNAs, miRNAs and mRNAs, which may play a key role in the LuAC in Xuanwei with bone metastasis.

Figure 4

The lung adenocarcinoma bone metastasis-related ceRNA network. (A) Regulatory network analysis of lncRNAs-miRNAs-mRNAs in lung adenocarcinoma bone metastasis patients. Global view of the ceRNA network. This network involves 27 DElncRNAs, 12 DEmiRNAs and 67 DEmRNAs. (B) Regulatory network analysis of circRNAs-miRNAs-mRNAs in lung adenocarcinoma bone metastasis patients. Global view of the ceRNA network. This network involves 40 DEcircRNAs, 17 DEmiRNAs and 71 DEmRNAs. Elliptic circular nodes represent mRNAs, rectangular nodes represent miRNAs, triangular nodes stand for lncRNAs and rhombus nodes stand for circRNAs. Red means up regulate, green means down regulate.

QRT-PCR confirmed the ncRNAs and mRNA expression

Through bioinformatics analysis and literature review, we finally selected 12 significantly different candidate genes to verify by qRT-PCR in Xuanwei LuAC samples. We verified the expression of 2 DEcircRNAs (circ_0000053 and circ_00001495), 2 DElncRNAs (ADAMTS9-AS2 and TEX41), 4 DEmiRNAs (miR-196b-5p, miR-223-3p, miR-135a-5p and miR-9) and 4 DEmRNAs (MEIS1, HOXA5, SCARB1 and RUNX2) (). Compared with Xuanwei patients with LuAC and normal, the expression of circ_0000053, circ_00001495, TEX41, miR-223-3p, miR-9, RUNX2, SCARB1 and ADAMTS9-AS2, miR-135a-5p, MEIS1, HOXA5 were detected significantly upregulated and down regulated in Xuanwei patients with LuAC bone metastasis (). Meanwhile, in this study, the expression of these eleven genes was united with the RNA-seq results (Table S6).

Figure 5

qRT-PCR confirmed the ncRNAs and mRNA expression. qRT-PCR validate the expression level of key ncRNAs and mRNAs in Normal (N) (n=15) vs. lung adenocarcinoma patients (P) (n=15) vs. lung adenocarcinoma bone metastasis patients (M) (n=14). *, P<0.05; **, P<0.01; ***, P<0.001.

However, in our study, the expression of miR-196b-5p was inconsistent with that in RNA-seq results (, Table S6).

Discussion

Recent study shows that, Xuanwei has the highest morbidity and mortality of LuAC in China. And the rates are also rather high among non-smoking women (36). Bone is a common transfer site in cancer patients with LuAC. Around 30–40% of patients with advanced lung adenocarcinoma cancer occur bone metastases during their onset, this will lead a remarkable negative influence on morbidity and survival (37). The presence of bone metastases suggests a poor prognosis, with a 1-year survival rate of 40% to 50% and a median survival time of only 6 to 10 months (38,39). In recent years, non-coding RNAs have been found in multiple cancers and are connected with comprehensive biological regulatory functions (40). Circulating non-coding RNAs in blood (serum, plasma) as important potential of biomarkers in liquid biopsies for diagnosis and follow-up of patients with all kinds of tumors (41,42).This study used a next-generation sequencing technique to quantify and comprehensively analyze the coding and non-coding transcriptomes of primary LuAC and LuAC bone metastasis in Xuanwei. Our analyses proved remarkable differences in the expression patterns of mRNA, lncRNA, circRNA and miRNA in primary LuAC and LuAC bone metastasis in Xuanwei, and dynamic changes in DEmRNA, DElncRNA, DEmiRNA and DEcircRNAs. We first found that the expression patterns of mRNAs and ncRNAs and are more appropriate for the identification of primary LuAC and LuAC bone metastasis samples in Xuanwei.

Overall, our result has proved that mRNAs, miRNAs, circRNAs and lncRNAs different populations participated in the pathogenesis of LuAC bone metastasis in Xuanwei. Unsupervised hierarchical clustering analysis the expression profiles of DElncRNAs, DEcircRNAs, DEmiRNAs and DEmRNAs, it was found that mRNA and ncRNAs expression profiles can largely distinguish between the primary LuAC and LuAC bone metastasis in Xuanwei. Using GO and KEGG pathway analyze the potential functions of DEmRNAs and DEncRNAs in the Xuanwei LuAC bone metastasis. In addition, according to these findings, we suggested that ceRNAs play an important role in Xuanwei LuAC bone metastasis and the ceRNA regulatory networks were established. We found that most of the DEmiRNAs may be associated with Xuanwei LuAC metastasis mechanism. However, due to the lack of research, majority of the DElncRNAs/DEcircRNAs were with unknown function.

A total of 2,141 DEmRNAs, 43 DEmiRNAs, 136 DElncRNAs and 706 DEcircRNAs were identified in both primary LuAC and LuAC bone metastasis in Xuanwei. DEmRNAs and DEncRNAs were significantly enriched in ECM-receptor interaction, PI3K/Akt signaling pathway and TGF-β signaling pathway (43). Previous studies also indicated that the dysregulation of ncRNAs and mRNA may contribute to Xuanwei LuAC metastasis. MiR-196b as a tumor suppressor target down-regulate Runx2 and regulate the PI3K/AKT/GSK3β, JNK and Smad signal pathways to suppress lung cancer cells growth and metastasis (44). MiR-130a-5p directly targeted RUNX2 regulates tumor invasive and metastatic potential in NSCLC (45). MiR-449a as a tumor suppressor in NSCLC and its function is through targeting HMGB1-mediated NF-κB signaling pathway, induced suppression of growth, migration and invasion (46). MiR-370 may bind to the 3’UTR of EGFR to restrain EGFR expression, inhibition of growth, angiogenesis and metastasis of NSCLC by down-regulating ERK1/2 and AKT signaling pathways (47). MiR-196a as a key mediator, through regulating FoxO1, CDKN1B and HOXA9, mediated the proliferative, pro-migratory and tumorigenic activity of lung cancer cells (48). MiR-223-3p directly inhibited the expression ofNLRP3 and IкB kinaseα activity and these two important mediators are thought to be associated with I/R-induced cell necroptosis and inflammation in the treatment of ischemic heart disease (49). MiR-223-3p regulated RIPK3 expression by directly targeting the 3’UTR of RIPK3, up-regulated 3-MCPD-dipalmitate-induced AKI and activated RIPK1/RIPK3/MLKL associated necrotic pathway in acute kidney injury (50). MiR-223-3p targets FOXO3 and promotes osteogenic differentiation of bone marrow mesenchymal stem cells by enhancing autophagy (42). MiR-9 regulates ferroptosis in melanoma cells by targeting GOT1, indicating the crucial role of miRNA in ferroptosis (51). In our current study, we showed that miR-196a-5p, miR-9-3p and miR-223-3p were upregulated in Xuanwei LuAC bone metastasis. Therefore, these miRNAs may act as important roles in Xuanwei LuAC metastasis mechanism, ferroptosis and necroptosis. MiR-196a can enhances cell migration, invasion and EMT by regulating the expression of HOXA5 in osteosarcoma (52). However, the expression of miR-196b-5p was inconsistent with that in RNA-seq results in our study (45). The reason may be related to the expansion of the sample size in the qRT-PCR and the selection of serum as the test sample. However, the function of miR-196b-5p in Xuanwei LuAC bone metastasis and its mechanism need further experimental verification.

Most of the DElncRNAs/DEcircRNAs were few studied. High expression of lncRNA ADAMTS9-AS2 in lung cancer cells obviously could reduce proliferation ability and inhibit migration, as well as promoting their apoptosis rate (43). The lncRNATEX41 was the variants of AS and affect the pathogenesis of aortic valve stenosis (53). TEX41-miR-340-COMMD6 ceRNA network to play an important role in head and neck squamous cell carcinoma (54). The non-coding RNA Linc00152 is a newly discovered lncRNA that has been reported to participate in tumor initiation and progression (55,56). LncRNA DLEU2 influences the proliferation, migration, and invasion of laryngeal cancer (57). In the breast cancer, hsa_circ_0002113 expression upregulated regulated breast cancer tumorigenesis and progression (58). However, lncRNA and circRNA of LuAC bone metastasis have little been reported until now.

Pseudogenes, circRNAs, lncRNAs and mRNAs can be used as ceRNAs because of RNA containing miRNA response elements. As natural miRNA sponges, they inhibit miRNA function by sharing miRNA response elements and coordinate with each other (59). In the process of LuAC bone metastasis, the underlying mechanism by which ncRNAs interact with mRNAs is unclear. Therefore, lncRNA/circRNA-miRNA-mRNA of LuAC bone metastasis was constructed that basing on the data analysis from next generation sequencing. These breakthrough discoveries will increase our understanding of the function of ncRNAs in the metastasis mechanism of LuAC bone metastasis. CeRNA network prediction and bioinformatics analysis are able to fully understand the ncRNA function, which may involve in the initiation and progression of the disease (60).

In our ceRNA network, we also found lncRNA ADAMTS9-AS2, TEX41, DLEU2 and LINC00152 were ceRNAs of miR-223-3p targeting SCARB1 (alias SR-B1). The SR-B1is a high-affinity HDL receptor that has been shown to mediate HDL-dependent lipid transport signaling in a variety of different cell types (61). The SCARB1 targeting of HPPS (HDL-mimicking peptide-phospholipid nanoscaffold (HPPS) nanoparticle composed of the cholesteryl oleate) plays an important role inefficient delivery in lysosomes (62). The research has identified a role of SR-B1 in regulating the protective effects of HDL against necrosis in myocardial necrosis (63) The circRNA-associated-ceRNA networks involved in the LuAC bone metastasis was selected. Hsa_circ_0000053 was ceRNAs of miR-196a-5p/miR-196b-5p targeting HOXA5. HOXA5 belongs to the Hox gene family encoding transcription factors known for its role in the formation of skeletal patterns (64). The dimer between HOX and PBX inhibits the necrotic pathway, possibly through inhibition of RIP kinase. HOX/PBX dimer also inhibit the expression of p21, which might in turn block apoptosis in acute myeloid leukemia cells (65). In summary, we speculated that lncRNA, circRNA and miRNA molecules may be pivotal regulators of varying Xuanwei LuAC course, such as metastasis and necroptosis.

Conclusions

In conclusion, our current research confirmed the extensive expression profile of ncRNAs in Xuanwei LuAC bone metastasis and provided a breakthrough data incorporation analysis of ceRNAs in Xuanwei LuAC with bone metastasis. A catalog of predicted ncRNAs in Xuanwei LuAC with bone metastasis was also defined. These rich data suggested that ceRNAs may regulate the expression of genes encoding related proteins and play a crucial role in Xuanwei LuAC with bone metastasis. Taken together, the findings further extended our knowledge of ceRNA biology and contributed to the understanding of their regulation roles in Xuanwei LuAC with metastasis mechanism and necroptosis. Regardless, our results offer worthy information for ncRNAs studies and clinical diagnosis, treatment, and prevention of Xuanwei LuAC.