BIRC5 promotes cancer progression and predicts prognosis in laryngeal squamous cell carcinoma

Background Laryngeal squamous cell carcinoma (LSCC) remains one of the most common respiratory tumors worldwide. Baculoviral inhibitor of apoptosis repeat containing 5 (BIRC5) is a member of the inhibitor-of-apoptosis protein family. BIRC5 plays an important role in various types of cell proliferation, differentiation, migration and invasion. However, the specific role of BIRC5 in LSCC remains unclear. Methods To provide a prognostic biomarker for LSCC, we screened the prognostic genes of LSCC via bioinformatics. PPI network and KEGG pathways were used to select hub genes. Clinical prognoses were performed using a Kaplan–Meier plotter and Cox proportional-hazard analysis. BIRC5 expression in LSCC tissues and cell lines were detected by RT-PCR, Western blot and Immunohistochemistry (IHC). Cell proliferation, cell cycle and cell apoptosis were detected with Cell Counting Kit-8 (CCK-8) and Flow Cytometry assay, respectively. Results Here, BIRC5 was strongly correlated with higher tumor grade and differentiation. BIRC5 was highly expressed in LSCC tissues when compared with normal tissues and increased expression of BIRC5 was associated with overall survival in LSCC patients. The suppression of BIRC5 induced cell apoptosis and cell cycle arrest, thereby inhibiting the proliferation of LSCC cells. The survival analysis confirmed that higher level of BIRC5 expression predicted poor prognosis of LSCC patients. BIRC5 may act as an oncogene of LSCC development and was suggested as a promising prognostic biomarker for LSCC.


INTRODUCTION
Laryngeal squamous cell carcinoma (LSCC) is the second most common form of head and neck carcinoma (HNC) with high incidence rate (177,422 new cases) in 2018 (Bray et al., 2018). LSCC includes malignant squamous lesions in the larynx, maxillofacial and pharynx, of which more than 90% are squamous cell carcinoma. Currently, the occurrence of laryngeal carcinoma in China is linearly increasing and approximately four times higher than in the United States (Torre et al., 2015). Despite considerable advances

