This article is mentioned in:

Introduction

Cellular components of bone marrow have important roles in pre-metastatic niche (PMN) formation (1). In lung cancer, distant metastases are common and this type of cancer usually spreads to the bone (39%), liver (35%) and central nervous system (47%) (2). Patients with lung cancer and metastasis have a poor prognosis with a shortened median survival time following diagnosis (3). In order to metastasize, tumor cells need an organ with a suitable environment for their growth and proliferation, which is defined as the metastatic niche (4).

Cancer cells initiate and establish the environmental surroundings required for future metastasis through various mechanisms, including cancer cell intravasation, immune evasion and arrival at designated site, extravasation, colonization and tumor growth (5). The bone marrow microenvironment has been described as fertile ground for dormant and proliferating tumor cells. For example, bone marrow and tumor cells can modify the activity of osteoclasts (6), and pro-tumorigenic cells, including mesenchymal stem cells, have been reported to serve a crucial role in promoting osteolytic bone metastasis and tumor cell proliferation in the tumor microenvironment (7). Additional tumor-derived factors have been reported to promote tumor progression. These factors can stimulate the differentiation of immature myeloid cells into strong immune response suppressors and therefore inhibit the activation of antitumor T cells (8).

Numerous factors, including tumor-derived secreted factors and extracellular vesicles, are involved in PMN establishment (5). In addition, other cell types, including bone marrow-derived cells (BMDCs) such as mesenchymal stem cells and regulatory T cells, are directed to the secondary organs. Once these cells have reached the PMN, they modify its local microenvironment through inflammatory cytokines, growth factors and proangiogenic molecules to facilitate tumor cell colonization and proliferation, and therefore promote tumor metastasis (4,5,9).

Notably, a PMN is established through the combination of various tumor-derived factors, tumor-mobilized BMDCs and the local environment (5,9). However, the role of BMDCs in PMN formation is not yet fully understood. In the present study, it was hypothesized that lung cancer cells can remotely modify BMDCs, which could therefore become actively involved in PMN establishment in target organs. The present study aimed to investigate the role of BMDCs in lung cancer metastasis from a macroscopic perspective (Fig. 1). To do so, bone marrow tissue samples were examined by next generation sequencing (NGS).

Figure 1. Study flowchart. To simulate a pre-metastatic niche in lung cancer, LLC cells were injected into the tail veins of C57BL/6J mice for 10 days. The bone marrow of the mice was harvested for further NGS analysis. NGS data were presented in different ways, including volcano plot, and underwent GO and KEGG analyses, and meta-analysis. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; LLC, Lewis lung carcinoma; LLC-BMDCs, LLC-bone marrow-derived cells; NGS, next generation sequencing.

Materials and methods

Cell culture

The LL/2 mouse Lewis lung carcinoma (LLC) cell line [LLC1; American Type Culture Collection (ATCC)® CRL-1642™] was purchased from the ATCC (Manassas, VA, USA) and cultured in Dulbecco's modified Eagle's medium (DMEM) containing 10% fetal bovine serum (Gibco; Thermo Fisher Scientific, Inc., Waltham, MA, USA) in an incubator containing 5% CO2 at 37°C.

BMDCs isolation

Male C57BL/6J mice (age, 8 weeks; weight, 20–25 g) were obtained from the National Laboratory Animal Center (Taipei, Taiwan). These mice were housed at a constant temperature (21±1°C) and humidity (55–65%) under a 12-h light/dark cycle. Mice had free access to food and water. Animal experiments were performed according to the Institutional Animal Care and Use Committee, with the approval of the Animal Care and Use Committee of the School of Kaohsiung Medical University (Kaohsiung, Taiwan). The study group contained six mice that were injected with 1×106 LLC cells, whereas the control group comprised six mice that received sham injection with DMEM. All injections were administered into the tail vein. Mice were sacrificed after 10 days. The LLC-BMDCs and normal-BMDCs were obtained after femoral bone removal, followed by flushing the bone marrow out from the femurs bilaterally with PBS. The collected BMDCs were suspended in DMEM (10). The BMDCs from mice contained in each group were pooled prior to further experiments.

RNA sequencing

