Bioinformatic Analysis Identifying S100a9 Associated With Cell Proliferation and Metastasis in Neuroblastoma

Background: More than half of neuroblastoma (NB) patients presented with distant metastases and the mortality of patients suffering from metastatic relapse was about 90%. It is urgent to nd a biomarker that can facilitate the prediction of metastasis in NB patients. Methods and Results: In the present study, we systematically analyzed Gene Expression Omnibus (GEO) datasets and focused on identifying the critical molecular networks and novel key hub genes implicated in NB. We found totally, 176 up-regulated and 19 down-regulated differentially expressed genes (DEGs) were identied. Based on these DEGs, a PPI network composed of 150 nodes and 452 interactions was established. PPI network identication combined with qRT-PCR, ELISA and IHC, S100A9 as was screened as an outstanding gene. Furthermore, in vitro tumorigenesis assays demonstrated that S100A9 overexpression enhanced the proliferation, migration, and invasion of NB cells. Conclusions: Taken together, our ndings suggested that S100A9 could participate in NB tumorigenesis and metastasis and that S100A9 has the potential to be used as a biomarker in the prediction of NB metastasis. apoptosis[33], [35],


Introduction
Neuroblastoma (NB), the most aggressive form of solid tumor of infants, account for 15% of cancer deaths in children. [1,2] For NB patients in high-risk, the main preferred treatment choices remain chemotherapy and radiotherapy, which lead to tremendous toxicity and drug resistance inevitably, , the cases ratio to recurrence and progression is about 50% [3]. Thus, it is urgent for us to identify new effective biomarker for early diagnosis, metastasis prediction and ideal therapeutic target for NB patients. S100A9, which is also called calgranulin B or migration inhibitory factor-related protein 14 , belongs to the low-molecular-weight calcium-binding S100 protein family. This family comprises more than 20 members. [4] Many evidences show S100A9 is correlated with tumorigenesis, immuno uorescence analysis has shown that S100A9 protein is elevated in metastatic melanoma and prostate cancer [5,6]. in these tumors, increased expression of S100A9 was correlated with poor differentiation. In both malignant and surrounding non-malignant tissues, a huge number of in ammatory cells expressing S100A9 were associated with signi cantly shorter cancer survival. [7]. while downregulated the expression of S100 family members is correlated with tumor proliferation, in ammation invasion, and angiogenesis,.
In the present study, based on the Geo database [14], bioinformatic analysis was further performed. After screening of DEGs and functional analysis, PPIs of DEGs were analyzed. And we found that the expression of S100A9 is pretty high in NB patients. In addition, we investigated the biological functions of S100A9 in NB cell line, as well as in primary human NB samples. These results may identify new effective biomarker for early diagnosis, metastasis prediction and ideal therapeutic target for NB patients, and provide valuable biological information for further investigation of NB.

Bioinformatic analysis
Data processing and analysis of gene expression pro le The public microarray dataset GSE90121, which was obtained based on the Affymetrix GPL570 platform Affymetrix Human Genome U133 Plus 2.0 Array), was downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/), [14]. This dataset was deposited by Kaplan D et al., containing information from human NB SK-N-AS metastatic subpopulations isolated after in vivo selection, aimed to identify genomic signatures that regulate metastasis and candidate therapeutics for NB patients. A total of 16 samples were included in the current dataset, including 12 metastatic samples and 4 primary samples. Robust Multi-array Average (RMA) affy package of Bioconductor was used to adjust the raw data. The processed gene expression data was then ltered to include those probe sets with annotations which reference the new version annotation les. To identify DEGs, we used the Linear Models (Microarray Data package in Bioconductor) to compare the expression levels of genes between the metastatic group and the localized tumor group [15].
An adjusted p-value of <0.05 and a |log2FC (fold change) | of ≥2 was used as the threshold.

