Schwann cell reprogramming and lung cancer progression: a meta-analysis of transcriptome data

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.


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 Meta-Analysis patient through pain and paralysis, associated with a high risk of local recurrence and decreasing survival patients rates [1][2][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 redifferentiate into subtypes with different phenotypes. After the nerve damage, Schwann cells become dedifferentiated, lose their capacity to myelinate and, promote neuronal orientation during repair [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) [1,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 [5,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 [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 [2,7,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 [2,3,7,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 overrepresentation 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 posttranscriptional regulatory network. Moreover, microarraybased 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.

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 (Tables 1-4).
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 (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 (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 ( Figure 1A-B).
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 ( 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 (Supplementary  Tables 7-8). In contrast, upregulated miRNAs had target genes associated with the "axon guidance" pathway in both lung cancers (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).

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 ( 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 ( 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 (Supplementary Figure 25). Similarly, in LUSC samples, Ccne1, Cdh1, Cdkn2a, Gap43, and Perp had a higher mRNA expression associated with gain or amplification alterations (Figure 2).

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 [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. 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". www.oncotarget.com  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).

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 (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 (Supplementary Figure 26).

Prognostic analyses
TIMER survival analyses showed that Ccne1, Mki67, and Gap43 mRNA higher expression are associated with poor LUAD prognostic (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 ( Figure  3). Figure 4 illustrates the possible association between Schwann cells and lung cancer.

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 [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 [17,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 [18,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 [20,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 [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 [4,23]. Cellcell contact between Schwann cells and tumor tissue is necessary to potentiate the ability of neoplastic cells to penetrate into the underlying tissue [1,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 [1,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 [8]. The loss of E-cadherin during the epithelialmesenchymal transition in cancer is associated with a positive regulation of NCAM1 and CDH2 [1,3,4,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 [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 [2,24,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 [24].
Previous studies have shown that Schwann cells induce cellular aggressiveness in lung cancer [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 [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 [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 [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 [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 www.oncotarget.com 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.

Collection and inclusion criteria of studies
We searched the GEO database (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 (GSE4030) was collected for our analysis.

Microarray data and Data processing
One gene expression profiles (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 (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 (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 [9] were combined and identified with a Venn Diagram 2.1.0 (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 [10].
The Kyoto Encyclopedia of Genes and Genomes (KEGG) is an integrated database resource for biological interpretation of genome sequences and other highthroughput data [11]. KEGG analyses were available at the DAVID database (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 [12].

RNA-seq and clinical information data from The Cancer Genome Atlas (TCGA)
We used TCGAbiolinks, an R/Bioconductor software (http://bioconductor.org/packages/release/bioc/ html/TCGAbiolinks.html) [13] and the interphase TCGAbiolinksGUI [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 (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 (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 [15]. Copy number alteration analysis was performed using the cBio Cancer Genomics Portal (http://cbioportal.org).

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

Lung cancer expression analyses
To obtain individual gene-protein data, relevant information from The Cancer Proteome Atlas (TCPA) website (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 (https://www.proteinatlas.org/). LOCATE database for protein subcellular location was included on the analysis (http://locate.imb.uq.edu.au/cgibin/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 twoway ANOVA and a two-tailed unpaired t test were used to compare the means between groups.

Note
LOCATE database for protein subcellular location was included in the analysis (http://locate.imb.uq.edu.au/ cgi-bin/sort_search.cgi), that was used in the study design has been removed by the owner and is no longer available.
Novelty and Significance: 1. What Is New? This is the first study showing that Schwann cells is associated with lung cancer, by alteration in mRNA, miRNA, methylation and protein expressions, deregulating several pathways such as axon guidance and neuroactive ligand receptor-interaction.