Abstract

   Schwann cells were identified in the tumor surrounding area prior to
   initiate the invasion process underlying connective tissue. These cells
   promote cancer invasion through direct contact, while paracrine
   signaling and matrix remodeling are not sufficient to proceed.
   Considering the intertwined structure of signaling, regulatory, and
   metabolic processes within a cell, we employed a genome-scale
   biomolecular network. Accordingly, a meta-analysis of Schwann cells
   associated transcriptomic datasets was performed, and the core
   information on differentially expressed genes (DEGs) was obtained by
   statistical analyses. Gene set over-representation analyses was
   performed on core DEGs to identify significantly functional and pathway
   enrichment analysis between Schwann cells and, lung adenocarcinoma
   (LUAD) and lung squamous cell carcinoma (LUSC). DEGs were further
   integrated with genome-scale human biomolecular networks. miRNAs were
   proposed by the reconstruction of a transcriptional and
   post-transcriptional regulatory network. Moreover, microarray-based
   transcriptome profiling was performed, and the prognostic power of
   selected dedifferentiated Schwann cell biomolecules was predicted. We
   observed that pathways associated with Schwann cells dedifferentiation
   was overexpressed in lung cancer samples. However, genes associated
   with Schwann cells migration inhibition system were downregulated.
   Besides, miRNA targeting those pathways were also deregulated. In this
   study, we report valuable data for further experimental and clinical
   analysis, because the proposed biomolecules have significant potential
   as systems biomarkers for screening or for therapeutic purposes in
   perineural invasion of lung cancer.

   Keywords: bioinformatic, lung squamous cell carcinoma, lung
   adenocarcinoma, neuroactive ligand-receptor interaction

INTRODUCTION

   Peripheral nervous system plays an important role in neoplastic
   invasion. The glial cells present in this system, specifically Schwann
   cells, have been associated with the process known as perineural
   invasion. Perineural invasion occurs mainly in pancreas, prostate,
   breast, lung and head and neck cancers; and causes morbidity of patient
   through pain and paralysis, associated with a high risk of local
   recurrence and decreasing survival patients rates [[49]1–[50]3].
   Understanding how cancer invades nerves is an essential step in
   developing more effective treatment strategies. The main questions are
   how the neoplastic cells interact with the cells of the nervous system
   and how they acquire mobile and invasive characteristics from these
   interactions.

   Schwann cells performs several roles, depends on their ability to
   reversibly dedifferentiate and re-differentiate into subtypes with
   different phenotypes. After the nerve damage, Schwann cells become
   dedifferentiated, lose their capacity to myelinate and, promote
   neuronal orientation during repair [[51]3]. This is accompanied by the
   reexpression of proteins lost during the myelinizing differentiation
   program, such as the glial fibrillary acidic protein (GFAP) and the
   neural cell adhesion molecule 1 (NCAM1) [[52]1, [53]4]. Paracrine
   signaling has been implicated in perineural invasion, with factors
   secreted by nerves, including glial cell line derived neurotrophic
   factor (GDNF), increasing the invasion of tumor cells along the nerves
   [[54]5, [55]6].

   Schwann cells were identified in the tumor surrounding area prior to
   initiate the invasion process underlying connective tissue. These cells
   promote cancer invasion through direct contact, while paracrine
   signaling and matrix remodeling are not sufficient to proceed [[56]4].
   These cells stimulate tumor cells at the cell-cell contact sites and
   promote detachment and spread of individual cancer cells from
   neighboring tissue. Physical contact between Schwann cells and tumor
   cells is necessary to potentiate the process of tumor invasion, by
   induction of specific chemokines, such as CX3CL1 and its receptor
   CX3CL1 [[57]2, [58]7, [59]8]. Schwann cells are also able to degrade
   the extracellular matrix, forming tunnels or bands coated with laminin.
   Tumor cells use these structures to migrate during tissue invasion.
   These activities strongly promote tumor spread and are dependent on
   NCAM1 expression by Schwann cells [[60]2, [61]3, [62]7, [63]8].

   Considering the intertwined structure of signaling, regulatory, and
   metabolic processes within a cell, we employed a genome-scale
   biomolecular network. Accordingly, a meta-analysis of Schwann cells
   associated transcriptomic datasets was performed, and the core
   information on differentially expressed genes (DEGs) was obtained by
   statistical analyses. Gene set over-representation analyses was
   performed on core DEGs to identify significantly functional and pathway
   enrichment analysis. DEGs were further integrated with genome-scale
   human biomolecular networks. miRNAs were proposed by the reconstruction
   of a transcriptional and post-transcriptional regulatory network.
   Moreover, microarray-based transcriptome profiling was performed, and
   the prognostic power of selected dedifferentiated Schwann cell
   biomolecules was predicted. Consequently, this systematic study reports
   candidate biomolecules that can be considered as biomarkers or
   potential targets for further experimental and clinical trials for
   perineural invasion risk in lung cancer patients.

RESULTS

Identification of differentially expressed genes, gene ontology enrichment,
and functional classification

   We obtained gene expression data belonging to specimens across two
   cancer types (LUAD and LUSC) from TCGA; these data were preprocessed
   using standard methods. After the identification of DEGs, a DAVID
   analysis was performed using them. Downregulated genes in both lung
   cancers belong to “complement and coagulation cascades” and “cell
   adhesion molecules (CAMs)” pathways. While “cell cycle” and “DNA
   replication” were the most active pathways in both lung cancers
   ([64]Tables 1–[65]4).

Table 1. Functional annotation analysis of downregulated differentially
expressed genes (DEGs) in lung adenocarcinoma (LUAD) datasets using the DAVID
tool.

   Term Count P Value Genes Pop Hits Fold Enrichment
   hsa04610:Complement and coagulation cascades 17 2,02E-05 F11, C7, A2M,
   F10, C5AR1, MASP1, C6, F8, SERPING1, C1QA, C8B, VWF, C1QB, THBD, CFD,
   PROS1, CPB2 69 3,423022096
   hsa04270:Vascular smooth muscle contraction 21 1,10E-04 RAMP3, ADCY4,
   RAMP2, NPR1, PRKCE, PRKG1, MYL9, AGTR1, ACTG2, PRKCQ, PTGIR, GNAQ,
   ADCY9, GUCY1A2, PLA2G1B, MYH11, CALCRL, PLA2G3, PLA2G5, MYLK, PPP1R14A
   112 2,605020492
   hsa04514:Cell adhesion molecules (CAMs) 21 0,00103288 SELP, ITGAL,
   CLDN18, PTPRM, CADM1, ICAM2, CLDN5, NFASC, CLDN11, HLA-E, CDH5,
   SIGLEC1, CD34, ITGA8, ESAM, JAM2, JAM3, NEGR1, SELE, SELPLG, SPN 132
   2,210320417
   hsa04080:Neuroactive ligand-receptor interaction 32 0,00245857 F2RL3,
   THRB, LEPR, PTH1R, FPR1, GRIN3B, FPR2, VIPR1, APLNR, AGTR1, EDNRB,
   AGTR2, PTGIR, S1PR1, NMUR1, S1PR4, CNR1, P2RY1, S1PR5, CALCRL, GHR,
   GRID1, C5AR1, PTGER4, PTGFR, ADRB2, ADRB1, P2RX1, SSTR1, GRIA1, P2RY14,
   CTSG 256 1,736680328
   hsa04614:Renin-angiotensin system 6 0,0055758 AGTR1, ACE, AGTR2, MME,
   CPA3, CTSG 17 4,903567985
   hsa02010:ABC transporters 9 0,0117617 ABCA8, ABCC9, ABCA9, ABCB1, CFTR,
   ABCA3, ABCG1, ABCA6, ABCG2 44 2,841840537
   hsa04670:Leukocyte transendothelial migration 16 0,02073563 ITGAL,
   CLDN18, NCF2, NCF1, CLDN5, CLDN11, CXCL12, CDH5, MYL9, CYBB, ESAM,
   PIK3R5, RAPGEF4, JAM2, JAM3, PIK3R1 118 1,883856627
   hsa03320:PPAR signaling pathway 11 0,02420575 LPL, CD36, CYP27A1,
   SORBS1, OLR1, PPARG, RXRG, FABP4, SCD5, ACADL, FABP5 69 2,21489665
   hsa04062:Chemokine signaling pathway 22 0,02578677 ADCY4, CCL2, FGR,
   NCF1, PREX1, HCK, CXCL3, CXCL2, CXCR1, GNG11, CXCR2, CXCL12, CCL14,
   CCL23, PPBP, ADCY9, ARRB1, RASGRP2, CX3CR1, PIK3R5, GRK5, PIK3R1 187
   1,634522662
   hsa04060:Cytokine-cytokine receptor interaction 28 0,03410983 CSF3,
   ACVRL1, CCL2, PDGFB, CXCL3, LEPR, CXCL2, BMPR2, CXCR1, CXCR2, TNFSF13,
   TNFSF12, IL7R, CXCL12, CCL23, PLEKHO2, FIGF, GHR, IL18R1, IL6, BMP2,
   FLT4, TGFBR2, LIFR, CCL14, PPBP, IL20RA, CX3CR1, IL3RA 262 1,484795395
   hsa05020:Prion diseases 7 0,0361229 C1QA, EGR1, C1QB, C8B, C7, IL6, C6
   35 2,778688525
   hsa00590:Arachidonic acid metabolism 9 0,04474724 ALOX15, PTGIS, PTGDS,
   GPX3, PLA2G1B, LTC4S, ALOX5, PLA2G3, HPGDS, PLA2G5 56 2,232874707
   hsa04350:TGF-beta signaling pathway 12 0,04546484 BMP2, SMAD9, ACVRL1,
   CDKN2B, ZFYVE9, SMAD6, LEFTY2, TGFBR2, BMPR2, ID4, DCN, ID3 87
   1,916336914
   [66]Open in a new tab

   The results revealed that downregulated consistent DEGs indicated the
   genes that were associated with the top four pathways: “Complement and
   coagulation cascades”, “Vascular smooth muscle contraction”, “Cell
   adhesion molecules (CAMs)” and, “Neuroactive ligand-receptor
   interaction”.

