Bioinformatics analysis of laryngeal squamous cell carcinoma: seeking key candidate genes and pathways

DEGs identification in LSCC

We obtained the gene expression profiles of LSCC from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) (Barrett et al., 2013). The raw gene expression profile datasets GSE59102 and GSE107591 were obtained from that database. The platform for GSE59102 is GPL6480, Agilent-014850 Human Genome Microarray 4 × 44 K G4112F (Probe Name version), including 13 normal tissues and 29 tumour tissues. The platform for GSE107591 is GPL6244, HuGene-1_0-st Affymetrix Human Gene 1.0 ST Array, containing 23 normal and 24 tumour tissues. The DEGs between cancerous and noncancerous tissues were detected by GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r) with the criteria of logFC (fold change) > 1 and adjusted p-value < 0.01, which were thought to indicate clearly differential expression (Davis & Meltzer, 2007). The next step was to use the online Venn software to identify the intersection of DEGs between the two datasets.

KEGG and GO enrichment analyses of DEGs

We used the Database for Annotation, Visualization and Integrated Discovery (DAVID, http://david.abcc.ncifcrf.gov/) online tool to study the roles of the identified DEGs and put into effect functional and pathway enrichment analysis (Huang et al., 2007). KEGG (http://www.kegg.jp), an integrated database, was applied to determine the high-level functions and the biological interpretation of genome sequences and other high-throughput data (Kanehisa et al., 2016). The GO (http://www.geneontology.org) project is a primary bioinformatics tool for annotating genes and analysing the biological processes (BPs) of the genes and includes the BP, cellular component (CC) and molecular function (MF) (Ashburner et al., 2000). p-value < 0.05 was thought to refer a important difference in statistic.

PPI network construction and module analysis

A PPI network of DEGs was foreseen by the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) online database (version 11.0; http://string-db.org) (Szklarczyk et al., 2015). Based on the STRING online tool, PPIs of the DEGs were constructed with a confidence score > 0.4, we exported the results as a simple tabular text. Cytoscape software version 3.7.1 (www.cytoscape.org) is an open source bioinformatics software platform for visualizing molecular interaction networks. Subsequently, the PPI network was visualized by means of Cytoscape software (Holmås et al., 2019). The most significant module in that network was detected by the Molecular Complex Detection (MCODE) (version 1.4.2) plug-in of Cytoscape (Bandettini et al., 2012). The selection standards: MCODE score > 5; degree cut-off = 2; node score cut-off = 0.2; max depth = 100; and k-score = 2. Afterward, KEGG and GO analyses of the genes in the PPI were performed by DAVID.

Hub gene selection and analysis

The CytoHubba plug-in of Cytoscape is an application that contains several topological algorithms for ranking nodes in a PPI network by the network characteristics. In this work, maximal clique centrality (MCC) was selected to explore the top nine hub genes among the 12 available computing methods (Chin et al., 2014). Moreover, the UALCAN (http://ualcan.path.uab.edu/) and GEPIA (http://gepia.cancer-pku.cn/) online databases were resorted to validate gene expression. UALCAN and GEPIA are internet tools that provide a straightforward way to investigate publicly available cancer transcriptome data, including The Cancer Genome Atlas (Chandrashekar et al., 2017; Tang et al., 2017). Moreover, the GEPIA database was applied to survival analysis of the hub genes. p-Value < 0.05 was considered as inportant in statistic.

Cell culture and cell transfection

The human LSCC cell lines Hep2 and LSC-1 (Bluefbio, China) and human laryngeal epithelial cells (HLECs; Lifeline, Framingham, MA, USA) were cultured in Dulbecco’s modified Eagle’s medium (DMEM) containing 10% foetal bovine serum, 100 μ/mL streptomycin and 100 μ/mL penicillin, the culture conditions were 37 °C and 5% CO₂. All cells were passaged 3–4 times before use (Wang et al., 2018).

The siRNA-SPP1, siRNA-SERPINE1 and siRNA-MMP1 that were transfected into Hep2 and LSC-1 cell lines by Lipofectamine 2000 were obtained from GeneCreate (China) in accordance with the instuctions of manufacturer. Hep2 and LSC-1 cells were cultured in DMEM at least 24 h and before transfection, they also need to washed with phosphate-buffered saline (PBS, pH 7.4) before transfection (Liu, Ren & Song, 2019).

Western blot

Total protein isolated from Hep2 and LSC-1 cells and HLECs were loaded into sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE) gels and transferred to polyvinylidene fluoride (PVDF) membranes. The PVDF was hindered with 5% skim milk at 37 °C for 90 min. Next, the PVDF membrane was incubated with SPP1, SERPINE1 or MMP1 rabbit polyclonal antibodies (Proteintech, Rosemont, IL, USA) and β-actin rabbit polyclonal antibodies (ABclonal, Woburn, MA, USA) at 4 °C overnight. And the PVDF was rinsed with PBST for 10 min and repeated three times. Then, the PVDF was incubated with horseradish peroxidase-labelled goat anti-rabbit IgG secondary antibody (Jackson, MS, USA) at 37 °C for 60 min. The PVDF was rinsed with PBST for 10 min and repeated three times. Then, the immunoactivity was detected by optical luminometry (Mishra, Tiwari & Gomes, 2017).

qRT-PCR

To isolate total RNA from Hep2 cells, LSC-1 cells and HLECs, we used TRIzol reagent (TAKARA, Japan) according to the instructions. The conditions of reverse transcription reaction can be noticed: 25 °C for 10 min, 50 °C for 30 min and 85 °C for 5 min. Detection using the fluorescence quantitative PCR kit. The conditions of fluorescence quantitative PCR were as follows: 95 °C for 5 min, 95 °C for 10 s, 60 °C for 30 s, a total of 40 cycles. The solubility curve temperature range was set at 60–95 °C, and three replicate wells were prepared for each specimen. SPP1, SERPINE1 and MMP1 used β-actin as internal factors (Wilhelm & Pingoud, 2003). The qRT-PCR results were determined by the 2^−ΔΔCt approach, and the primer sequences are listed in Table 1.

CCK-8 assay

To determine the proliferation of Hep2 and LSC-1 cells, we used a CCK-8 assay. We incubated the cells in 96-well culture plates and the inoculation density was 3 × 10⁵. For cell transfection, the cells were cultured overnight. After transfection for 48 h, 10 mL of CCK-8 solution was added to each well, and the cells were incubated at 37 °C for another 60 min. A total of 490 nm was used as the absorbance of the solutions, and the absorbance of the solutions was detected by a Smart Microplate Reader 16.1 (Ma et al., 2017).

Transwell migration assay

For cell migration tests of Hep2 and LSC-1 cells, Transwell chambers with a polycarbonate membrane were used. In the upper chambers, Hep2 and LSC-1 cells were incubated in serum-free DMEM, and 10% FBS was out to the lower chambers. After 10 h, the Hep2 and LSC-1 cells in the upper chambers were removed, and the cells in the lower chambers were dyed with crystal violet at 25 °C for 1 min. A light microscope was used to observe and count the cells in five randomly selected fields.

Wound healing assay

Hep2 and LSC-1 cells were inoculated into the 6-well plates. When cells were at 90–100% confluence, wounds of uniform width were created by slowly pulling a 10 μL pipette perpendicular to the bottom of the 6-well plate (three wounds/well). Then, the cells were rinsed with PBS three times to remove floating cells, and incubated in serum-free DMEM to inhibit cell proliferation and division. At 0 and 24 h after wounding and culture, the migration distance of cells in the wound area was observed under a microscope, and several different fields of view were randomly selected for photographing (He et al., 2020). The experiment was repeated three times.

Statistical analysis

Statistical Product and Service Solutions 25.0 software was ran for analysis in statistic. The t-test was applied to numerical data, and p < 0.05 indicated the significance level.

DEGs identification in LSCC

Before identifying DEGs in LSCC, the microarray was standardized. We identified 3,359 DEGs in GSE59102 and 444 DEGs in GSE107591. The overlap between the two datasets included 341 genes (Fig. 1), including 205 downregulated genes and 136 upregulated ones between LSCC tissues and control tissues.

KEGG and GO enrichment analyses of DEGs

Gene Ontology and KEGG pathway enrichment analyses were performed in DAVID (Table 2) to study the biological classification of the DEGs. The BP category outcomes of the GO analysis outcomes expressed that the DEGs were markedly associated with extracellular matrix (ECM) organization, cell adhesion, the collagen catabolism pathway and positive regulation of cell-substrate adhesion. Within the CC categories, the DEGs were markedly linked with the extracellular region, ECM and proteinaceous ECM. Among the MF categories, the DEGs were significantly enriched in heparin binding and ECM structural constituents. In addition, KEGG signalling pathway analysis illustrated that the DEGs were pivotal in the following pathways: ‘ECM-receptor interaction,’ ‘drug metabolism-cytochrome P450,’ ‘focal adhesion,’ ‘PI3K-Akt signaling pathway’ and ‘complement and coagulation cascades.’

PPI network construction and remarkable module identification

The STRING database was used to predict the interactions among the DEGs with a combined score >0.4 at the protein level. The PPI network was built with Cytoscape software and included 176 nodes and 371 edges (Fig. 2A). Additionally, the hubs of the PPI network module were obtained by MCODE and consisted of 9 nodes and 34 edges (Fig. 2B). The functional analyses of genes contained with the modules were identified by running DAVID. The outcomes indicated that the genes of this module were highly abundant in condensed nuclear chromosome outer kinetochore, ATP binding and ATP-dependent microtubule motor activity, plus-end-directed (Table 3). To more deeply investigate the top nine hub genes, we used the cytoHubba plug-in of Cytoscape in the above PPI using the MCC method. The hub genes MMP9, SPP1, SERPINE1, MMP1, MMP13, MMP3, CXCL8, OSM and COL1A1 were used for further analysis. The names, abbreviations and synonyms for these hub genes can be seen in Table 4.

High expression of SPP1, SERPINE1 and MMP1 in LSCC and their association with poor prognosis

Tools in the UALCAN and GEPIA databases were used to evaluate the expression levels of these nine genes in LSCC. Both databases confirmed that the expression of MMP9, SPP1, SERPINE1, MMP1, MMP13, MMP3, CXCL8, OSM and COL1A1 showed obvious differences in LSCC samples and healthy samples, which was in accordance with the results discussed above (Fig. 3). Using the GEPIA database to study the relationship between gene expression and patient survival, we found that only SPP1, SERPINE1 and MMP1 were obviously connected with the overall survival of patients with LSCC (Fig. 4). Therefore, we verified the expression of SPP1, SERPINE1 and MMP1 in LSCC cells by western blot (Figs. 5A–5D) and qRT-PCR (Fig. 5E). The outcomes indicated that SPP1, SERPINE1 and MMP1 were highly expressed in HEP2 and LSC-1 cells. These results were consistent with the bioinformatics analysis.

Downregulation of SPP1, SERPINE1 and MMP1 expression inhibited LSCC cell proliferation in vitro

To further verify the function of SPP1, SERPINE1 and MMP1 in HEp-2 and LSC-1 cells, SPP1, SERPINE1 and MMP1 were knocked down in vitro. CCK-8 was performed to verify the effects of SPP1, SERPINE1 and MMP1 on HEp-2 and LSC-1 cell proliferation in vitro. HEp2 and LSC-1 cells treated with siRNA-SSP1, siRNA-SERPINE1, siRNA-MMP1 and siRNA-NC were cultured and then subjected to CCK-8 assays. The results showed that the growth rate of HEp-2 and LSC-1 cells was significantly inhibited after knockdown of SPP1, SERPINE1 and MMP1 (Fig. 6).

Downregulation of SPP1, SERPINE1 and MMP1 expression inhibited LSCC cell migration and invasion in vitro

Invasion and metastasis are important features of tumours. Therefore, the effects of SPP1, SERPINE1 and MMP1 on migration and invasion of HEp-2 and LSC-1 cells were analysed in depth. First, transwell migration assay was performed, and the results were consistent with those of the wound healing assay. After knockdown of any of these genes, the amount of cells passing through the filter membrane to the lower chamber was decreased significantly (p < 0.05, Figs. 7A–7J). Next, a wound healing assay was used to detect the cell migration capacities in the control group and the SPP1, SERPINE1 and MMP1 knockdown groups. The results showed that the migration distance of cells in the SPP1, SERPINE1 and MMP1 knockdown groups was shorter when comparing with that in control set after 24 h (p < 0.05, Figs. 7K and 7L).