OL Oncology Letters 1792-1074 1792-1082 D.A. Spandidos 10.3892/ol.2019.10085 OL-0-0-10085 Articles Differential expression profiles of the transcriptome in bone marrow-derived cells in lung cancer revealed by next generation sequencing and bioinformatics ChangWei-An 1 2 TsaiYing-Ming 1 2 3 TsaiYu-Chen 2 4 WuCheng-Ying 1 ChangKuo-Feng 5 LienChi-Tun 6 HungJen-Yu 2 3 HsuYa-Ling 4 KuoPo-Lin 1 7 Graduate Institute of Clinical Medicine, College of Medicine, Kaohsiung Medical University, Kaohsiung 807, Taiwan, R.O.C. Division of Pulmonary and Critical Care Medicine, Department of Internal Medicine, Kaohsiung Medical University Hospital, Kaohsiung 807, Taiwan, R.O.C. School of Medicine, College of Medicine, Kaohsiung Medical University, Kaohsiung 807, Taiwan, R.O.C. Graduate Institute of Medicine, College of Medicine, Kaohsiung Medical University, Kaohsiung 807, Taiwan, R.O.C. Welgene Biotech, Inc., Taipei 115, Taiwan, R.O.C. Department of Internal Medicine, Kaohsiung Municipal United Hospital, Kaohsiung 804, Taiwan, R.O.C. Center for Infectious Disease and Cancer Research Kaohsiung Medical University, Kaohsiung 807, Taiwan, R.O.C. Correspondence to: Dr Jen-Yu Hung, Division of Pulmonary and Critical Care Medicine, Department of Internal Medicine, Kaohsiung Medical University Hospital, 16ES, 100 Tz-You 1st Road, Kaohsiung 807, Taiwan, R.O.C., E-mail: jenyuhung@gmail.com 05 2019 28 02 2019 17 5 4341 4350 08062018 04022019 Copyright: © Chang et al. 2019 This is an open access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.