Table 4. Functional annotation analysis of upregulated differentially
expressed genes (DEGs) in Lung squamous cell carcinoma (LUSC) datasets using
the DAVID tool.

   Term Count P Value Genes Pop Hits Fold Enrichment
   hsa04110:Cell cycle 54 7,45E-13 E2F1, E2F2, E2F3, DBF4, PKMYT1, TTK,
   PTTG1, CCNE2, CCNE1, CDC45, MCM7, CDKN2A, CCNA1, CCNA2, CDC7, CDK1,
   CDC6, RBL1, SKP2, CDK6, ESPL1, MCM2, MCM3, CDK4, MCM4, MCM5, ORC1L,
   CDK2, MCM6, MAD2L1, BUB1B, ORC5L, ANAPC7, MAD2L2, YWHAZ, PRKDC, CHEK1,
   CHEK2, SFN, ORC6L, BUB1, TFDP2, TFDP1, CDC20, ATR, CDC25C, CDC25A,
   CCNB1, HDAC2, CCNB2, HDAC1, PLK1, PCNA, SMC1B 125 2,715352287
   hsa03030:DNA replication 25 9,11E-12 RNASEH1, POLA2, RPA3, PRIM1, MCM7,
   POLE2, POLE3, PRIM2, FEN1, LIG1, POLE, MCM2, RNASEH2A, MCM3, MCM4,
   MCM5, MCM6, RFC5, DNA2, RFC3, RFC4, RFC2, POLD1, POLD2, PCNA 36
   4,364956737
   hsa04115:p53 signaling pathway 27 6,14E-06 BID, LRDD, RPRM, CHEK1, SFN,
   CHEK2, PMAIP1, GTSE1, SESN3, CCNE2, CCNE1, CDKN2A, BAI1, TP53AIP1,
   CDK1, CYCS, CDK6, ATR, CDK4, CDK2, TP73, CCNB1, CCNB2, SERPINB5, RRM2,
   PERP, IGFBP3 68 2,495728205
   hsa00480:Glutathione metabolism 21 3,38E-05 GSTA1, ODC1, GCLC, GGCT,
   PGD, GCLM, GSS, GSTM1, GPX2, GSTM2, GSR, GSTM3, GSTM4, G6PD, OPLAH,
   RRM2, RRM1, IDH2, GSTO2, GPX7, SMS 50 2,639925834
   hsa03410:Base excision repair 16 1,29E-04 APEX2, UNG, NEIL3, LIG1,
   POLE, LIG3, POLB, SMUG1, POLE2, POLE3, POLD1, POLD2, PCNA, TDG, PARP1,
   FEN1 35 2,873388663
   hsa05217:Basal cell carcinoma 21 1,64E-04 FZD9, WNT5A, DVL3, WNT10A,
   WNT16, WNT10B, WNT5B, LEF1, FZD3, GLI2, FZD7, FZD6, GLI1, WNT2B, SMO,
   FZD10, WNT7B, WNT4, WNT3, WNT11, PTCH2 55 2,399932577
   hsa03430:Mismatch repair 12 3,26E-04 RFC5, EXO1, MSH6, RFC3, RFC4,
   RFC2, MSH2, POLD1, LIG1, POLD2, PCNA, RPA3 23 3,279410974
   hsa00010:Glycolysis / Gluconeogenesis 21 6,17E-04 ALDOA, LDHA, ALDOC,
   PGAM1, HK2, PFKP, ADH1C, ALDH3B2, PGAM2, ADH7, PCK1, ALDH3A1, PGM2,
   GPI, TPI1, PKM2, ENO2, ENO3, PGK1, GAPDH, ENO1 60 2,199938195
   hsa05200:Pathways in cancer 75 6,70E-04 FGF19, E2F1, E2F2, E2F3, HRAS,
   PGF, MMP9, FGF11, FGF12, GLI2, MMP1, ACVR1C, GLI1, CCNE2, CCNE1, WNT4,
   WNT3, CDKN2A, SLC2A1, TGFA, NOS2, CCNA1, FGF3, EGFR, WNT10A, RET, PLD1,
   WNT10B, CYCS, SKP2, LEF1, FADD, CDK6, CDK4, CDK2, RAD51, JUP, SMO,
   PIAS3, LAMC2, WNT11, WNT5A, BID, CKS1B, WNT16, WNT5B, FGFR3, EGLN3,
   TFG, CDH1, LAMB3, RAC3, EGF, TRAF4, PIK3R2, FZD9, MSH6, DVL3, MSH2,
   BRCA2, ITGA2, FZD3, BIRC5, COL4A6, FZD7, FZD6, WNT2B, LAMA1, CBLC,
   FZD10, WNT7B, HDAC2, ITGA6, HDAC1, PTCH2 328 1,437241852
   hsa04114:Oocyte meiosis 32 6,81E-04 YWHAZ, ADCY2, PKMYT1, AURKA, PTTG1,
   PRKX, CCNE2, CCNE1, CALML3, FBXO43, BUB1, STAG3, FBXO5, CAMK2B, CALML5,
   CDK1, SGOL1, CHP2, IGF2, ESPL1, CDC20, CDC25C, CDK2, CCNB1, RPS6KA6,
   MAD2L1, CCNB2, MAPK12, PLK1, ANAPC7, MAD2L2, SMC1B 110 1,828520058
   hsa00240:Pyrimidine metabolism 28 0,00129063 POLR2H, CTPS, PNPT1,
   DTYMK, CAD, POLA2, POLR2D, TK1, PRIM1, TYMS, POLE2, NT5M, POLE3, PRIM2,
   ENTPD3, UCK2, POLR3G, POLE, POLR1C, POLR2J3, NME4, UMPS, NME2,
   NME1-NME2, NME1, RRM2, POLD1, POLD2, RRM1, TXNRD1 95 1,852579533
   hsa00250:Alanine, aspartate and glutamate metabolism 13 0,00179493
   GLS2, GOT2, ADSSL1, ALDH4A1, ADSL, IL4I1, GPT, CAD, ASNS, CPS1, GAD1,
   GPT2, PPAT 31 2,635870649
   hsa04340:Hedgehog signaling pathway 19 0,00182644 WNT5A, WNT10A, WNT16,
   WNT10B, WNT5B, GLI2, PRKX, ZIC2, GLI1, WNT2B, SMO, WNT7B, WNT4, WNT3,
   WNT11, PTCH2, BMP7, BMP8B, BMP8A 56 2,132593149
   hsa04512:ECM-receptor interaction 25 0,0021767 TNC, COL3A1, ITGB4,
   ITGA11, ITGA2, COL2A1, COL5A2, COL5A1, COL4A6, HMMR, LAMA1, LAMB3,
   SDC1, ITGA6, ITGB8, COMP, COL1A2, LAMC2, SV2B, SV2A, COL1A1, COL11A2,
   THBS2, COL11A1, SPP1 84 1,870695744
   hsa03440:Homologous recombination 12 0,00241416 XRCC3, XRCC2, BLM,
   POLD1, POLD2, EME1, SHFM1, BRCA2, RAD54B, RAD54L, RAD51, RPA3 28
   2,693801872
   hsa00360:Phenylalanine metabolism 10 0,00433744 GOT2, DDC, NAA20,
   LCLAT1, ALDH3B2, IL4I1, PAH, PNPLA3, MIF, ALDH3A1 22 2,857062591
   hsa00230:Purine metabolism 37 0,00801809 XDH, POLR2H, GDA, ADCY2,
   PNPT1, POLA2, POLR2D, HPRT1, AK3L1, ADA, PPAT, PRIM1, POLE2, ATIC,
   POLE3, NT5M, PRIM2, ENTPD3, ADCY10, ENTPD2, POLR3G, ADSSL1, POLE,
   POLR1C, GMPS, POLR2J3, GART, NME4, NME2, NME1-NME2, NME1, RRM2, POLD1,
   PKM2, POLD2, RRM1, ADSL, PAICS, PRPS2 153 1,520031993
   hsa00980:Metabolism of xenobiotics by cytochrome P450 18 0,01003269
   GSTA1, CYP2C18, CYP2S1, ADH1C, ALDH3B2, ADH7, CYP2E1, UGT1A1, ALDH3A1,
   AKR1C3, GSTM1, UGT1A7, GSTM2, AKR1C2, UGT1A10, UGT1A6, UGT1A9, GSTM3,
   GSTM4, UGT1A8, UGT2A1, GSTO2, AKR1C1 60 1,88566131
   hsa00030:Pentose phosphate pathway 10 0,01140598 PGM2, ALDOA, GPI,
   TALDO1, G6PD, ALDOC, PGD, PFKP, TKTL1, PRPS2 25 2,51421508
   hsa00330:Arginine and proline metabolism 16 0,01525979 PYCRL, ODC1,
   ALDH18A1, CKMT1B, AGMAT, CPS1, GLS2, GOT2, ABP1, PYCR1, CKMT1A, P4HA1,
   P4HA3, ALDH4A1, GAMT, NOS2, SMS 53 1,897520815
   hsa04310:Wnt signaling pathway 35 0,0189717 WNT5A, WNT16, WNT5B, MMP7,
   PRKX, WNT4, CSNK2A1, WNT3, RAC3, CACYBP, CAMK2B, FZD9, WNT10A, TBL1XR1,
   DVL3, WNT10B, VANGL1, VANGL2, CHP2, LEF1, FZD3, CSNK2A1P, PORCN, FZD7,
   WNT2B, FZD6, DKK4, SENP2, WNT7B, FZD10, DKK1, SFRP1, SFRP2, WNT11,
   RUVBL1, TBL1X 151 1,456912712
   hsa04916: Melanogenesis 25 0,01931079 WNT5A, HRAS, WNT16, ADCY2, WNT5B,
   POMC, PRKX, WNT4, WNT3, MC1R, CALML3, CAMK2B, CALML5, TUBB3, FZD9,
   DVL3, WNT10A, WNT10B, LEF1, FZD3, FZD7, FZD6, WNT2B, WNT7B, FZD10,
   WNT11 99 1,587256995
   hsa05219:Bladder cancer 13 0,02646666 E2F1, EGFR, E2F2, E2F3, HRAS,
   FGFR3, PGF, MMP9, CDH1, CDK4, MMP1, CDKN2A, EGF 42 1,945523574
   hsa00670:One carbon pool by folate 7 0,0303469 TYMS, MTHFD2, SHMT2,
   ALDH1L1, DHFR, ATIC, GART 16 2,749922744
   hsa05222:Small cell lung cancer 21 0,03694956 E2F1, E2F2, CKS1B, E2F3,
   CYCS, SKP2, ITGA2, CDK6, CDK4, CDK2, COL4A6, CCNE2, LAMA1, CCNE1,
   LAMB3, ITGA6, PIAS3, LAMC2, NOS2, TRAF4, PIK3R2 84 1,571384425
   [67]Open in a new tab

   The results revealed that upregulated consistent DEGs indicated the
   genes that were associated with the top three pathways: “Cell cycle”,
   “DNA replication” and, “p53 signaling pathway”.