Functional enrichment analysis of DEGs
Database for annotation, visualization and integrated discovery (DAVID) integrates a set of functional enrichment tools to distinguish functional genes underlying diseases processing (http://david.abcc.ncifcrf.gov/) [16]. GO (Gene Ontology) function and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analyses of DEGs were performed based on DAVID. P value < 0.05 and count ≥ 2 was regarded as statistically signi cant differences.

Protein-protein interaction network construction by STRING
We used the Search Tool for the Retrieval of Interacting Genes (STRING), an online tool and biological database for prediction of interactions between proteins, to construct the PPI network [17]. According to our analysis based on STRING, score median con dence > 0.15 was the standard of PPIs of DEGs selections. The Cytoscape software was used to visualize the PPI network. [18]. The proteins that have many interaction partners constitute the extremely important nodes in the PPI network. These proteins were de ned as the hub proteins in this study. To identify such hub proteins in the PPI network, we utilized six bioinformatic tools, namely Closeness, Degree, EPC, MNC, Radiality, and Stress centrality. Sub-network analysis was then conducted to help us discover the outstanding genes. NB patients, tissue samples and cell lines Primary tumor tissues were obtained from 9 NB patients with bone marrow metastasis and 10 NB patients without bone marrow metastasis who had undergone tumor resection surgery at the A liated Hospital of Qingdao University. None of the included patients was treated with chemotherapy, hormonal therapy or radiotherapy before the tumor resection surgery. Written informed consent was obtained from all the participants. The current research was conducted with the permission of the Medical Ethical Committee of A liated Hospital of Qingdao University (Qingdao, China).
Human NB cell line SH-SY5Y was kindly provided by Professor Xiao from the Guizhou Medical University. The cells were cultured in Dulbecco's modi ed Eagle's medium (DMEM) contained 10% fetal bovine serum (FBS, Hyclone, USA) under the conditions of 37°C, 5% CO 2 .
NB serum samples and Quantitative real-time PCR The expression of TAC1, PTGS2 and FGF1 was examined by quantitative real-time reverse transcription polymerase chain reaction (qRT-PCR). The collected tissues were immediately frozen at -80°C after surgery. Total RNA was extracted from cancer tissues by TRIzol reagent (Invitrogen Life Technologies). Then the extracted RNA was transcribed into cDNA using random primers and analyzed with an ABI 7000 Real-Time PCR System (Applied Biosystems). PCR primers were as follows: TAC1 primers: CCA CCC TGT-3′. Reactions were performed in triplicate using SYBR Green master mix (TaKaRa, Japan) and normalized to GAPDH mRNA level using the ΔΔCt method.

ELISA
To quantify levels of S100A9 in plasma, ELISA was performed as previously described [19] By using human S100A9 (JYM0539Hu, JYM, China) ELISA kits, S100A9 in plasma of the NB patients were detected according to the manufacturer's recommended procedure.

IHC staining
In brief, the formalin xed, para n-embedded tissues sections were depara nized, rehydrated and boiled in 0.01 M citrate buffer for 10 min, then incubated with 0.3% H 2 O 2 in methanol to block endogenous peroxidase activity. Then the sections were incubated with the anti-S100A9 antibody (Cell Signaling Technology, USA), followed by incubation with secondary antibody tagged with the peroxidase enzyme for 30min at room temperature, nally visualized with 0.05% DAB (3,3'-diaminobenzidine) until the desired brown reaction product was obtained. All slides were observed under a OLYMPUS BX41 Microscope, and representative photographs were taken.
Construction of plasmids and establishment of stably transfected cells SBI-piggyBac vector, GST-S100A9, were kindly provided by Professor T.C. He from the University of Chicago. To construct an S100A9 overexpression plasmid, the complete coding sequence of human S100A9 gene was subcloned into the SBI-piggyBac vector. For S100A9 silencing, siRNAs targeting human S100A9 with the sequences of 5'-GCAAGACGAUGACUUGCAA -3' and 5'-UUGCAAGUCAUCGUCUUGC -3'were synthesized and assembled into the SBI-piggyBac vector, resulting in SBI-siS100A9. After transfecting SH-SY-5Y NB cells with the constructed plasmids, the stably transfected cells were selected by incubation with puromycin for one week. The stably transfected cell lines, namely Control (SH-SY5Y transfected with SBI empty vector), S100A9(SH-SY5Y transfected with SBI-S100A9) and siS100A9(SH-SY5Y transfected with SBI-siS100A9) brie y.

Cell viability assay
The viability of SH-SY5Y cells was assessed by 3-(4, 5-dimethylthiazol-2-yl)-2, 5-diphe-nyltrazolium bromide (MTT) assay. Briefly, stably transfected SH-SY5Y cells were seeded in 96-well plates (1000 cells/well). The cells were incubated in DMEM supplemented with 1% FBS for 24, 48, 72, 96 and 120 h, then incubated with MTT reagent (Progema, Madison, WI, USA, 20 µL /well) for another 4 h at 37°C to allow the formation of formazan. After that, 100 µL of dimethyl sulfoxide was added into the cell culture medium for another 10-min incubation at room temperature. At last, in every day of the next ve days, a microplate reader (Bio-rad, iMark) was used to measure the absorbance at 492 nm of each well. The experiment included three independent replicates for each sample.

Colony formation assay
Exponentially growing stably transfected SH-SY5Y cells were seeded at a low density (100 cells/well) in cell culture medium containing 1% FBS in 6-well plates. The cells were allowed to grow for about 10 days to form colonies. The culture medium was refreshed every 3-4 days. Crystal violet was used to stained the colonies. Colony numbers from 3 wells were used to calculated the average colony number.

Scratch wound healing assay
The scratch wound healing assay was performed as described previously [20,21]. Brie y, stably transfected SH-SY5Y cells were seeded in 6-well plates and grown to ~90% con uency. Then, sterile micro-pipette tips were used to scratch the monolayer formed by SH-SY5Y cells to create a wound. After that, the medium (DMEM containing 1% FBS) was refreshed every day to remove the oating cells. Bright eld microscopy was used to monitor the wound healing status at 24 h, 48 h, and 72 h after the wound was created. Each assay was repeated three times. imageJ software was used to calculate the wound healing ratio.

Transwell invasion assay
A chamber coated with non-type I-collagen (Millipore, USA) was used for the Transwell assay. The upper chamber coated with ECM gel (Sigma, USA) was lled with 400 µL of serum-free DMEM and seeded with exponentially growing stably transfected SH-SY5Y cells (1×10 4 cells). The lower chamber was lled with 500 µL of DMEM supplemented with 20% FBS, which served as a chemoattractant. After 24 h of incubation, the cells migrated across the Transwell membrane were dried, xed with methanol, and then stained with hematoxylin-eosin (H&E). Cotton swabs were used to remove the cells on the upper surface of the Transwell membrane. At last, an inverted microscope (×100 magni cation) was used to count the number of cells migrated across the Transwell membrane. Five randomly-selected elds were examined to obtain the mean value of the number of cells migrated across the Transwell membrane. The experiment was repeated three times.

Statistical analysis
All data are presented as means ± standard deviations. T test was performed in the GraphPad Prism software to determine the statistical signi cance of differences between groups. A P value of less than 0.05 was considered statistically signi cant.

Enrichment analyses of DEGs
In general, 176 up-regulated DEGs and 19 down-regulated DEGs were identi ed based on the cut-off criteria(P-value < 0.05 and count ≥ 2), which were used to generate a heatmap for the NB patients with and without metastasis (Fig. 1A). Functional enrichment analysis revealed that the top 18 mostly enriched GO terms and top 2 mostly enriched KEGG pathways were associated with regulation of nucleosome assembly, innate immune response in mucosa, extracellular matrix organization, angiogenesis and innate immune response, etc. (Table 1). In addition, volcano plot was able to quick identify the expression changes within the gene sets by combination of statistical tests (adjusted p-value) and magnitude of change. (Fig. 1B) PPI network analysis To identify the potential biomarkers that can predict metastasis in NB patients, the DEGs were delineated to construct a PPI network. A DEG with a combined score (median con dence) > 0.15 was regarded as a signi cantly differentially expressed gene. As shown in Fig. 2A, the PPI network constructed by STRING is composed of 150 nodes and 452 interactions. The Cystoscope software was used to visualize the PPI network. The top 10 hub proteins in the PPI network determined by Closeness, Degree, EPC, MNC, Radiality, and Stress centrality are listed in Table 2. The results showed that TAC1, PTGS2 and FGF1 were the most outstanding genes and maybe play an important role in the NB metastasis. Sub-network analysis suggested that S100A9 proteins were the outstanding hub proteins (Fig. 2B).

Validation of the outstand gene in NB patients
To verify the expression of TAC1, PTGS2, FGF1 in microarray data, qPCR was performed to detect the expression of those three genes, nally, the qRT-PCR results showed that the expression of three gene there were no signi cant differences between the 9 NB patients with bone marrow metastasis and 10 NB patients without bone marrow metastasis (Fig. 3A). Next, S100A9 identi ed in the microarray data, ELISA and IHC were performed to detect the protein levels of S100A9 in NB patients (9 NB patients with bone marrow metastasis and 10 NB patients without bone marrow metastasis). The serum levels of S100A9 were higher from NB patients with bone marrow metastasis than that from NB patients without bone marrow metastasis. (Fig. 3B) Further, we also examined the expression of S100A9 in sections from NB patients with bone marrow metastasis tissues and NB patients without bone marrow metastasis using IHC staining (Fig. 3C). The results showed that the expression levels of S100A9 were signi cantly higher in NB patients with bone marrow metastasis than from NB patients without bone marrow metastasis. S100A9 overexpression promoted the proliferation, migration and invasion of SH-SY5Y cells To investigate the effects of S100A9 on the proliferation of SH-SY5Y cells, the coding sequence of human S100A9 gene expressed in GST-S100A9 vector was subcloned into the SBI-piggyBac plasmid to overexpress S100A9 in SH-SY5Y cells. According to the MTT assay results, the SH-SY5Y cells overexpressing S100A9 (de ned as the S100A9 group) exhibited higher proliferation ability than the SH-SY5Y cells transfected with empty vector (de ned as the control group) at days 3, 4, and 5 (P < 0.05). The S100A9-knockdown SH-SY5Y cells (de ned as the siS100A9 group) exhibited signi cantly slower proliferation when compared with the control group at days 3, 4, and 5 (all p values < 0.001) (Fig. 4A). Colony formation assay showed that the S100A9 group formed more colonies compared with the control group. Quantitatively, the number of colonies formed in the S100A9 group was approximately double than that of the control group (Fig. 4B). These results indicated that S100A9 overexpression accelerated the proliferation of SH-SY5Y cells.In addition, the results showed that the migration, and invasion of SH-SY5Y cells were signi cantly active by S100A9, as revealed by Transwell assays (Fig. 5A), wound healing (Fig. 5B).

Discussion
In the current study, the GSE90121 dataset which was deposited by Kaplan D. were downloaded and analyzed by bioinformatics method to identify potential crucial genes associated with NB metastasis. A total of 195 genes including 19 down-regulated and 176 up-regulated genes were obtained. Besides, the signi cantly enriched GO terms were mainly focused in regulation of cell growth, in ammatory response, positive regulation of cell migration. Hub genes of the regulatory network were then selected and conducted with PPI network module. Dysregulated TAC1, PTGS2 and FGF1 were the top 3 outstanding genes based on both six methods (Closeness, Degree, EPC, MNC, Radiality, and Stress centrality) evaluation. The expression levels of TAC1, PTGS2 and FGF1 in resected specimen of NB patients with or without metastasis were then validated by qRT-PCR, although the expression of TAC1, PTGS2 and FGF1 were related with many kinds of tumorigenesis, such as non-small cell lung cancer [22], pancreatic ductal cancer [23], colorectal cancer [24][25][26], squamous cell carcinoma [27], gastric cancer[28] and clear cell renal cell carcinoma [29] ,but these genes expression did not exhibit signi cant change as expected as the microarray results, indicating that these three genes may not the pivotal gene that participate in the metastasis of NB.
After go through the differentially expressed genes list and reviewing the relevant literatures, we found that S100A9 exhibits a broad range of biological functions involving in various cancer progression [30][31][32],speci cally, S100A9 proteins take part in a broad range of intracellular and extracellular functions by regulating calcium balance, cell apoptosis [33], migration [34], proliferation [35], differentiation[36, 37,9], and in ammation[38]But fewer research investigated the expression and biological function of S100A9 in NB cells. Here, the promoted expression levels of S100A9 were veri ed to be consistent with the microarray results. To further validate the results, we observed that elevated S100A9 promoted the proliferation and migration of NB cancer cells. These ndings suggested that S100A9 may be an important carcinogenic factor in the occurrence and progression of NB and may serve as an attractive therapeutic target for NB patients.

Conclusions
Collectively, the novel potential biomarker S100A9 was identi ed by bioinformatics methods in the present study. And our results suggested that S100A9 promote the proliferation and migration of NB cell and may play a pivotal role in NB metastasis and expected to a get a further insight into the molecular mechanisms of NB patients and would help us discriminate the metastatic possibility of NB patients.

Declarations
Con ict of interest All authors declare that they have no con ict of interest.
Ethical approval This study was approved by the Medical Ethics Committee of Tongji Medical College, Huazhong University of Science and Technology.
Consent to participate Informed consent were obtained from the participants or their parents. Consent for publication All authors had gone over the manuscript and consent to submit Funding This study is fully supported by the National Natural Science Foundation of China (No. 81601821 to YYZ), the Natural Science Foundation of Shandong Province (No. ZR2020MH314, to XC) Author Contributions YYZ conceived the study. XC, YKX, JF, QWT and QW perform the study. Statistical analysis was undertaken by XC. All authors read and approved the nal version of the manuscript. QW and YYZ wrote the manuscript.

Data Availability
The following information was supplied regarding data availability: The data is available at NCBI GEO: GSE90121 Figure 2 Protein-protein interaction network constructed with the up-and down-regulated differentially expressed genes.

Figure 3
PPI network was visualized by Cystoscope software. Figure 5 (A). MTT analysis for Blank/Control/S100A9/siS100A9 SH-SY-5Y cells for sequential 5 days. *P<0.05, **p < 0.01 (B). Colony formation assay for Blank/Control/S100A9/siS100A9 SH-SY-5Y cells. The representative images of transmembrane cells are shown in the right panel, the mean numbers of transmembrane cells ± SD per microscopic eld of three independent experiments are quanti ed in the right panel. **p < 0.01, ***p<0.001