Data acquisition
The GEPIA (http://gepia.cancer-pku.cn/index.html) was used to confirm the relative levels of gene expression in normal and tumor samples. The selection criteria are as follows: Log 2 FC cutoff: 1, P-value cutoff: 0.01. The Ualcan (http://ualcan.path.uab.edu/) was used to validate the relative level of BIRC5 in tumor subgroups according to different grades (Chandrashekar et al., 2017). Genes were pasted into the text and chose the cancer type. The relative expression level and survival information of the gene in normal vs. tumor samples and across cancer subgroups were output.

Kaplan-Meier plotter database analysis
The Kaplan-Meier Plotter Database (http://kmplot.com/analysis/) was used to analyze the correlation between BIRC5 expression and survival in head and neck squamous cancer (HNSC). We analyzed the overall survival of HNSC patients by using a Kaplan-Meier survival plot. First, chose Pan-cancer RNA-seq, then genes were uploaded into the database, and the inclusion criteria was: Head-neck squamous cell carcinoma (n = 500), Follow up threshold: 240 months, and number-at-risk were also shown. The hazard ratio (HR) with 95% confidence intervals and long rank P-value were calculated and displayed on the webpage.

Microarrays in GEO datasets and profiles
The Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo) is an open database of microarray/gene profiles (Barrett et al., 2011). In order to understand the diagnostic value of BIRC5 in LSCC, two mRNA expression datasets containing microarray data from normal tissues and LSCC tissues, GSE51985 (Lian et al., 2013) and GSE59102 (Placa et al., 2015), were selected. The GSE51985 microarray data including 10 adjacent non-neoplastic tissues and 10 LSCC tissues and were generated using Illumina HumanHT-12 V4.0 expression beadchip. The GSE59102 microarray data including 13 normal tissues and 29 LSCC tissues and were generated using Agilent-014850 Whole Human Genome Microarray 4 × 44K G4112F.

Cell cycle analysis
Cells were grown in six-well plates until they achieved 30-40% confluence. After transfected with shBIRC5 or shNC for 48 h, the cells were fixed with 75% ethanol at 4 C overnight, then centrifuged at 300×g for 5 min. RNase A (50 µg/ml) and propidium iodide (50 µg/ml) (Sigma, St. Louis, MO, USA) were added into cells at room temperature for 15 min. Cell cycle was analyzed by BD Accuri TM C6 flow cytometry (BD Biosciences, Franklin Lakes, NJ, USA) according to the manufacturer's protocol and determined by ModFit LT 4.1 software (BD Biosciences, Franklin Lakes, NJ, USA), CXP software was used to determine the cell percentage in different phases. The linearity was set to two and internal standards were also selected.

Cell proliferation and colony-formation assays
Cells were seeded in 96-well plates at 2 × 10 3 per well for 24 h. Cell proliferation was evaluated using CCK8 Assay kit according to the manufacturer's instructions. CCK8 reagent was added to cell culture medium after 24, 48, 72, and 96 h, and incubated for 2 h. The absorbance at 490 nm was measured using a Spark Ò multimode microplate reader (Tecan, Männedorf, Switzerland). Each experiment was repeated three times. For colony-formation assays, 500 cells were seeded in six-well plate and incubated in complete growth medium for 7-14 days. Cell colonies were fixed with methanol and stained with 0.1% crystal violet (Sigma-Aldrich, St. Louis, MO, USA) in order to be counted manually.

Bioinformatics analysis
To screen the DEGs between LSCC tissues and normal tissues, two LSCC gene expression series were analyzed by GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/). The cutoff value criteria were: P value < 0.05 and |log 2 FC| > 2. For the overlapping genes, the DAVID database (http://david.ncifcrf.gov/) was used to perform GO enrichment and KEGG pathway analysis (Huang, Sherman & Lempicki, 2009). Venn diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/) was used to calculate the intersections of the list of elements. Further, we used mcode to analyse the GO enrichment pathway including biological process, cellular component and molecular function (Shannon et al., 2003). Meanwhile, the STRING database (http://string-db.org/) was designed to provide a protein-protein interaction network analysis (Szklarczyk et al., 2015). Hub genes in LSCC were obtained from the network. Univariable and multivariable Cox regression were performed to analyses survival-associated genes.

Western blot
The cells were transfected with shBIRC5 and shNC for 48 h and harvested with lysis buffer (Solarbio, Beijing, China) including protease inhibitors (Roche, Basel, Switzerland). Overall, 30 µg total protein was separated by SDS-PAGE and transferred to PVDF membranes (Millipore, Burlington, MA, USA) and blocked with 5% bovine serum albumin, then the membranes were incubated with primary antibodies overnight at 4 C. Immunoblots were detected using ECL detection reagent (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer's protocol. Antibodies used were: BIRC5 (sc-17779; Santa Cruz), GAPDH (#5174; Cell Signaling Technology).

Lentivirus transduction
The lentiviral shRNA plasmid targeted BIRC5 were from Ribo Bio Company (Guangzhou, China). The BIRC5 shRNA sequence was: CCGGCCTTTCTGTCAAGAAGCAGTTC TCGAGAACTGCTTCTTGACAGAAAGGTTTTTG. Briefly, the recombinant lentivirus was harvested by co-transfecting 293T cells. After an additional 72 h, culture media were collected and the aliquot supernatant was stocked at −80 C for infection in future. The efficiency of BIRC5 knockdown were evaluated by qPCR and western blot.

Immunohistochemistry (IHC) and antibodies
The Human Protein Atlas (HPA) is an online antibody-based proteomics database of human disease and normal tissues (https://www.proteinatlas.org/). In this study, representative BIRC5 expression was obtained from HPA in normal tissues and LSCC tissues to demonstrate the potential expression pattern in clinical practice. Our cohort tissues from 88 cases of LSCC fixed in the formalin and embedded in the paraffin for BIRC5 immunohistochemistry staining. BIRC5 antibody for IHC was obtained from Santa Cruz (1:1,000, sc-17779).

Statistical analysis
Stata 12.0 was used to perform a GEO microarray meta-analysis in this study. Statistical analyses were performed using GraphPad Prism 7.0a and SPSS 20.0 software. Two-tailed Student's t-test or ANOVA were used to compare experimental groups as indicated. Cox proportional hazards regression model were used for multivariable analyses. The data were presented as mean ± SD. P < 0.05 was considered to be statistically significant.

Bioinformatics analysis of the DEGs in LSCC
A total of 1,287 LSCC DEGs were screened from two gene expression profiles (GSE51985 and GSE59102) in the GEO database. Then, the upregulated and downregulated DEGs of GSE51985 and GSE59102 showing an overlap region (Figs. 1A and 1B). Moreover, compared with normal tissues, 41 DEGs with absolute P < 0.05 and log 2 FC > 2 were deemed to be upregulated and 84 DEGs with P < 0.05, log 2 FC < -2 were downregulated in LSCC tissues (Tables S1 and S2). The DEGs enrichment of GO terms (biological process, molecular function and cell component) and KEGG pathways analysis were performed by DAVID. Fig. 1C indicates the top 10 highly enriched GO items and the significantly enriched KEGG pathways were shown in Fig. 1D. The DEGs were enriched in the Human papillomavirus infection, cellular senescence, cell cycle and pathways in cancer.

PPI network analysis and KEGG reanalysis
The PPI analysis among DEGs was constructed by STRING database. The network displayed 79 upregulated genes and 78 downregulated genes ( Fig. 2A). These genes were further visualized by Cytoscape MCODE, of which MELK, MCM2, TPX2, FOXM1, CEP55, CHEK1, CDCA5, TK1 and BIRC5 genes were all upregulated ( Fig. 2B and Table 1). DAVID was applied to reanalysis the biological functions of the above nine DEGs. Additionally, 3 out of 9 genes associated with the apoptosis pathway (Fig. 3).

Core genes expression and prognosis in LSCC
Based on the TCGA and GTEx, the level of core genes expression in normal and LSCC patients were validated by GEPIA. Fig. 4 indicated the expression of nine core genes in LSCC patients and normal patients. Then, we investigated the prognostic value of core genes using Kaplan-Meier analysis. As shown in Fig. 5A, among the above genes only high BIRC5 showed shorter overall survival when compared with other genes. Further, we evaluated BIRC5 expression in LSCC by using the RNA-seq data in TCGA (Fig. 5B). By mining previously reported data, we found that high level of BIRC5 was significantly associated with higher tumor grade (Fig. 5C).

BIRC5 is overexpressed in LSCC patients
To further determine the expression level of BIRC5 in LSCC patients, we analyzed 88 primary LSCC patients and matched adjacent normal tissues by immunochemistry in our cohort. The histological types of LSCC were determined according to the World Health Organization system. BIRC5 was significantly upregulated in 88 LSCC tissues (Figs. 6A and 6B). Kaplan-Meier survival analysis of our cohort showed that higher BIRC5 patients had shorter overall survival, which was consistent with the online database results (Fig. 6C). We then validated the protein level of BIRC5 from HPA and found BIRC5 was also increased in LSCC tissues (Fig. 6D). In addition, RT-PCR analysis showed that the expression level of BIRC5 in LSCC was higher than that of matched adjacent normal tissues (Fig. 6E). These findings suggested a positive correlation between BIRC5 upregulation and LSCC progression.

BIRC5 serves as an independent prognostic marker in LSCC patients
To explore the relationship between BIRC5 and LSCC clinicopathologic parameters, we analyzed the corresponding clinical data in 88 pairs of LSCC tissues. The clinical characteristics were showed in Table 2. The results indicated that BIRC5 expression was related with advanced tumor grade and differentiation. Then, we performed univariable and multivariable cox regression analysis to select independent clinicopathological prognostic factors for LSCC from TCGA cohort. Multivariable analysis indicated that BIRC5 expression and differentiation grade were independent factors for OS (Table 3). Taken together, BIRC5 can serve as an independent prognostic factor in patients with LSCC.

BIRC5 promotes cell proliferation and controls cells cycle and apoptosis in vitro
To further verify the above bioinformatics analysis, BIRC5-specific shRNA was exploited to detect with RT-PCR and western blot analysis. BIRC5 was effectively knocked down in  Tu212 and Tu686 cells (Figs. 7A and 7B). To investigate the role of BIRC5 on LSCC cell proliferation, CCK8 assay and colony formation assay were performed. CCK8 assay revealed that the proliferation ability was reduced after BIRC5 depletion (Figs. 7C and 7D). Silencing BIRC5 significantly inhibited the colony-forming capabilities of Tu212 and Tu686 cells (Fig. 7E). Next, we further explored the functional impact of BIRC5 expression on cell growth behavior of LSCC cells. Flow cytometry was performed to analyze the changes in cell cycle, as shown in Fig. 7F, cells with BIRC5 suppression revealed that cell cycle was arrest in G1 phase, confirming the role of BIRC5 in controlling cell cycle progression. In addition, the percentage of apoptotic cells in the shBIRC5-transfected group was shown higher than the shNC group (Fig. 7G).

DISCUSSION
In recent decades, laryngeal squamous cell carcinoma (LSCC) is the most common malignancies in laryngeal cancer. Despite comprehensive therapeutic advances for LSCC have improved, the outcomes for patients with advanced LSCC have not improved in the last 30 years (Jenckel & Knecht, 2013). LSCC is prone to local and distant invasion, and lymph node metastasis, which are the main factors for poor prognosis (Gao et al., 2019). Therefore, discovering more effective biomarkers that can accurately signify the prognosis of patients would be essential for LSCC diagnosis and treatment. The present study screened two gene expression profiles from GEO database. Using differential gene expression analysis, A total of 125 overlapping DEGs were found. These DEGs were enriched in the Human papillomavirus infection, cellular senescence, cell cycle and pathways in cancer. Also, the Cytoscape MCODE was used to screen network modules and identify hub genes. The hub genes were MELK, MCM2, TPX2, FOXM1, CEP55, CHEK1, CDCA5, TK1 and BIRC5. The survival analysis showed that only the high expression of BIRC5 was associated with poor overall survival. BIRC5, also known as survivin, was a member of the inhibitor-of-apoptosis proteins family identified and characterized in 1997 (Ambrosini, Adida & Altieri, 1997). BIRC5 is located in the nucleus and regulated cell cycle-related factors. It is essential for mitotic checkpoint and participates in cell differentiation, proliferation and invasion (Wang et al., 2012;Yamamoto, Ngan & Monden, 2008;Zwerts et al., 2007). However, the significance of BIRC5 expression in the prognosis of LSCC remains unclear. This is the first study to evaluate BIRC5 as a potential predictive biomarker for LSCC. We performed integrated bioinformatics analysis to explore the expression and prognostic value of BIRC5 in LSCC. Here, we found that the expression of BIRC5 was elevated in patients with LSCC by analyzing TCGA database. Kaplan-Meier plotter and multivariate cox regression analysis validated that the increased expression of BIRC5 was associated with overall survival in LSCC patients. With the data from UALCAN dataset showed that BIRC5 had higher expression levels in advanced LSCC tumor grade.
The results of western blot showed that BIRC5 was notably higher in LSCC cells compared with the HOK cell (data not show), which was consistent with mRNA level. These results collectively elucidated that the expression of BIRC5 could be a predictive biomarker for the prognosis of LSCC. Dysregulated BIRC5 have been linked to therapeutic resistance in stage II/III breast cancer (Hamy et al., 2016). Moreover, upregulated BIRC5 was associated with poor prognosis of ovarian cancer (Zhao et al., 2017). Dong et al. previously demonstrated that the expression of survivin is associated with unfavorable clinical pathological parameters and could serve as a biomarker for prognosis of LSCC patients (Dong et al., 2002). Another study by Pizem, Cor & Gale (2004) a high level of survivin expression predicts poor survival of LSCC. Consistently, we found that BIRC5 was tightly correlated with LSCC progression and tumorigenesis. Of note, our original prediction was based on TCGA and GEO datasets analyses, we found BIRC5 was strongly correlated with higher tumor grade of LSCC. More importantly, our Immunohistochemical (IHC) results validated that BIRC5 was highly expressed in LSCC tissues when compared with normal tissues. RT-PCR, western blot, colony-formation assay, and flowcytometry were designed to reveal the relationship between BIRC5 expression and LSCC progression. Collectively, as observed in some other cancer types (Cao et al., 2019;Liu, Lin & He, 2019), our findings highlight the potential of BIRC5 as a novel biomarker for LSCC patients.
Although our current study have screened and validated the potential prognostic biomarker of LSCC, some limitations still exist. Such as the number of samples from the database is insufficient to avoid bias. The clinical characteristics of human tissues used in this study are incomplete. Our results might provide valuable information and pathways concerned with LSCC for better clarifying the molecular mechanism of LSCC carcinogenesis.

CONCLUSIONS
In summary, the present study identified nine hub genes by data mining from TCGA and GEO datasets, and further analysis demonstrated that BIRC5 might be a prognostic biomarker in LSCC. Our findings revealed that BIRC5 could promote laryngeal cancer tumorigenesis and contribute to the poor prognosis of LSCC patients. BIRC5 may serve as an independent novel prognostic factor for LSCC.