Table 2. Functional annotation analysis of upregulated differentially
expressed genes (DEGs) in lung adenocarcinoma (LUAD) datasets using the DAVID
tool.

   Term Count P Value Genes Pop Hits Fold Enrichment
   hsa04110:Cell cycle 43 2,82E-09 E2F1, E2F2, E2F3, E2F5, DBF4, PRKDC,
   TTK, PKMYT1, CHEK1, CHEK2, SFN, PTTG1, CCNE2, CCNE1, CDC45, MCM7,
   CDKN2A, ORC6L, BUB1, CCNA2, CDC7, CDC6, CDK1, RBL1, SKP2, ESPL1, CDC20,
   MCM2, CDC25C, MCM3, CDK4, MCM4, ORC1L, CDC25A, MCM6, CCNB1, MAD2L1,
   CCNB2, PLK1, PCNA, BUB1B, MAD2L2, SMC1B 125 2,591466667
   hsa03030:DNA replication 17 4,23E-06 LIG1, POLE, MCM2, MCM3, RNASEH2A,
   MCM4, MCM6, RFC5, PRIM1, DNA2, RFC3, RFC4, MCM7, POLE2, PRIM2, PCNA,
   FEN1 36 3,557407407
   hsa04512:ECM-receptor interaction 26 5,15E-05 IBSP, TNC, COL3A1,
   ITGA11, ITGB4, COL2A1, VTN, ITGB3, CHAD, HMMR, LAMB3, ITGB8, COMP,
   COL6A3, SV2A, COL11A2, COL11A1, THBS2, SPP1, THBS4, ITGA2, COL5A2,
   COL5A1, LAMA1, COL1A2, COL1A1 84 2,331746032
   hsa00270:Cysteine and methionine metabolism 13 9,00E-04 DNMT3A, LDHB,
   LDHA, AHCY, SRM, IL4I1, CTH, GOT1, SDS, MAT1A, DNMT1, DNMT3B, CBS 34
   2,880392157
   hsa00330:Arginine and proline metabolism 17 9,64E-04 PYCRL, ODC1,
   ALDH18A1, CKMT1B, SRM, AGMAT, CPS1, PYCR1, ABP1, CKMT1A, GOT1, CKM,
   ALDH1B1, P4HA1, P4HA3, GAMT, OAT, ADC 53 2,416352201
   hsa00250:Alanine, aspartate and glutamate metabolism 12 0,00139536
   ADSSL1, GOT1, ACY3, GFPT1, IL4I1, GPT, CAD, ASNS, CPS1, GAD1, GPT2,
   PPAT 31 2,916129032
   hsa00010:Glycolysis / Gluconeogenesis 18 0,00148413 ALDOA, LDHB, LDHA,
   ALDOB, PFKP, ADH1C, ALDH3B2, PGAM2, PCK1, ALDH3A1, GPI, TPI1, G6PC,
   ALDH1B1, PKM2, ENO3, GAPDH, ENO1 60 2,26
   hsa00601:Glycosphingolipid biosynthesis 10 0,0033785 B4GALT3, B4GALT2,
   FUT9, B3GALT5, B3GNT4, FUT6, FUT3, B3GNT3, FUT2, B4GALT4 25 3,013333333
   hsa04115:p53 signaling pathway 18 0,00619724 CDK1, LRDD, CHEK1, CHEK2,
   PMAIP1, SFN, CDK4, GTSE1, CCNB1, CCNE2, CCNE1, CCNB2, CDKN2A, SERPINB5,
   RRM2, BAI1, PERP, IGFBP3 68 1,994117647
   hsa00051:Fructose and mannose metabolism 11 0,01033738 ALDOA, TPI1,
   AKR1B15, SORD, GMDS, AKR1B10, GMPPA, ALDOB, PFKP, TSTA3, MTMR7, PMM2 34
   2,437254902
   hsa04950:Maturity onset diabetes of the young 9 0,0124429 HNF1A, HNF4A,
   BHLHA15, FOXA3, MNX1, NEUROD1, IGF2, HNF4G, PDX1 25 2,712
   hsa00260:Glycine, serine and threonine metabolism 10 0,01604608 CTH,
   SHMT2, SDS, GCAT, DMGDH, GAMT, PSAT1, PSPH, CBS, GLDC 31 2,430107527
   hsa03440:Homologous recombination 9 0,02497133 XRCC3, XRCC2, BLM, EME1,
   SHFM1, BRCA2, RAD54B, RAD54L, RAD51 28 2,421428571
   hsa00533:Keratan sulfate biosynthesis 6 0,02885381 B4GALT3, B4GALT2,
   FUT8, CHST6, CHST4, B4GALT4 14 3,228571429
   hsa00140:Steroid hormone biosynthesis 12 0,0343766 AKR1C3, AKR1C2,
   UGT1A6, UGT1A10, UGT1A9, AKR1C4, HSD17B2, HSD17B1, CYP11B1, SRD5A3,
   UGT2B4, SRD5A1, UGT2B15, AKR1C1 46 1,965217391
   hsa00512:O-Glycan biosynthesis 9 0,03707018 GALNT3, GCNT3, GALNT7,
   GALNT6, GALNT4, B3GNT6, C1GALT1, GALNT14, ST6GALNAC1 30 2,26
   hsa05219:Bladder cancer 11 0,04377571 E2F1, E2F2, E2F3, CDKN2A, PGF,
   MMP9, ERBB2, CDH1, EGF, CDK4, MMP1 42 1,973015873
   hsa00830:Retinol metabolism 13 0,04676524 BCMO1, CYP2C18, CYP2B6,
   ADH1C, RDH5, UGT1A10, UGT1A6, UGT1A9, RDH10, LRAT, DGAT1, DGAT2,
   CYP2A6, UGT2B4, UGT2B15 54 1,813580247
   hsa00980:Metabolism of xenobiotics by cytochrome P450 14 0,04737766
   CYP2C18, CYP2B6, ADH1C, ALDH3B2, CYP2E1, ALDH3A1, AKR1C3, UGT1A10,
   UGT1A6, AKR1C2, UGT1A9, AKR1C4, UGT2B4, UGT2B15, AKR1C1, MGST1 60
   1,757777778
   [68]Open in a new tab

   The results revealed that upregulated consistent DEGs indicated the
   genes that were associated with the top three pathways: “Cell cycle”,
   “DNA replication” and, “ECM-receptor interaction”.