A pre-metastatic niche (PMN) facilitates cancer metastasis through mobilization and recruitment of bone marrow-derived cells (BMDCs) and associated factors. In bone marrow, hematogenous cells, including osteoclasts, macrophages and lymphocytes, and mesenchymal cells, including mesenchymal stem cells, osteoblasts and adipocytes, are involved in PMN formation. Patients with lung cancer and metastasis have a poor prognosis and shortened median survival time. Bone marrow has been considered fertile ground for dormant and proliferating tumor cells, and mobilizing and recruiting BMDCs and immune cells can establish a PMN. However, the role of BMDCs in PMN formation is not yet fully understood. The present study aimed to investigate the association between BMDCs and PMN in bone marrow tissue samples. The results demonstrated that bone marrow served an important role in lung cancer progression and that eight pathways were potentially involved, including ‘T-cell receptor signaling pathway’, ‘osteoclast differentiation’, ‘MAPK signaling pathway’, ‘VEGF signaling pathway’, ‘leukocyte transendothelial migration’, ‘signaling pathways regulating the pluripotency of stem cells’, ‘oxytocin signaling pathway’ and ‘cell adhesion molecules (CAMs)’. In addition, the present study investigated the role of BMDCs in facilitating lung cancer metastasis. In conclusion, the results from the present study suggested that molecular alterations in gene expression may provide a novel signature in lung cancer, which may aid in the development of novel diagnostic and therapeutic strategies for patients with lung cancer and bone metastasis.

 bone marrow-derived cells next generation sequencing bioinformatics lung cancer 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).

 Materials and methods <sec> <title>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 <sec> <title>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.

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).

 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).

 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).

 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).

 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.

 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.

 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 KaplanRNRibaRDZacharoulisSBramleyAHVincentLCostaCMacDonaldDDJinDKShidoKKernsSAVEGFR1-positive haematopoietic bone marrow progenitors initiate the pre-metastatic nicheNature438820827200510.1038/nature0418616341007 RiihimäkiMHemminkiAFallahMThomsenHSundquistKSundquistJHemminkiKMetastatic sites and survival in lung cancerLung Cancer867884201410.1016/j.lungcan.2014.07.02025130083 ColemanREClinical features of metastatic bone disease and risk of skeletal morbidityClin Cancer Res126243s6249s200610.1158/1078-0432.CCR-06-093117062708 PeinadoHZhangHMateiIRCosta-SilvaBHoshinoARodriguesGPsailaBKaplanRNBrombergJFKangYPre-metastatic niches: Organ-specific homes for metastasesNat Rev Cancer17302317201710.1038/nrc.2017.628303905 KaplanRNPsailaBLydenDBone marrow cells in the ‘pre-metastatic niche’: Within bone and beyondCancer Metastasis Rev25521529200610.1007/s10555-006-9036-917186383 RoatoIBone metastases: When and how lung cancer interacts with boneWorld J Clin Oncol5149155201410.5306/wjco.v5.i2.14924829862 BergfeldSADeClerckYABone marrow-derived mesenchymal stem cells and the tumor microenvironmentCancer Metastasis Rev29249261201010.1007/s10555-010-9222-720411303 RutkowskiMRSvoronosNPerales-PuchaltAConejo-GarciaJRThe tumor macroenvironment: Cancer-promoting networks beyond tumor bedsAdv Cancer Res128235262201510.1016/bs.acr.2015.04.01126216635 LiuYCaoXCharacteristics and significance of the pre-metastatic nicheCancer Cell30668681201610.1016/j.ccell.2016.09.01127846389 LiuXQuanNImmune cell isolation from mouse femur bone marrowBio Protoc5e1631201510.21769/BioProtoc.163127441207 SheuCCTsaiMJChenFWChangKFChangWAChongIWKuoPLHsuYLIdentification of novel genetic regulations associated with airway epithelial homeostasis using next-generation sequencing data and bioinformatics approachesOncotarget88267482688201710.18632/oncotarget.1975229137293 ChenSCChenFWHsuYLKuoPLSystematic analysis of transcriptomic profile of renal cell carcinoma under long-term hypoxia using next-generation sequencing and bioinformaticsInt J Mol Sci18E2657201710.3390/ijms1812265729215599 BolgerAMLohseMUsadelBTrimmomatic: A flexible trimmer for illumina sequence dataBioinformatics3021142120201410.1093/bioinformatics/btu17024695404 FriedlanderMRMackowiakSDLiNChenWRajewskyNmiRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal cladesNucleic Acids Res403752201210.1093/nar/gkr68821911355 TrapnellCRobertsAGoffLPerteaGKimDKelleyDRPimentelHSalzbergSLRinnJLPachterLDifferential gene and transcript expression analysis of RNA-seq experiments with TopHat and CufflinksNat Protoc7562578201210.1038/nprot.2012.01622383036 VejnarCEZdobnovEMMiRmap: Comprehensive prediction of microRNA target repression strengthNucleic Acids Res401167311683201210.1093/nar/gks90123034802 BrownlieRJZamoyskaRT cell receptor signalling networks: Branched, diversified and boundedNat Rev Immunol13257269201310.1038/nri340323524462 HuseMThe T-cell-receptor signaling networkJ Cell Sci12212691273200910.1242/jcs.04276219386893 NiyongereSSaltosAGrayJEImmunotherapy combination strategies (non-chemotherapy) in non-small cell lung cancerJ Thorac Dis10S433S450201810.21037/jtd.2017.12.12029593889 ReungwetwattanaTDyGKTargeted therapies in development for non-small cell lung cancerJ Carcinog1222201310.4103/1477-3163.12397224574860 ZhuZJHeJKTINCR facilitates non-small cell lung cancer progression through BRAF-activated MAPK pathwayBiochem Biophys Res Commun497971977201810.1016/j.bbrc.2018.02.05929427662 BurottoMChiouVLLeeJMKohnECThe MAPK pathway across different malignancies: A new perspectiveCancer12034463456201410.1002/cncr.2886424948110 FarahaniEPatraHKJangamreddyJRRashediIKawalecMRao ParitiRKBatakisPWiechecECell adhesion molecules and their relation to (cancer) cell stemnessCarcinogenesis35747759201410.1093/carcin/bgu04524531939 SchneiderJGAmendSRWeilbaecherKNIntegrins and bone metastasis: Integrating tumor cell and stromal cell interactionsBone485465201110.1016/j.bone.2010.09.01620850578 SchwankhausNGathmannCWickleinDRieckenKSchumacherUValentinerUCell adhesion molecules in metastatic neuroblastoma modelsClin Exp Metastasis31483496201410.1007/s10585-014-9643-824549749 IhnHJKimJABaeYCShinHIBaekMCParkEKAfatinib ameliorates osteoclast differentiation and function through downregulation of RANK signaling pathwaysBMB Rep50150155201710.5483/BMBRep.2017.50.3.22328256196 NiuGChenXVascular endothelial growth factor as an anti-angiogenic target for cancer therapyCurr Drug Targets1110001017201010.2174/13894501079159139520426765 JinGHXuWShiYWangLBCelecoxib exhibits an anti-gastric cancer effect by targeting focal adhesion and leukocyte transendothelial migration-associated genesOncol Lett1223452350201610.3892/ol.2016.497627698798 KashyapVRezendeNCScotlandKBShafferSMPerssonJLGudasLJMonganNPRegulation of stem cell pluripotency and differentiation involves a mutual regulatory circuit of the NANOG, OCT4, and SOX2 pluripotency transcription factors with polycomb repressive complexes and stem cell microRNAsStem Cells Dev1810931108200910.1089/scd.2009.011319480567 CassoniPSapinoAMarroccoTChiniBBussolatiGOxytocin and oxytocin receptors in cancer cells and proliferationJ Neuroendocrinol16362364200410.1111/j.0953-8194.2004.01165.x15089975

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.

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.

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.

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 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. 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.

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.

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

miRNA, microRNA.

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

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

miRNA, microRNA.