Total RNA was isolated from BMDCs using TRIzol® reagent (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol. RNA purity was measured by assessing optical density (OD)260nm/OD280nm absorbance ratio (1.96 for normal-BMDCs and 2.03 for LLC-BMDCs) with an ND-1000 spectrophotometer (NanoDrop Technologies; Thermo Fisher Scientific, Inc., Wilmington, DE, USA). RNA concentration was determined according to the RNA integrity number (RIN; RIN, 10 for normal-BMDCs and LLC-BMDCs) with an Agilent Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA).

Library preparation and deep sequencing were performed using Illumina Solexa system (Illumina, Inc., San Diego, CA, USA) according to the manufacturer's protocol, as previously described (11,12). For small RNA sequencing, total RNA was reverse transcribed using TruSeq Small RNA Sample Prep kit (Illumina, Inc.) following the manufacturer's instructions, and cDNA (18–40-nucleotide RNA fragments; 140–155 nucleotides in length with both adapters) was sequenced on the Illumina system (75 single-end cycles). After trimming or removing low-quality data using Trimmomatic software version 0.36 (13), the qualified reads were analyzed using miRDeep2 software (miRBase 21) (14) and the human genome from the University of California Santa Cruz database (https://genome.ucsc.edu/). Micro (mi)RNAs with low levels [<1 normalized read per million (RPM)] in normal- and LLC-BMDCs were excluded. For transcriptome sequencing, the library constructed with SureSelect Strand Specific RNA Library Preparation kit (Agilent Technologies, Inc.) was sequenced using a TruSeq SBS kit on the Solexa platform (Illumina NextSeq; 75 cycles, single-end or paired-end). After trimming or removing low-quality data using Trimmomatic software (13), the qualified reads were analyzed using TopHat/Cufflinks (15) and Ensembl (https://www.ensembl.org/index.html) databases. Genes with low expression levels (<0.3 fragment per kilobase of transcript per million mapped reads) in normal- and LLC-BMDCs were excluded. The criteria of differentially expressed genes (DEGs) were -log10 (P-value) >1.3 and >2 fold change (FC), and -log10 (q-value) >0.6 and q<0.25 for stringent criteria.

miRNA database analysis

Putative genes targeted by candidate miRNAs were predicted using the miRmap database (http://cegg.unige.ch/mirmap) (16). The search criteria of the putative targeted genes by miRNAs were ‘mice species’ and a miRmap score of >99.0. Potential miRNA interactions were also searched using the miRmap TargetScan (http://www.targetscan.org/vert_71/) and miRDB (http://www.mirdb.org/) databases.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database analysis

The functions of DEGs were studied by analyzing the genes of interest according to previous studies (11,12). The Database for Annotation, Visualization and Integrated Discovery (https://david.ncifcrf.gov/) was used to conduct GO and KEGG database analyses. Briefly, a panel of potential genes was classified into clusters of associated biological functions, signaling pathways and diseases by calculating the similarity of global annotation profiles using the agglomeration algorithm method.

Results

Gene expression profiling and miRNA alterations in LLC-BMDCs versus normal BMDCs

The differentially downregulated (green spots, left panel) and upregulated (red spots, right panel) genes in the LLC-BMDCs versus normal-BMDCs were presented as a volcano plot (Fig. 2). A total of 820 genes with -log10 (P-value) >1.3 and >2 FC were chosen for further analyses (Fig. 2A). However, when using the following more stringent criteria: -log10 [q-value] >0.6 and FDR<0.25, only 139 genes were detected (Fig. 2B). Using the stringent criteria, too few genes were obtained to be analyzed.

Figure 2. Gene expression profiling in LLC-BMDCs versus normal-BMDCs. Differentially downregulated (green spots, left panel) and upregulated (red spots, right panel) genes expression in LLC-BMDCs versus normal-BMDCs were identified using the volcano plot. (A) Genes with -log10 (P-value) >1.3 and >2 FC were selected for further analyses. (B) Genes with -log10 (q-value) >0.6 and >2 FC. BMDCs, -bone marrow-derived cells; Ctrl, control; FC, fold change; FDR, false discovery rate; LLC, Lewis lung carcinoma.

The 820 genes expressed on the RNA sequencing heatmap (Fig. 3; left panel) revealed DEGs with FC>2 (either increase or decrease). The ‘gene’ Venn diagram demonstrated that 309 genes were upregulated and 511 were downregulated in the LLC-BMDCs compared with normal-BMDCs (Fig. 3). NGS miRNA heatmap (Fig. 3; right panel) analysis highlighted 161 differentially expressed miRNAs with FC>2 and P<0.05. Results from the ‘miRNA’ Venn diagram revealed the miRNAs with a threshold of RPM>1, which included 55 precursor miRNAs that were downregulated and 106 precursor miRNAs that were upregulated in BDMCs. The ‘targets’ Venn diagram exhibited the predicted genes of the miRNAs from the ‘miRNA’ Venn diagram using the miRmap online database. The selection threshold was a miRmap score ≥99.0. The intersection in the Venn diagram between ‘genes’ and ‘targets’ revealed 37 potential miRNA-mRNA interactions (Fig. 3), of which 34 included downregulated and three included upregulated genes (Table I).

Figure 3. There are 37 potential miRNA-mRNA interactions between ‘genes’ and ‘targets’. RNA sequencing heatmap (left panel) analysis presented 820 differentially expressed genes (both up/downregulated in LLC-BMDCs versus normal-BDMCs). The ‘genes’ Venn diagram revealed that 309 genes were upregulated and 511 genes were downregulated in the BMDCs. Next generation sequencing heatmap (right panel) analysis detected 161 differentially expressed miRNAs with FC>2. The ‘miRNA’ Venn diagram exposed the selected miRNAs by threshold of reads per million >1, with 55 microRNAs being downregulated and 106 miRNAs being upregulated in the BMDCs. The ‘targets’ Venn diagram revealed the predicted gene targets of miRNAs from the ‘miRNA’ Venn diagram using the miRmap online database. The selection threshold was a miRmap score ≥99.0. BMDCs, bone marrow-derived cells; Ctrl, control; FC, fold change; LLC, Lewis lung carcinoma; miRNA, microRNA.

Table I. Potential miRNA-mRNA interactions in LLC-BMDCs.

Table I.

Potential miRNA-mRNA interactions in LLC-BMDCs.

Gene symbol	Associated miRNA	Gene	log2 ratio (LLC/Ctrl)	Expression
BMPR1A	mmu-miR-135a-5p	Bone morphogenetic protein receptor, type 1A	2.80	Up
LRRC75B	mmu-miR-150-5p	Leucine rich repeat containing 75B	5.41	Up
FAM217B	mmu-miR-195a-5p	Family with sequence similarity 217, member B	2.93	Up
AXL	mmu-miR-1249-5p	AXL receptor tyrosine kinase	−2.89	Down
PHYHIP	mmu-miR-1249-5p	Phytanoyl-CoA hydroxylase interacting protein	−1.73	Down
FOSB	mmu-miR-1249-5p	FosB proto-oncogene, AP-1 transcription factor subunit	−1.75	Down
CACNA1E	mmu-miR-1249-5p	Calcium channel, voltage-dependent, R type, α1E subunit	−1.92	Down
MPP2	mmu-miR-1249-5p	Membrane protein, palmitoylated 2 (MAGUK p55 subfamily member 2)	−1.79	Down
ZBTB4	mmu-miR-1249-5p	Zinc finger and BTB domain containing 4	−2.02	Down
ADAM11	mmu-miR-1249-5p	A disintegrin and metallopeptidase domain 11	−1.69	Down
OTUB2	mmu-miR-3085-3p	OTU domain, ubiquitin aldehyde binding 2	−1.87	Down
NR1D2	mmu-miR-148a-5p	Nuclear receptor subfamily 1, group D, member 2	−1.93	Down
LYNX1	mmu-miR-1249-5p	Ly6/neurotoxin 1	−2.04	Down
Fhl1	mmu-miR-1249-5p	Four and a half LIM domains 1	−2.91	Down
Cd4	mmu-miR-1249-5p	CD4 antigen	−2.09	Down
EHD3	mmu-miR-1249-5p	EH-domain containing 3	−1.68	Down
DNASE1L3	mmu-miR-1249-5p	Deoxyribonuclease 1-like 3	−3.31	Down
SDC3	mmu-miR-1249-5p	Syndecan 3	−2.92	Down
NFATC2	mmu-miR-1249-5p	Nuclear factor of activated T cells 2	−1.33	Down
TNIK	mmu-miR-1198-3p	TRAF2 and NCK interacting kinase	−1.71	Down
PDE4B	mmu-miR-1249-5p	Phosphodiesterase 4B	−1.55	Down
BCL7A	mmu-miR-1249-5p	BAF chromatin remodeling complex subunit BCL7A	−1.50	Dpwn
DUSP16	mmu-miR-1249-5p	Dual specificity phosphatase 16	−1.64	Down
SORBS2	mmu-miR-21a-3p	Sorbin and SH3 domain containing 2	−4.21	Down
THY1	mmu-miR-1249-5p	Thy-1 cell surface antigen	−1.48	Down
CAPN5	mmu-miR-1249-5p	Calpain 5	−1.53	Down
TRIM58	mmu-miR-1249-5p	Tripartite motif-containing 58	−1.95	Down
CAMK1D	mmu-miR-1249-5p	Calcium/calmodulin-dependent protein kinase ID	−1.15	Down
FHDC1	mmu-miR-193b-3p	FH2 domain containing 1	−2.42	Down
SBK1	mmu-miR-1249-5p	SH3-binding kinase 1	−2.25	Down
DUSP18	mmu-miR-1249-5p	Dual specificity phosphatase 18	−1.65	Down
2900026A02RIK	mmu-miR-1249-5p	RIKEN cDNA 2900026A02 gene	−1.90	Down
MAPK11	mmu-miR-1249-5p	Mitogen-activated protein kinase 11	−2.78	Down
SH3PXD2A	mmu-miR-1249-5p	SH3 and PX domains 2A	−2.08	Down
CECR2	mmu-miR-1249-5p	Cat eye syndrome chromosome region, candidate 2	−2.43	Down
SOX4	mmu-miR-129-5p	SRY-box 4	−1.36	Down
ACSF2	mmu-miR-1249-5p	Acyl-CoA synthetase family member 2	−1.39	Down

[i] miRNA, microRNA.

KEGG analysis of DEGs in LLC-BMDCs

The three upregulated and 34 downregulated genes in the ‘genes’ and ‘targets’ diagram (Fig. 3) were mapped to KEGG pathways. The results identified eight potential pathways that may involve the LLC-BMDCs, including the ‘T-cell receptor signaling pathway’, ‘osteoclast differentiation’, ‘MAPK signaling pathway’ ‘VEGF signaling pathway’, ‘leukocyte transendothelial migration’, ‘signaling pathways regulating pluripotency of stem cells’, ‘oxytocin signaling pathway’ and ‘cell adhesion molecules (CAMs)’ (Table II).

Table II. KEGG analysis of differentially expressed genes predicted by miRmap in LLC-BMDCs.

Table II.

KEGG analysis of differentially expressed genes predicted by miRmap in LLC-BMDCs.

KEGG pathway	Count	P-value	Upregulated genes	Downregulated genes	Fold enrichment
T-cell receptor signaling pathway	3	0.01		Cd4, Mapk11, Nfatc2	18.56
Osteoclast differentiation	3	0.01		Fosb, Mapk11, Nfatc2	15.32
MAPK signaling pathway	3	0.05		Dusp16, Mapk11, Cacna1e	7.63
VEGF signaling pathway	2	0.08		Mapk11, Nfatc2	21.44
Leukocyte transendothelial migration	2	0.16		Mapk11, Thy1	10.63
Signaling pathways regulating pluripotency of stem cells	2	0.18	BMPR1A	Mapk11	9.32
Oxytocin signaling pathway	2	0.20		Camk1d, Nfatc2	8.14
Cell adhesion molecules (CAMs)	2	0.21		Cd4, Sdc3	7.94

Analysis of DEGs in LLC-BMDCs versus normal-BMDCs

DEGs were analyzed using GO in order to investigate the enriched functions of the 820 potential miRNA-mRNA interactions. The top 20 cellular components, top 20 biological processes and top 20 molecular functions of the DEGs in the LLC-BMDCs are presented in Fig. 4A-C. All top 20 results, including molecular functions, biological processes and cellular components, were statistically significant (P<0.05).

Figure 4. Analysis of differentially expressed genes in LLC-BMDCs. For enrichment analysis, the 820 potential mRNAs were analyzed using Gene Ontology. (A) Top 20 cellular components of the differentially expressed genes in the LLC-BMDCs. (B) Top 20 biological processes of the differentially expressed genes in the LLC-BMDCs. (C) Top 20 molecular functions of the differentially expressed genes in LLC-BMDCs. LLC, Lewis lung carcinoma, LLC-BMDCs, LLC-bone marrow-derived cells.

KEGG pathway enrichment analysis of 820 differentially expressed mRNAs in LLC-BMDCs versus normal-BMDCs

The top 20 pathways with the most significant enrichment are presented in Fig. 5, where each marker represents a KEGG pathway. The fold enrichment represents the ratio of the proportion of DEGs annotated to the pathway to all genes. The larger the fold enrichment score, the more significant the enrichment level of the DEGs in the pathway. In addition, the greater -log10(P-value), the more reliable the enrichment significance of the DEGs in the pathway. Two pathways (‘Hematopoietic cell lineage’ and ‘Primary immunodeficiency’) were most significantly enriched in LLC-BMDCs (Fig. 5).

Figure 5. KEGG pathway enrichment analysis of differentially expressed mRNAs. Top 20 pathways with the most significant enrichment are presented. Each marker represents a KEGG pathway (name captioned on the right). The fold enrichment represents the ratio of the proportion of differentially expressed genes annotated to the pathway to all genes. The larger the fold enrichment score, the more significant the enrichment level of the differentially expressed genes in the pathway. The y-coordinate [-log10 (P-value)] represents the significance of the study. KEGG, Kyoto Encyclopedia of Genes and Genomes.

Potential miRNA-mRNA interactions from miRmap, TargetScan and miRDB databases

The miRmap (selection criteria: miRmap score >97.0), TargetScan (selection criteria: TargetScan score >99.0) and miRDB (selection criteria: miRDB score >95.0) databases were used to explore the potential miRNA-mRNA interactions of the miRNAs with >2 FC (Fig. 6). By combining the results from the three databases, 23 miRNA-mRNA interactions were identified (Table III). Due to the insufficiency of gene interactions, the miRmap and miRDB scores were modified to >99 and >95.0, respectively.

Figure 6. Potential miRNA-mRNA interactions from miRmap, TargetScan and miRDB databases. Using miRmap (selection criteria: miRmap score >97.0), TargetScan (selection criteria: TargetScan score >99.0) and miRDB (selection criteria: miRDB score >95.0), potential miRNA-mRNA interactions of the differentially expressed miRNAs with 2-fold changes (downregulated and upregulated) in expression were explored. The interactions among these three databases reported 23 potential miRNA-mRNA interactions. miRNA, microRNA.

Table III. Potential miRNA-mRNA interactions of miRNAs with 2-fold change validated in miRmap, TargetScan and miRDB databases.

Table III.

Potential miRNA-mRNA interactions of miRNAs with 2-fold change validated in miRmap, TargetScan and miRDB databases.

miRNAs	Gene symbol
mmu-miR-33-5p	Abca1
mmu-miR-205-5p	Lrrk2
mmu-miR-129-5p	Cdc42ep3
mmu-miR-151-3p	Ago2
mmu-miR-205-5p	Cdh11
mmu-miR-342-3p	Fam208a
mmu-miR-150-5p	Prkar1a
mmu-miR-342-3p	Fam53c
mmu-miR-136-5p	Mtmr4
mmu-miR-411-5p	Pxdc1
mmu-miR-342-3p	Rictor
mmu-miR-150-5p	Prickle2
mmu-miR-129-5p	Glcci1
mmu-miR-129-5p	Chmp2b
mmu-miR-455-5p	Yipf6
mmu-miR-150-5p	Myb
mmu-miR-129-5p	ZFP36L1
mmu-miR-129-5p	BACH2
mmu-miR-455-5p	ADD3
mmu-miR-147-3p	NDUFA4
mmu-miR-33-5p	CNTN4
mmu-miR-409-3p	ZDHHC20
mmu-miR-33-5p	SLC12A5

[i] miRNA, microRNA.

Discussion

Patients with lung cancer and metastasis have a poor prognosis. Notably, PMN establishment is crucial for the development of distant metastases. The predicted mRNAs of specific miRNAs imply potential regulatory mechanisms of miRNA targets in the formation of a PMN. In the present study, ‘gene’ and miRNA-predicted ‘target’ genes were matched to strengthen the significance of the selected genes that mediate PMN. Eight potential pathways that may be involved in the bone marrow modification observed during lung cancer progression were revealed. These eight pathways were derived from the common genes of the ‘genes’ and ‘targets’, as follows: ‘T-cell receptor signaling pathway’, ‘MAPK signaling pathway’, ‘osteoclast differentiation’, ‘VEGF signaling pathway’, ‘leukocyte transendothelial migration’, ‘signaling pathways regulating pluripotency of stem cells’, ‘oxytocin signaling pathway’ and ‘cellular adhesion molecules (CAMs)’. In addition, the present study analyzed the KEGG pathways derived from mRNAs with >2 FC in the LLC-BMDCs compared with the control-BDMCs. Results demonstrated that three pathways were identical, including the ‘T-cell receptor signaling pathway’, ‘MAPK signaling pathway’ and ‘cellular adhesion molecules (CAMs)’.

The T-cell receptor (TCR) signaling pathway represents a branching network, which is initiated by cognate peptide-major histocompatibility complex molecules, and is associated with the mitogen-activated protein kinase (MAPK) and nuclear factor-κB (NF-κB) signaling pathways (17). It has also been reported that the TCR signaling pathway involves the mobilization of transcription factors. Subsequently, it serves a crucial role in T-cell gene expression and is essential for T-cell growth and differentiation (17). Furthermore, it induces expression of the inducible costimulatory molecule (ICOS) and programmed cell death 1 (PD1), which are members of the CD28 family. ICOS and PD1 can control the sustained phase of T-cell signaling (18), whereas PD1 inhibitors can target the binding between T cell PD1 receptors and PD-ligands 1 and 2. PD1 inhibitors have also been reported to block inhibitory signaling, which results in the activation of T-cell effector function, which provides an anticancer potential for activated T cells (19).

The Ras/Raf/mitogen-activated protein kinase kinase/extracellular signal-regulated kinase (ERK) signaling pathway, which is one of the major intracellular axes that regulate intracellular signaling trafficking, is associated with cell proliferation, growth, invasion, metastasis, resistance to apoptosis and angiogenesis in lung cancer (20). It has been reported that the MAPK pathway interferes in non-small cell lung cancer tumorigenesis through terminal differentiation-induced noncoding RNA (21), although the role of the MAPK/ERK pathway in lung cancer treatment remains unclear. The therapeutic inhibition of elements from this pathway has been reported to be beneficial in cancer treatment, including colorectal and ovarian cancer (22). However, treatment with small inhibitors that specifically target proteins from the MAPK/ERK pathway can lead to the development of secondary malignancies. For example, Raf inhibitors have been reported to potentially induce abnormal skin cell proliferation and cause secondary squamous cell cancers (22).

CAMs, which comprise cadherins, integrins, selectins and members of the immunoglobulin family, are important components involved in cell-to-cell and cell-to-extracellular matrix anchoring. Their main roles are to maintain cell and tissue structure, cell signaling, tissue repair and wound healing (23). In addition, integrins serve an important role in platelet aggregation, hematopoietic cell mobilization, neoangiogenesis and stromal function. They are also associated with bone metastasis, which suggests that they may represent potential novel therapeutic targets in the prevention and treatment of bone metastasis (24). For example, CAMs, including CD44, N-cadherin, neural cell adhesion molecule and integrins, are involved in the metastatic cascade observed in a metastatic neuroblastoma model (25).

Osteoclast differentiation is associated with bone metastasis in lung cancer. Non-small cell lung cancer can metastasize to the bone, which results in bone osteolytic lesions via osteoclast lineage activation. Through the downregulation of the receptor activator of NF-κB signaling pathway, osteoclast differentiation attenuates osteoclastogenesis and reduces osteolysis following bone metastasis (26).

Additional pathways predicted in this study, such as VEGF, leukocyte transendothelial migration, cancer stem cells, etc., are associated with lung cancer metastasis. To support cellular function, tumor cells release VEGF that stimulates angiogenesis through vascular sprouting, intussusception and incorporation of bone marrow-derived endothelial precursors (27). In gastric cancer, celecoxib has been used as an inhibitor of cysteine and glycine rich protein 1, thrombospondin 1, myosin light chain 9, filamin A, actinin α1, vinculin, laminin subunit γ2 and claudin 1 expression, and can also suppress leukocyte transendothelial migration and focal adhesion, which highlights its anti-gastric cancer effect (28). Furthermore, cancer stem cells possess self-renewal ability and are involved in the initiation of cancer, malignant transformation and metastatic progression, and in the post-treatment recurrence of various types of human cancer (29). Understanding the underlying mechanisms of transcription factors associated with pluripotency is therefore crucial to clarify human carcinogenesis. In addition, the hypothalamic nonapeptide, oxytocin, serves as a tumor growth regulator by activating specific G-coupled transmembrane receptor through oxytocin receptor. Oxytocin attenuates proliferation of cancer cells that originate from epithelium, nerves and bone (30).

In conclusion, results from the present study demonstrated that bone marrow may serve a crucial role in mediating lung cancer metastasis. Data revealed that numerous pathways were involved in the LLC-BMDC microenvironment, including the ‘T-cell receptor signaling pathway’, ‘MAPK signaling pathway’, ‘osteoclast differentiation’, ‘VEGF signaling pathway’, ‘leukocyte transendothelial migration’, ‘signaling pathways regulating pluripotency of stem cells’, ‘oxytocin signaling pathway’ and ‘cellular adhesion molecules (CAMs)’. These eight pathways were all associated with cancer progression. In addition, the present study provided macroscopic evidence of the association between lung cancer cells and BMDCs (Fig. 7). Present findings suggested that these genetic interactions may be associated with lung cancer cell progression. Molecular alterations in gene expression may therefore represent a novel signature in lung cancer, which may be used to develop diagnostic and therapeutic strategies for patients with lung cancer and bone metastasis. Further investigation is required to study the role of BMDCs in the bone microenvironment in lung cancer metastasis, and to identify the main cells that mediate bone metastasis in lung cancer.

Figure 7. Potential pathways involved in the microenvironment of lung cancer and bone marrow-derived cells. The injection of Lewis lung carcinoma cells into the mice modified bone marrow-derived cells and initiated a pre-metastatic niche. Through next generation sequencing, eight potential pathways, including the ‘T-cell receptor signaling pathway’, ‘osteoclast differentiation’, ‘MAPK signaling pathway’ (red font, P<0.05), ‘VEGF signaling pathway’, ‘leukocyte transendothelial migration’, ‘signaling pathways regulating pluripotency of stem cells’, ‘oxytocin signaling pathway’, and ‘cell adhesion molecules (CAMs)’ (blue font, P>0.05) were identified. These pathways may mediate the pre-metastatic niche in lung cancer metastasis.

Acknowledgements

The authors would like to thank Mr. Chi-Hsien Chou and other staff members of the Center for Research Resources and Development of Kaohsiung Medical University for their assistance in bioinformatics analysis.

Funding

This study was supported by the Ministry of Science and Technology (grants nos. MOST 107-2314-B-037-107-MY3, MOST 106-2314-B-037-016-MY2 and MOST 106-2314-B-037-064), the Kaohsiung Medical University Hospital (grants nos. KMUHS10701, KMUHS10712 and KMUH103-3R09) and the Kaohsiung Medical University (grants nos. 105KMUOR05, KMU-DK108003 and KMU-DK108008).

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Authors' contributions

WAC and JYH designed the study. WAC, YMT, YCT, CYW, KFC, CTL, JYH, YLH, and PLK analyzed the data and interpreted the results. WAC, YMT and JYH wrote the manuscript. The final version of the manuscript has been read and approved by all authors.

Ethical approval and consent to participate

Animal experiments were performed according to the Institutional Animal Care and Use Committee, and with the approval of the Animal Care and Use Committee of the School of Kaohsiung Medical University.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References