Table 3. Functional annotation analysis of downregulated differentially
expressed genes (DEGs) in Lung squamous cell carcinoma (LUSC) datasets using
the DAVID tool.

   Term Count P Value Genes Pop Hits Fold Enrichment
   hsa04610:Complement and coagulation cascades 31 1,15E-10 C3AR1, C7,
   A2M, MASP1, C3, F13A1, C6, C5, C1QC, FGG, SERPINA1, C2, CFI, CFD, F11,
   CR1, F10, C5AR1, C4A, F8, SERPING1, C4BPA, C1QA, C1QB, VWF, C8B, CD55,
   TFPI, SERPIND1, CPB2, PROS1 69 3,586444611
   hsa04514:Cell adhesion molecules (CAMs) 44 8,51E-10 HLA-DQB1, ITGAL,
   CADM3, CLDN18, CADM1, HLA-DRB1, CLDN5, ITGB2, HLA-DMB, SDC4, HLA-DMA,
   CDH5, ITGAM, CD22, HLA-DRB5, ESAM, CD4, HLA-DPB1, HLA-DOA, SELPLG,
   NEGR1, SPN, ICAM1, SELP, PTPRC, PTPRM, ICAM2, HLA-E, HLA-DQA2, CLDN23,
   HLA-DQA1, SIGLEC1, NCAM2, ITGA9, CD86, CD34, ITGA8, CLDN2, HLA-DPA1,
   JAM2, JAM3, SELE, CD226, HLA-DRA 132 2,660910518
   hsa04640:Hematopoietic cell lineage 30 2,22E-07 CSF3, IL1R1, HLA-DRB1,
   CSF1, MME, ANPEP, IL7R, ITGAM, FCGR1C, FCGR1A, CSF3R, CD22, HLA-DRB5,
   CD4, CSF2RA, CSF1R, IL6, CR1, CD1C, ITGA1, IL6R, IL11RA, CD1E, CD55,
   CD37, CD36, CD34, CD33, EPOR, IL3RA, HLA-DRA 86 2,784673798
   hsa05416:Viral myocarditis 25 2,37E-06 HLA-DQB1, PRF1, ITGAL, CAV1,
   HLA-DRB1, ITGB2, HLA-DMB, HLA-DMA, HLA-DRB5, HLA-DPB1, HLA-DOA, ICAM1,
   HLA-E, HLA-DQA2, HLA-DQA1, LAMA2, CD86, CD55, SGCG, MYH11, SGCD,
   HLA-DPA1, SGCA, HLA-DRA, MYH10 71 2,81082097
   hsa04060:Cytokine-cytokine receptor interaction 58 1,18E-05 ACVRL1,
   PDGFB, IL6ST, LEPR, TNFSF15, CXCR1, CXCR2, TNFSF13, TNFSF12, CXCL12,
   TGFB2, CSF3R, CSF2RB, CSF2RA, IL18RAP, LIFR, IL6R, IL11RA, TNFRSF10C,
   PPBP, TNFRSF10D, CCR2, CX3CR1, TNFSF12-TNFSF13, CSF3, CCL3, IL1R1,
   CCL2, CXCL5, CSF1, CCR1, CXCL3, CXCL2, BMPR2, IL7R, CCL4, LIF,
   TNFRSF1B, IL12RB1, CCL21, IL10RA, PLEKHO2, FIGF, CSF1R, IL18R1, IL6,
   BMP2, FLT1, FLT4, TGFBR2, EDA2R, HGF, CCL18, KDR, CCL17, CCL13, CCL14,
   CXCL16, EPOR, IL3RA 262 1,767169581
   hsa04062:Chemokine signaling pathway 45 1,61E-05 ADCY4, PRKCZ, CCL3,
   CCL2, GNAI2, CXCL5, FGR, PREX1, CXCL3, CCR1, CXCL2, NCF1C, NFKBIA,
   CXCR1, GNG11, CXCR2, CCL4, CXCL12, DOCK2, CCL21, RASGRP2, GNG2, PIK3R5,
   SHC3, PLCB2, GNG7, PIK3CG, ITK, NCF1, HCK, WAS, CCL18, ELMO1, PRKCB,
   CCL17, GNGT2, CCL13, CCL14, PPBP, ADCY9, ARRB2, ARRB1, CXCL16, CCR2,
   CX3CR1, GRK5 187 1,920978181
   hsa05310:Asthma 14 1,62E-05 FCER1A, HLA-DQB1, HLA-DRB1, HLA-DMB,
   HLA-DMA, HLA-DQA2, HLA-DQA1, HLA-DRB5, MS4A2, FCER1G, HLA-DPA1,
   HLA-DPB1, HLA-DOA, HLA-DRA 29 3,853732474
   hsa04270:Vascular smooth muscle contraction 31 2,82E-05 ADCY4, ADORA2A,
   PPP1R12B, MRVI1, PRKG1, MYL9, EDNRA, AGTR1, ACTG2, PTGIR, GUCY1A2,
   PLA2G1B, GUCY1A3, CALCRL, PLCB2, PPP1R14A, RAMP3, RAMP2, PLA2G10, NPR1,
   PRKCE, ITPR1, PRKCB, PRKCQ, ADCY9, GNAQ, MYH11, CACNA1C, CACNA1D,
   PLA2G5, MYLK 112 2,209506055
   hsa05332:Graft-versus-host disease 16 3,38E-05 HLA-DQB1, PRF1, IL6,
   HLA-DRB1, HLA-DMB, HLA-E, HLA-DMA, HLA-DQA2, HLA-DQA1, CD86, HLA-DRB5,
   HLA-DPA1, HLA-DPB1, HLA-DOA, KLRD1, HLA-DRA 39 3,274966791
   hsa05322:Systemic lupus erythematosus 27 1,37E-04 HLA-DQB1, C7,
   HLA-DRB1, C3, C6, C5, HLA-DMB, C1QC, HLA-DMA, FCGR1C, FCGR1A, HLA-DRB5,
   HLA-DPB1, C2, FCGR3A, HLA-DOA, FCGR3B, C4A, HLA-DQA2, HLA-DQA1, C1QA,
   C8B, C1QB, CD86, HLA-DPA1, FCGR2A, CTSG, HLA-DRA 99 2,177108606
   hsa04672:Intestinal immune network for IgA production 17 1,84E-04
   HLA-DQB1, IL6, HLA-DRB1, TNFSF13, TNFSF12, PIGR, HLA-DMB, HLA-DMA,
   CXCL12, HLA-DQA2, HLA-DQA1, TGFB2, CD86, HLA-DRB5, HLA-DPA1, HLA-DPB1,
   TNFSF12-TNFSF13, HLA-DOA, HLA-DRA 49 2,769519111
   hsa05330:Allograft rejection 14 2,39E-04 HLA-DQB1, PRF1, HLA-DRB1,
   HLA-DMB, HLA-E, HLA-DMA, HLA-DQA2, HLA-DQA1, CD86, HLA-DRB5, HLA-DPA1,
   HLA-DPB1, HLA-DOA, HLA-DRA 36 3,104395604
   hsa04940:Type I diabetes mellitus 15 3,63E-04 HLA-DQB1, PRF1, HLA-DRB1,
   PTPRN2, HLA-DMB, HLA-E, HLA-DMA, HLA-DQA2, HLA-DQA1, CD86, HLA-DRB5,
   HLA-DPA1, HLA-DPB1, HLA-DOA, HLA-DRA 42 2,850975555
   hsa05414:Dilated cardiomyopathy 24 7,06E-04 ADCY4, SLC8A1, TNNC1,
   ITGA1, CACNG4, ITGA10, CACNB4, TTN, CACNA2D2, TGFB2, LAMA2, ITGA9, DES,
   ADRB1, SGCG, ADCY9, PLN, ITGA8, ITGA7, RYR2, SGCD, CACNA1C, CACNA1D,
   SGCA 92 2,08245171
   hsa05410:Hypertrophic cardiomyopathy (HCM) 22 0,00139878 SLC8A1, IL6,
   TNNC1, ITGA1, CACNG4, ITGA10, CACNB4, TTN, CACNA2D2, TGFB2, LAMA2,
   ITGA9, ACE, DES, SGCG, ITGA8, ITGA7, RYR2, SGCD, CACNA1C, CACNA1D, SGCA
   85 2,066118755
   hsa04614:Renin-angiotensin system 8 0,00290177 AGTR1, ACE, AGTR2, MME,
   CPA3, ANPEP, ENPEP, CTSG 17 3,756579555
   hsa05320:Autoimmune thyroid disease 14 0,0083188 HLA-DQB1, PRF1,
   HLA-DRB1, HLA-DMB, HLA-E, HLA-DMA, HLA-DQA2, HLA-DQA1, CD86, HLA-DRB5,
   HLA-DPA1, HLA-DPB1, HLA-DOA, HLA-DRA 51 2,191338074
   hsa04670:Leukocyte transendothelial migration 25 0,00994356 ITGAL,
   CLDN18, GNAI2, CLDN5, NCF1C, ITGB2, CXCL12, CDH5, ITGAM, MYL9, ESAM,
   PIK3R5, RAPGEF4, RAPGEF3, PIK3CG, ICAM1, ITK, NCF2, NCF1, NCF4, CLDN23,
   PRKCB, CYBB, CLDN2, JAM2, JAM3 118 1,691256685
   hsa04612:Antigen processing and presentation 19 0,01258302 HLA-DQB1,
   CIITA, HLA-DRB1, IFI30, CTSS, HLA-DMB, HLA-E, HLA-DMA, HLA-DQA2,
   HLA-DQA1, CD74, B2M, HLA-DRB5, CD4, HLA-DPA1, HLA-DPB1, HLA-DOA, KLRD1,
   HLA-DRA 83 1,827372283
   hsa05412:Arrhythmogenic right ventricular cardiomyopathy (ARVC) 17
   0,02357362 SLC8A1, ITGA1, CACNG4, ITGA10, CACNB4, CACNA2D2, LAMA2,
   ITGA9, DES, SGCG, ITGA8, ITGA7, RYR2, SGCD, CACNA1C, CACNA1D, SGCA 76
   1,785611006
   hsa04530:Tight junction 26 0,02458398 PRKCZ, CLDN18, GNAI2, MRAS,
   CLDN5, MYL9, RRAS, PARD6B, MAGI3, MAGI1, MPDZ, PRKCE, CLDN23, PRKCB,
   EPB41L2, EPB41L3, PRKCQ, TJP1, CGN, MYH11, CLDN2, TJP3, JAM2, JAM3,
   MYH10, SPTAN1 134 1,548888212
   hsa05020:Prion diseases 10 0,02477695 C1QA, EGR1, C1QB, NCAM2, C8B, C7,
   IL6, C6, C5, C1QC 35 2,280780444
   hsa04960:Aldosterone-regulated sodium reabsorption 11 0,0263487 PIK3CG,
   ATP1B2, HSD11B1, NR3C2, PIK3R5, NEDD4L, ATP1A2, SCNN1G, SCNN1B,
   SLC9A3R2, PRKCB 41 2,141708466
   hsa03320:PPAR signaling pathway 15 0,04364393 ACOX2, LPL, OLR1, PPARG,
   RXRG, ACADL, ACSL1, CD36, SORBS1, CYP27A1, HMGCS2, FABP3, FABP4, ACSL4,
   ACSL5 69 1,735376425
   hsa04666:Fc gamma R-mediated phagocytosis 19 0,04492025 PIK3CG, DNM3,
   PTPRC, WASF3, NCF1, HCK, NCF1C, PIP5K1B, PRKCE, WAS, PRKCB, DOCK2,
   GAB2, FCGR1C, CFL2, FCGR1A, PLA2G4F, PIK3R5, FCGR2A, FCGR3A, PPAP2B 95
   1,596546311
   hsa04020:Calcium signaling pathway 31 0,04546262 GNA14, ADCY4, ERBB4,
   CYSLTR1, ADORA2A, TNNC1, EDNRA, AGTR1, EDNRB, PDE1B, PLCB2, SLC8A1,
   BST1, PTGFR, ITPR1, PRKCB, ADRB2, P2RX7, PLCE1, ADRB1, P2RX1, ADCY9,
   GNAQ, PLN, TBXA2R, CACNA1H, RYR2, CACNA1C, CACNA1D, MYLK, PTAFR 176
   1,406049308
   [69]Open in a new tab

   The results revealed that downregulated consistent DEGs indicated the
   genes that were associated with the top three pathways: “Complement and
   coagulation cascades”, “Cell adhesion molecules (CAMs)” and,
   “Hematopoietic cell lineage”.

   As Schwann cells act on both of the aforementioned upregulated
   pathways, we selected datasets from the GEO database containing gene
   expression profiles of both early and late passage human Schwann cells
   exposed to cancer growth factors heregulin and forskolin ([70]GSE4030).
   These expression profiles were used to identify DEGs with the aid of
   the online tool GEO2R. We found 638 upregulated and 1250 downregulated
   genes. Then, we performed a GO term enrichment and functional
   classification by KEGG analysis, using DAVID platform, to investigate
   the biological and functional roles of these DEGs. Upregulated DEG
   enrichment included “neuroactive ligand-receptor interaction,” while
   “pathways in cancer” and “focal adhesion” were enriched for
   downregulated DEGs ([71]Supplementary Tables 1-2).

   To facilitate the analysis of the large throughput of DEGs, a protein
   classification analysis was performed using the PANTHER classification
   system. According to the study, upregulated genes categories were
   identified, in which hydrolase (PC00121), enzyme modulator (PC00095),
   nucleic acid binding (PC00171), receptor (PC00197), transcription
   factor (PC00218), transporter (PC00227), and signaling molecule
   (PC00207) are the top seven most abundant protein classes. However,
   protease (PC00190) is the top one most abundant protein class among
   downregulated DEGs ([72]Figure 1A-B).

Figure 1. Molecular function and protein class terms of downregulated.

   [73]Figure 1
   [74]Open in a new tab

   (A) and upregulated (B) among the differentially regulated genes in
   Schwann cells. Figure 1C-D - Venn Diagrams of combined overrepresented
   differentially expressed genes in Schwann cells and lung adenocarcinoma
   (C and D), and Lung squamous cell carcinoma (E and F) cancer samples.
   The results are showing overlapping downregulated (Figure 1C) and
   upregulated (Figure 1D) genes between lung adenocarcinoma and Schwan
   cells. Overlapped downregulated and upregulated genes between Lung
   squamous cell carcinoma and Schwan cells are showed in Figure 1E and
   Figure 1F, respectively. Figure 1G-H - Venn Diagrams of downregulated
   miRNA target genes combined with upregulated genes in lung
   adenocarcinoma (G) and Lung squamous cell carcinoma (H). Figure 1I-J -
   Venn Diagrams of upregulated miRNA target genes combined with
   downregulated genes in lung adenocarcinoma (I) and Lung squamous cell
   carcinoma (J). The number in each intersecting region represents the
   number of overlapping genes.

Overview of the cancer transcriptomic analysis

   We conducted a systematic and integrative analysis to explore cancer
   type-specific and Schwann cell specific DEGs, to construct a cancer
   network. We first determined DEGs by comparing gene expression levels
   between tumor and normal samples. A Venn diagram was then constructed
   to visualize the overlap between DEG genes, both upregulated and
   downregulated, from both cancer types and from Schwann cells. Results
   are shown in [75]Figure 1C-F. Among upregulated genes from LUSC and
   LUAD, we identified 84 and 70 co-expressed DEGs with Schwann cells,
   respectively.

   We also ran a KEGG analysis to investigate pathways with major
   expression changes in studied cell lines and lung cancer, based on
   previously overlapping identified DEGs. The analysis revealed that
   “neuroactive ligand-receptor interaction” and “p53 signaling pathway”
   were among the most active pathways according to overlapped Schwann
   cells and LUSC DEGs. In contrast, overlapped Schwann cells and LUAD
   DEGs samples show “neuroactive ligand-receptor interaction” pathway as
   mainly downregulated, while “p53 signaling pathway” is upregulated
   ([76]Supplementary Tables 3-6). Grm1, Chrm3, Gabra5, Gpr50, and Pth2r
   were the common upregulated genes associated with the “neuroactive
   ligand-receptor interaction” pathway; Ccne1, Cdkn2a, and Perp were
   commonly associated to “p53 signaling pathway”.

   To confirm our results, we performed a differentially expressed miRNA
   analysis over cancer data. Data from identified miRNAs target genes
   were crossed with DEGs data from Schwann cells ([77]Figure 1G-J). We
   found that downregulated miRNAs had target genes associated with
   “neuroactive ligand-receptor interaction” in LUAD and “pathways in
   cancer” in both cancers ([78]Supplementary Tables 7-8). In contrast,
   upregulated miRNAs had target genes associated with the “axon guidance”
   pathway in both lung cancers ([79]Supplementary Tables 7-10). Genes
   from the “axon guidance” pathway, common to all our analyses, were
   Robo2 and Slit2 genes. ROBO2 and SLIT2 proteins were highly expressed
   in normal tissue compared to cancer samples (data not shown).

Hallmarks of cancer analysis

   In order to understand the mechanism by which Schwann cells aid in
   neoplastic development, we analyzed the behavior of genes associated
   with cell differentiation processes (Gata3, Cdh1, and Cdh2), apoptosis
   (Casp3, Casp9, Bax, and Bcl2), motility (Cxcr2, Cxcl5, Mmp9, and
   Ccl12), and cell proliferation (Mki67). While Cdh1 and Cdh2 showed an
   increased mRNA expression in LUSC samples, only CDH2 protein expression
   was decreased in LUAD when compared with LUSC (data not shown). Also,
   Mmp9 and Mki67 were overexpressed in both lung cancers
   ([80]Supplementary Figures 1-2).

Analysis of genes involved in dedifferentiation of Schwann cells

   Schwann cells produce cell differentiation maintenance proteins (Sox10,
   S100, Egr2, Mbp, and Mpz) that, after nerve damage, diminish their
   expression, provoking cellular dedifferentiation. Consequently, the
   expression of a new set of proteins that form non-myelinated cells,
   Sox10, Gap43, S100, Ncam1, Ngfr1, and Gfap, is augmented. We observed
   increased expression of Gap43 and Gfap genes in both lung cancers,
   whereas Ngfr1 was upregulated only in LUSC cells ([81]Supplementary
   Figures 1-2).

Methylation and protein analyses

   Methylation analysis of the Cdh1, Cdh2, Mmp9, Mki67, Gfap, Gap43,
   Robo2, and Slit2 genes from tumor samples demonstrated that all of them
   were methylated in their promoter regions unlike those from normal
   tissues. However, there was no significant correlation between
   methylation and Gfap gene expression in lung cancers. Whereas there was
   a positive correlation between methylation and Gap43 gene expression in
   LUSC samples ([82]Supplementary Figures 3-24).

Copy number alteration

   Copy number alteration data demonstrated that Cdh1, Cdh2, Gfap, Perp,
   and Robo2 had a higher mRNA expression than normal tissues; increased
   expression was associated with gain or amplification alterations in
   LUAD samples ([83]Supplementary Figure 25). Similarly, in LUSC samples,
   Ccne1, Cdh1, Cdkn2a, Gap43, and Perp had a higher mRNA expression
   associated with gain or amplification alterations ([84]Figure 2).

Figure 2. Correlation of copy number variation and expression of.

   [85]Figure 2
   [86]Open in a new tab

   Gap43 (A), Gfap (B), Robo2 (C) and Slit2 (D) genes in lung
   adenocarcinom a (LUAD) and lung squamous cell carcinoma (LUSC). Figures
   were generated using Cbioportal data.

Schwann cell differentiation protein expression in lung cancer samples

   For analysis of proteins expressed in dedifferentiated Schwann cells,
   we initially identified the gene list in Pubmed and GeneCards
   databases. We only analyzed genes with relevance score higher than
   seven. The relevance scores of elements related to genes are based on
   the analysis of co-occurrences of two elements in Medline documents.
   The observed are compared to an expected value based on a
   hypergeometric distribution. We identified 325 Schwann cell
   dedifferentiated related genes in both databases. Data were then
   cross-checked with previously published protein data expression
   analysis [[87]9], which the expression of normal lung tissue and lung
   cancer were compared. 10 proteins (GFAP, STAT3, SRC, CD36, CAV1, PCNA,
   HDAC9, AQP1, APOA1, RALA) associated with dedifferentiation of Schwann
   cells were increased in lung cancer, including GFAP. No downregulated
   protein expression was associated with Schwann cell dedifferentiation.

Cancer protein expression patterns correspond to pathway activation levels

   We also performed an RPPA protein analysis. Only CDH1 and CDH2 protein
   expression data were available for analysis. We observed that CDH2 was
   overexpressed in LUSC compared to LUAD. Additionally, no significant
   difference was found in CDH1 analysis ([88]Supplementary Figure 26).

   In order to analyze the pathway by which Schwann cells induce
   neoplastic and their own cell proliferation and migration, we evaluated
   the expression of MEK1 (MEK1 and MEK1_pS217S221), ERK2, AKT
   (PRAS40_pT246, AKT_pT308 and AKT_pS473), RAF (CRAF and CRAF_pS338), and
   GSK3Β (GSK3_pS9 and GSK3ALPHABETA_pS21S9) proteins. We observed higher
   expression of ERK2, AKT, CRAF and GSK3Β proteins in LUSC when compared
   with LUAD. Only MEK1_pS217S221 presented higher expression in LUAD
   ([89]Supplementary Figure 26).

Prognostic analyses

   TIMER survival analyses showed that Ccne1, Mki67, and Gap43 mRNA higher
   expression are associated with poor LUAD prognostic ([90]Supplementary
   Figure 27). Similarly, PRECOG analysis showed that Grm1, Gap43 and
   Mki67 higher mRNA expression is associated with poor prognostic in
   LUAD. However, PRECOG analysis demonstrated Slit2 and Robo2
   downregulation are associated with poor survival in LUAD samples
   ([91]Figure 3). [92]Figure 4 illustrates the possible association
   between Schwann cells and lung cancer.

Figure 3. PRECOG Survival analyses in lung adenocarcinoma for Grm1.

   [93]Figure 3
   [94]Open in a new tab

   (A), Gap43 (B), Slit2 (C) and Robo2 (D) mRNA expression. Note that Grm1
   and Gap43 upregulation and, Slit2 and Robo2 downregulation are
   associated with poor survival.

Figure 4. Mechanism of action of Schwann cells in the development of cancer.

   [95]Figure 4
   [96]Open in a new tab

   We suppose that Schwann cells assist the neoplastic lung cells to
   invade surrounding tissue during the early stages of carcinogenesis.
   This mechanism would be possible through the dedifferentiation of the
   Schwann cells and blocking axon guidance signaling pathways. Migration
   of Schwann cells to the peritumoral region would be associated with
   perineural invasion and, subsequently, central nervous system
   metastases.

DISCUSSION

   Cells surrounding a tumor form a molecular microenvironment known as
   stroma. Stroma formation is influenced by tumor cells and, in turn, it
   influences tumor growth, migration, and invasion [[97]17]. Depending on
   the characteristics of the primary tumor, the stroma, and the intrinsic
   ability of metastatic tumor cells to adapt to a new location, malignant
   cells use distinct mechanisms for proliferation, survival, and
   dissemination. These cells often reactivate the expression of genes
   employed during embryogenesis [[98]17, [99]18]. In order to leave the
   primary tumor and disseminate to distant organs, metastatic cells lose
   the ability to adhere to adjacent cells, enhancing their migratory and
   invasive capacity. This mechanism is accompanied by several
   modifications in the expression of genes, such as loss of epithelial
   receptor expression and increased expression of mesenchymal markers, a
   phenomenon also known as epithelial-mesenchymal transition [[100]18,
   [101]19].

   We investigated the association between Schwann cells and those lung
   cancer types often associated with perineural invasion. Initially, we
   used the GEO DataSets platform from the GEO repository to identify a
   database reporting gene expression in Schwann cells in a neoplastic
   context. Briefly, the database contains the expression results from
   experiments in which two factors produced by tumor cells were added
   into cell cultures. Comparisons were made between samples from the
   first and third passages. We then used these data to perform
   differential gene expression analysis and crossed data from upregulated
   genes with differential expression data from LUAD and LUSC. After
   identifying the genes in common, we ran several analysis tools to
   identify molecular pathways associated to these genes. Interestingly,
   we noted that the “neuroactive ligand-receptor interaction” pathway was
   active in LUSC. Studies have demonstrated the association between
   neuroactive ligand-receptors in the control of Schwann cell
   differentiation as well as in neoplastic cell proliferation. For
   example, GRM1 can activate the RAS/MEK1/ERK/AKT/RAF/GSK3Β
   phosphorylation cascade [[102]20, [103]21], leading to an increase in
   cell proliferation and invasion. Grm1 was upregulated in LUSC,
   suggesting that this gene participates in lung cancer development.

   In these context, it has also been demonstrated that neuroactive
   ligand-receptors are active during Wallerian degeneration of peripheral
   nerves. Together with macrophages, Schwann cells remove axon and myelin
   debris, and clear a path for subsequent axonal regrowth and nerve
   regeneration [[104]22]. Tumor cells benefit from nerve regeneration
   machinery to promote cell proliferation, migration, and invasion. It
   has been reported that Schwann cells are induced to migrate to the
   region close to the tumor at the beginning of the carcinogenic process.
   It was also suggested that Schwann cells promote neoplastic invasion by
   direct contact with cancer cells, since paracrine signaling and matrix
   remodeling are not yet sufficient to induce the migration process
   [[105]4, [106]23]. Cell-cell contact between Schwann cells and tumor
   tissue is necessary to potentiate the ability of neoplastic cells to
   penetrate into the underlying tissue [[107]1, [108]3]. After contact,
   the degradation of the extracellular matrix by Schwann cells provokes
   the formation of tunnels or bands coated with laminin, due to Schwann
   cells’ capacity to express matrix metalloproteins, especially MMP2 and
   MMP9. The mechanism of extracellular matrix degradation that promotes
   neoplastic migration also depends on the production of NCAM1 and
   N-cadherin (CDH2) by Schwann cells [[109]1, [110]7].

   In the present study, we analyzed the expression of genes associated
   with the hallmark of cancer. The importance of analyzing these genes
   derives from the hypothesis that there may be a correlation between the
   presence of Schwann cells and the aggressiveness of the neoplastic
   cell. Increased expression of genes related to cell differentiation
   (Cdh1 and Cdh2), motility (Mmp9), and proliferation (Mki67), mainly in
   LUSC, suggests this correlation. We observed that these genes are
   mutated and had higher expression in lung cancer samples when compared
   with normal tissue.

   Decreased E-cadherin protein expression after contact of Schwann cells
   with the tumor resembles the mechanism followed by cells during axonal
   repair process [[111]8]. The loss of E-cadherin during the
   epithelial-mesenchymal transition in cancer is associated with a
   positive regulation of NCAM1 and CDH2 [[112]1, [113]3, [114]4, [115]8].
   When E-cadherin is suppressed, NCAM1 and CDH2 are upregulated; they
   associate with the p59fyn protein, whose subsequent activation leads to
   inhibition of focal adhesion and an increase in cell migration. A study
   using oncogenic K-ras pancreatic cancer cell lines identified increased
   levels of polysialylated NCAM1 expression, which interacts with
   E-cadherin to create steric hindrance of homophilic binding and
   decrease cell adhesion [[116]1]. In our study, we observed upregulation
   of E-cadherin and CDH2 in both lung cancers. A CDH2 protein RPPA
   analysis showed that it is overexpressed in LUSC. Copy number
   alteration analysis also showed that amplification and gain are
   associated with E-cadherin and N-cadherin mRNA levels in both cancers.
   However, we showed that Cdh2 mRNA expression is higher in LUSC when
   compared to LUAD. No difference was observed neither in Ncam1 mRNA nor
   in E-cadherin protein analyses.

   Differentiation of myelinating Schwann cells may undergo interference
   from inhibitory pathways that negatively control the expression of
   genes responsible for myelin sheath formation. NOTCH1 and JUN stand out
   among the negative regulators of the myelination program, in the same
   way as SOX-2 and PAX-3 [[117]2, [118]24, [119]25]. The myelinating
   phenotype also involves the inactivation of a number of genes linked to
   the production of immature Schwann cell markers. Some transcription
   factors are responsible for ensuring proper maturation in Schwann
   cells. SOX-10 acts synergistically with a second factor, OCT-6,
   resulting in the expression of Krox-20. In turn, KROX-20 is a key
   inducer of expression of myelin genes, such as Mbp, Mpz and Prx. The
   maintenance of the myelinating phenotype therefore requires the
   continuous expression of KROX-20 and SOX-10, considering that
   inactivation of both proteins results in dedifferentiation of Schwann
   cells [[120]24].

   Previous studies have shown that Schwann cells induce cellular
   aggressiveness in lung cancer [[121]26]. During dedifferentiation,
   Schwann cells express proteins initially lost during the myelination
   process. Among these proteins are GAP43, NCAM1, P75NTR (Ngfr), GFAP and
   SOX-2 [[122]2]. Data from our study showed that Gfap, Ngfr and Gap43
   gene expression is increased in lung cancers, especially in LUSC. To
   confirm our data, we performed an analysis by using protein expression
   data published before [[123]9]. We selected only genes associated with
   Schwann cells dedifferentiation. The results showed that proteins
   associated with dedifferentiated Schwann cells are overexpressed in
   lung cancer, when compared with normal lung samples. Besides, Gfap and
   Gap43 copy number alterations are associated with higher mRNA levels in
   LUAD and LUSC, respectively. It is likely that Schwann cells in these
   tissues are dedifferentiated, aiding tumors in their mechanisms of cell
   proliferation, migration, and tissue invasion.

   During the process of carcinogenesis, there is an inactivation of the
   Slit2 and Robo2 genes, being therefore considered as tumor suppressor
   genes [[124]27]. Both genes are extremely important during the process
   of nerve formation and repair. Both SLIT2 and ROBO2 have been reported
   to inhibit migration of Schwann cells [[125]28]. In the present study,
   we observed that Slit2 and Robo2 mRNAs are decreased in lung cancer
   samples when compared to normal tissue. Also, there is a participation
   of miRNAs and of methylation of their promoter regions in the
   regulation of these genes. Similarly, both genes are targets for
   upregulated miRNAs. Immunohistochemical expression analyses revealed
   low expression of both genes in the lung. Thus, we suggest that
   decreased expression of Slit2 and Robo2 genes in lung cancer may favor
   the migration of Schwann cells; consequently, favoring invasion by
   neoplastic tissue.

   We next analyzed the RAS/ MEK1/ ERK/ AKT/ RAF/ GSK3Β phosphorylation
   cascade. This pathway is activated during the proliferation process of
   Schwann cells and lung cancer. The authors observed that SC-treated
   lung cancer cells exhibit an increased phosphorylation of Akt and
   GSK-3β. Besides, they blocked CXCL5 expression in SC significantly
   decreases the phosphorylation of AKT and GSK-3β in SC-treated lung
   cancer cells; PI3K inhibitor LY294002 blocks the activation of
   AKT/GSK-3β signaling in SC-conditioned lung cancer cells. In our study,
   we observed that the RAS/ MEK1/ ERK/ AKT/ RAF/ GSK3Β phosphorylation is
   increased in LUSC when compared to LUAD, suggesting that, in LUSC
   samples, this pathway could be more active.

   In summary, we observed that the “neuroactive ligand-receptor
   interaction” pathway was upregulated in lung squamous cell carcinoma
   and downregulated in lung adenocarcinoma. The “p53 signaling pathway”
   was active in both lung cancers, since Ccne1, Cdkn2a, and Perp were
   upregulated. Meanwhile, upregulated miRNAs inactivate the “axon
   guidance” pathway, targeting Robo2 and Slit2 genes. Both genes are also
   associated with Schwann cells migration inhibition. Also, Gfap and
   Gap43 are overexpressed, leading to Schwann cells dedifferentiation. We
   believe that Schwann cells’ dedifferentiation and proliferation are
   induced by neoplastic tissue; consequently, Schwann cells produce
   different factors that will participate in several processes of tumor
   progression. These processes may also be involved in tumor invasion
   into the perineural tissue in lung cancer. Our data stimulate efforts
   in news studies to achieve the experimental and clinical validation
   about these biomolecules.

MATERIALS AND METHODS

Collection and inclusion criteria of studies

   We searched the GEO database ([126]https://www.ncbi.nlm.nih.gov/geo/)
   for publicly available studies. After a systematic review, only one GSE
   studies were retrieved. The inclusion criteria for studies were as
   follows: (1) experiments involving both cancer and Schwann cells and,
   (2) gene expression profiling of mRNA. Then, one gene expression
   profile ([127]GSE4030) was collected for our analysis.

Microarray data and Data processing

   One gene expression profiles ([128]GSE4030) was downloaded from the GEO
   database. GEO2R was applied to compare gene expression profiles in
   early and late passage human Schwann cells exposed to the cancer growth
   factors heregulin and forskolin. Because both groups in this comparison
   have been exposed to mitogens, differences in gene expression profiles
   will be interpreted as indicative of changes caused by prolonged versus
   short term exposure to mitogens.

   GEO2R ([129]http://www.ncbi.nlm.nih.gov/geo/geo2r/) is an interactive
   web tool for comparing two groups of data that can analyze any GEO
   series. The adjusted p-values using Benjamini and Hochberg false
   discovery rate method by default were applied to correct the occurrence
   of false positive results. An adj. P < 0.05 and a logFC ≥ 1 were set as
   the cut-off criteria.

Gene list collection

   We used Entrez Gene from NCBI (www.ncbi.nlm.nih.gov/gene/) and
   GeneCards ([130]https://www.genecards.org/), which were used as the
   identifiers for Schwann cell dedifferentiation-related genes. The gene
   list and downregulated/ upregulated proteins of a previous published
   data [[131]9] were combined and identified with a Venn Diagram 2.1.0
   ([132]http://bioinfogp.cnb.csic.es/tools/venny/index.html).

Functional and pathway enrichment analysis

   Gene ontology (GO) analysis of the relevant biological processes,
   cellular components, and molecular functions was performed using the
   protein analysis through evolutionary relationships program (PANTHER,
   www.pantherdb.org), a curated database of protein families, functions,
   and pathways. GO terms assigned into identified molecules were
   classified according to their functions [[133]10].

   The Kyoto Encyclopedia of Genes and Genomes (KEGG) is an integrated
   database resource for biological interpretation of genome sequences and
   other high-throughput data [[134]11]. KEGG analyses were available at
   the DAVID database ([135]https://david.ncifcrf.gov/), a data resource
   composed of an integrated biology knowledge base and analysis tools to
   extract meaningful biological information from large quantities of
   genes and protein collections. A p-value < 0.05 was set as the cut-off
   criterion [[136]12].

RNA-seq and clinical information data from The Cancer Genome Atlas (TCGA)

   We used TCGAbiolinks, an R/Bioconductor software
   ([137]http://bioconductor.org/packages/release/bioc/
   html/TCGAbiolinks.html) [[138]13] and the interphase TCGAbiolinksGUI
   [[139]14] to download genomic and clinical data of both normal and
   solid tumor tissues for two different types of cancer from TCGA.
   Selected cancer types were lung squamous cell carcinoma (LUSC) and lung
   adenocarcinoma (LUAD). We retrieved level data for raw count mRNA and
   miRNA expression (Illumina HiSeq 2000). Co-expressed upregulated and
   downregulated DEGs from the gene expression profiles were combined and
   identified with a Venn Diagram 2.1.0
   ([140]http://bioinfogp.cnb.csic.es/tools/venny/index.html). An adj. P <
   0.05 and a logFC ≥ 1 were set as the cut-off criteria.

   DNA methylation analyses were performed using MEXPRESS dataset
   ([141]https://mexpress.be/?ref=labworm). The cancer dataset, consisting
   of DNA methylation data (Illumina Infinium Human Methylation 450 K Bead
   array, Illumina, USA) and the β value, was considered as significantly
   hypermethylated only if the value was found in more than 5% of the
   tumors [[142]15]. Copy number alteration analysis was performed using
   the cBio Cancer Genomics Portal ([143]http://cbioportal.org).

Prediction of miRNA targets

   The target miRNAs of the DEGs from TCGA data were predicted with miRDB
   ([144]http://mirdb.org/miRDB/), which is an online database for
   predicting microRNA targets [[145]16].

Lung cancer expression analyses

   To obtain individual gene-protein data, relevant information from The
   Cancer Proteome Atlas (TCPA) website
   ([146]https://tcpaportal.org/tcpa/analysis.html) was used as the
   primary source of information for reverse phase protein array (RPPA)
   analysis. Several other web resources were used as source of
   information while some more were used as analysis tools for lung cancer
   expression studies: Immunohistochemistry image-based protein data for
   both normal and cancer samples are available at the Human Protein Atlas
   ([147]https://www.proteinatlas.org/). LOCATE database for protein
   subcellular location was included on the analysis
   (http://locate.imb.uq.edu.au/cgi-bin/sort_search.cgi). Survival
   analysis of the TCGA data was performed using the Survival module of
   the Tumor Immune Estimation Resource (TIMER) and Prediction of clinical
   outcomes from genomics (PRECOG). Kaplan–Meier plots were drawn using
   TIMER to explore the association between clinical outcome and gene
   expression, and to visualize survival differences.

Statistical analysis

   Statistical analyses involving copy number variation were performed
   using GraphPad Prism 7 software. A two-way ANOVA and a two-tailed
   unpaired t test were used to compare the means between groups.

SUPPLEMENTARY MATERIALS FIGURES AND TABLES

   [148]oncotarget-10-7288-s001.pdf^ (19.7MB, pdf)

ACKNOWLEDGMENTS AND FUNDING