Overexpression of LIMA1 Indicates Poor Prognosis and Promotes Epithelial-Mesenchymal Transition in Head and Neck Squamous Cell Carcinoma.

Background: LIMA1 encodes LIM domain and actin binding 1, a cytoskeleton-associated protein whose loss has been linked to migration and invasion behavior of cancer cells. However, the roles of LIMA1 underlying the malignant behavior of tumors in head and neck squamous cell carcinoma (HNSC) are not fully understood.
Methods: We conducted a multi-omics study on the role of LIMA1 in HNSC based on The Cancer Genome Atlas data. Subsequent in vitro experiments were performed to validate the results of bioinformatic analysis. We first identified the correlation between LIMA1 and tumor cell functional states according to single-cell sequencing data in HNSC. The potential downstream effects of LIMA1 were explored for gene ontology and Kyoto Encyclopedia of Genes and Genomes pathways through functional enrichment analysis of the gene sets that correlated with LIMA1 in HNSC. The prognostic role of LIMA1 was assessed using the log rank test to compare difference in survival between LIMA1High and LIMA1Low patients. Univariate Cox regression and multivariate Cox regression were further carried out to identify the prognostic value of LIMA1 in HNSC.
Results: LIMA1 was identified as a prognostic biomarker and is associated with epithelial-mesenchymal transition (EMT) progress in HNSC. In vitro silencing of LIMA1 suppressed EMT and related pathways in HNSC. Conclusions: LIMA1 promotes EMT and further leads to tumor invasion and metastasis. Increased expression of LIMA1 indicates poor survival, identifying it as a prognostic biomarker in HNSC.

Introduction

Head and neck squamous cell carcinoma (HNSC) is one of the most common malignancies worldwide, which exhibits the characteristics of diagnosis at an advanced stage and high recurrence, appearing as a heterogeneous group of tumors derived from the squamous epithelium of the upper aerodigestive tract, including the oral cavity, oropharynx, nasopharynx, larynx, and hypopharynx. According to updated cancer statistics, there were approximately 377 713 new cases of oral cavity cancer, 184 615 cases of laryngeal cancers, 133 354 cases of nasopharyngeal cancers, 98 412 cases of oropharyngeal carcinoma, and 84 254 cases of hypopharyngeal cancers worldwide in 2020.

The high incidence of recurrence or metastatic disease in HNSC means that the prognosis of patients with HNSC stage III/IV is not ideal. With a 5-year survival rate limited to 40%, the management of patients with metastatic HNSC has received increased interest from the HNSC oncology community. Although huge strides have been made in the diagnosis and treatment of HNSC, therapy for HNSC is still unsatisfactory. The pattern of clinical behavior and response to treatment of recurrent and/or metastatic HNSC is heterogeneous, necessitating the identification of novel prognostic biomarkers to increase its diagnostic efficacy, and to avoid unnecessary toxicities, eventually leading to prolonged survival. Studies have screened for prognostic markers or therapeutic targets[7,8]; however, few have focused on genomic biomarkers related to epithelial-mesenchymal transition (EMT). Further study is needed to gain an insight into genomic profiles of HNSC targeting EMT. The EMT is a biological process (BP) implicated in tumor metastasis, in which epithelial cells lose their epithelial features and gain mesenchymal features. The EMT-derived tumor cells exhibit stem cell properties characterized by remarkable therapeutic resistance; therefore, targeting EMT pathways constitutes an attractive strategy for cancer treatment. However, the underlying molecular mechanisms of EMT are not fully understood.

LIMA1 (LIM domain and actin binding 1), also known as EPLIN (epithelial protein lost in neoplasm), is preferentially expressed in human epithelial cells. As an actin-binding protein, LIMA1 was defined as a suppressive factor that is frequently downregulated in epithelial tumors. LIMA1 plays a vital role in the adjustment of actin dynamics and has multiple correlations with epithelial cell junctions. Thus, downregulation of LIMA1 might significantly affect the migration and invasion of cancer cells, thereby increasing their metastatic potential. Accumulating studies have demonstrated that LIMA1 depletion promotes EMT and correlates with clinical lymph node metastasis in prostate cancer cells and melanoma cells. Recently, p53 was demonstrated to directly regulate LIMA1 at the transcription level, in which TP53 mutation caused downregulation of LIMA1 and poor survival of cancer patients. However, the role and potential prognostic significance of LIMA1 in HNSC is unclear.

In this study, LIMA1 was comprehensively studied via its expression, methylation, copy number alteration, and mutation profile followed by functional analysis of gene sets that correlated with its expression and genomic alterations. In addition, the prognostic significance of LIMA1 was clarified to determine its role in the malignant behavior of HNSC. Further cell line experiments were performed to verify the multi-omics analysis.

Materials and Methods

Expression and methylation data mining

The mRNA expression profile of LIMA1 in HNSC was analyzed online through GEPIA (http://gepia.cancer-pku.cn/). The relative protein expression level was assessed in The Human Protein Atlas, available from http://www.proteinatlas.org. The immunohistochemistry (IHC) images of LIMA1 staining are available at https://www.proteinatlas.org/ENSG00000050405-LIMA1/pathology/head+and+neck+cancer#img. The methylation data were extracted and analyzed using MEXPRESS (https://mexpress.be/).[18,19] Correlation analysis was conducted to identify the correlation between methylation and expression. The prognostic value of LIMA1 methylation was assessed by SurvivalMeth (http://bio-bigdata.hrbmu.edu.cn/survivalmeth/). UALCAN (http://ualcan.path.uab.edu/) was used to further investigate the relationship between LIMA1 mRNA expression or methylation with the corresponding clinical data.

Copy number variation and somatic mutation analysis

We explored and analyzed the copy number variation (CNV) data through cBioportal (http://www.cbioportal.org/), which allows the exploration of CNV events (amplification or deletion) of top gene sets that correlated with LIMA1 genomic alterations. The top 10 gene sets were selected for further protein-protein interaction (PPI) network functional analysis by GeneMANIA (http://genemania.org/). The HNSC genetic mutation data were downloaded from The Cancer Genome Atlas (TCGA) database. To reveal the somatic mutations profile of LIMA1, we used the “maftools” package in the R software, a visualized horizontal histogram displaying the genes that possess a higher mutation frequency in patients with high or low LIMA1 expression.[24,25]

Prognosis analysis

The RNA-sequencing data (Raw counts, level 3) and its corresponding clinical information in HNSC were obtained from the TCGA dataset (https://portal.gdc.cancer.gov/). A Kaplan-Meier survival analysis with the log-rank test was used to distinguish the survival difference between specific groups. A time-dependent receiver operating characteristic (ROC) analysis was used to compare the predictive accuracy of LIMA1 and the risk score.[26,27] Univariate and multivariate Cox regression analyses aimed to identify LIMA1-associated terms and thus to build the nomogram. The forest plot, along with the statistical P value, hazard ratio (HR), and 95% confidence interval (CI) of each variable was achieved through the “forestplot” R package. Based on the results of multivariate Cox proportional hazards analysis, a nomogram was developed to predict the total recurrence rate over X years. The nomogram provides a graphical representation of the factors that can be used to calculate the recurrence risk of a single patient by calculating the relevant points for each risk factor through the “RMS” R software package.[28 -30]

Gene function analysis

LinkedOmics (http://www.linkedomics.org/login.php) was used to analyze multi-omics data from the TCGA. One of the analysis modules, LinkFinder, compared genomic mRNA expression data to obtain the genes that correlated with LIMA1 in the TCGA HNSC cohort (n = 522). Comparisons were carried out statistically using the Spearman coefficient test. The LinkInterpreter module of LinkedOmics implemented pathway enrichment analysis. The results from LinkFinder module analysis were further subjected to gene set enrichment analysis (GSEA) of tumor-associated gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. The enriched results were based on the Molecular Signatures Database (MSigDB). The ranking criterion was a false discovery rate (FDR) < 0.05, and 500 simulations were performed. The TIMER 2.0 web server (http://timer.cistrome.org/)[32 -34] was used to evaluate the expression of LIMA1 and the corresponding immune cell infiltration levels. We investigated the functional states of HNSC cells with LIMA1 expression using CancerSEA (http://biocc.hrbmu.edu.cn/CancerSEA/), which included 14 cancer single-cell functional states from heterogeneous cancer cells.

Cell culture

The human oral squamous cell carcinoma (OSCC) cell lines CAL27 and HSC4 were obtained from the American Type Culture Collection (ATCC, Manassas, VA, USA). We cultured the cells in Dulbecco’s Modified Eagle Medium (DMEM) with 10% fetal bovine serum (FBS), 100 mg/mL streptomycin, and 100 U/mL penicillin (Beyotime, Jiangsu, China), at 37°C in 5% CO2.

Quantitative real-time reverse transcription–polymerase chain reaction

RNA was extracted from cells using TRIZOL (Takara, Shiga, Japan) and reverse transcribed using SYBRVR Premix Ex TaqTM Reverse Transcription-PCR (polymerase chain reaction) kit (Takara), according to the manufacturer’s instructions. Then, quantitative real-time PCR was accomplished using the ABI PRISMVR 7500 Fast Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). The relative mRNA levels were analyzed using the 2-△△CT method with GAPDH as the internal control.

Small interfering RNA transfection

We accomplished LIMA1 knockdown by constructing lentiviral vectors that contained the target gene and corresponding small interfering RNA (siRNA) sequences, constructed by GeneCopoeia Co., Ltd. (Guangzhou, China), following the manufacturer’s instructions. The efficiency of transfection was measured using quantitative real-time reverse transcription–polymerase chain reaction (qRT-PCR).

Wound healing assay

Transfected CAL27 and HSC4 cells were placed into a 12-well plate and grown to confluence. A sterile tip was used to scratch a linear wound, and the cells were incubated with serum-free medium for 24 hours. The wounds were digitally photographed at 0 and 24 hours, and the extent of wound healing was determined to assess the cell migration rate.

Migration assay

We performed a migration assay to assess the migration capacity of transfected CAL27 or HSC4 cells using 24-well Transwell chambers (Corning Inc., Corning, NY, USA). A total of 4 × 104 CAL27 or HSC4 cells were added into the upper chamber in serum-free medium. The bottom chamber was added with medium supplemented with 10% serum. Subsequently, cells that had invaded onto lower surface of the membrane were stained with crystal violet after 24 hours of incubation.

Immunofluorescence staining

We performed the immunofluorescence staining as described previously. All images were viewed on an Olympus IX71 fluorescent microscope and recorded on the screen using an Olympus DP71 camera (Olympus Optical Co. Ltd, Tokyo, Japan) under × 10 magnification.

Western blotting

Western blotting was accomplished as described previously. The antibodies used were LIMA1 (ab154530; 1:1000), phosphatidylinositol-4,5-bisphosphate 3-kinase (PI3K) (ab32089; 1:1000), phosphorylated (p)-PI3K (ab182651; 1:1000), signal transducer and activator of transcription 3 (STAT3) (ab68153; 1:1000), p-STAT3 (ab267373; 1:1000), protein kinase B (AKT) (ab8805; 1:500), p-AKT (ab38449; 1:500), Claudin-1 (ab211737; 1:2000), E-cadherin (ab1416; 1:50), fibronectin (ab268020; 1:1000), N-cadherin (ab76011; 1:5000), Vimentin (ab92547; 1:1000), zona occludens 1 (ZO-1) (ab276131; 1:1000), glyceraldehyde-3-phosphate dehydrogenase (GAPDH, control) (ab8245; 1:1000), and Twist family BHLH transcription factor 1 (TWIST1) (ab175430; 1:1000). All antibodies were obtained from Abcam (Cambridge, MA, USA)

Statistical analysis

All statistical analysis was performed under default option as described by the web resources. The method used for differential analysis in GEPIA was 1-way analysis of variance (ANOVA), using the disease state (Tumor or Normal) as the variable for calculating differential expression, with a P value threshold set at <.01. A χ2 test or Wilcoxon test was used to analyze clinical attributes. The significance of the difference was estimated using Student’s t-test in UALACN. A Spearman or Pearson coefficient test was applied for correlation analysis, with P < .05 considered as significant. Gene sets correlated with LIMA1 were calculated using a Spearman coefficient test with functional gene set enriched results with both P value and FDR < 0.05. All analytical methods involving R packages were performed using R software version v4.0.3 (The R Foundation for Statistical Computing, 2020). All experimental data were statistically analyzed with GraphPad Prism 9.0 software (GraphPad Inc., La Jolla, CA, USA), and the significance of the differences was analyzed using Student’s t-test or ANOVA. In this study, P < .05 was considered statistically significant. All results are shown as the mean ± SD.

Results

A pan-cancer analysis reveals the remarkable role of LIMA1 in HNSC tumor aggressiveness

A pan-cancer analysis was performed to preliminarily explore the expression of LIMA1 and cancer cell functional states in pan-cancers through CancerSEA. Tumor malignant behaviors, such as metastasis, hypoxia, angiogenesis, and EMT, were remarkably associated with higher LIMA1 expression, whereas the cell cycle, DNA repair, and stemness correlated with lower LIMA1 expression in HNSC (Figure 1A and B).

Figure 1.

The pan-cancer and proteomic analysis of LIMA1 in HNSC. LIMA1 associated tumor cell functional states in pan-cancer (A, B). LIMA1 expression and immune cell infiltration level in HNSC in pan-cancer (C). The medium expression of LIMA1 in lymph node of a 62-year-old male patient with HNSC, Patient ID: 1743 (D). The medium expression of LIMA1 in skeletal muscle of a 51-year-old male patient with HNSC, Patient ID: 2608 (E). The medium expression of LIMA1 in oral tissue of a 49-year-old male patient with HNSC, Patient ID: 4427 (F). The high expression of LIMA1 protein in tumor tissue of a 95-year-old female patient with HNSC, Patient ID: 3916 (G). HNSC indicates head and neck squamous cell carcinoma.

TIMER 2.0 estimated the correlation between LIMA1 with tumor immune cell infiltration levels in pan-cancer scales, in which poor B-cell infiltration was observed to be associated with LIMA1 in HNSC (marked as red) (Figure 1C).

Overexpression of LIMA1 in HNSC is associated with TP53 mutant status and human papillomavirus negative status

According to the IHC staining images obtained from The Human Protein Atlas, the proteomic characteristic of LIMA1 in HNSC was analyzed, and moderate or strong LIMA1 staining was observed in tumor tissues (Figure 1D to G).

To determine the differential expression of LIMA1 mRNA in HNSC and normal tissue, the TCGA cohort (TCGA tumors vs TCGA normal match GTEx normal) was assessed. LIMA1 expression was significantly elevated in HNSC tumors compared with that in normal tissues (Figure 2A, P < .01).

Figure 2.

The expression and methylation profiles of LIMA1 in HNSC. The mRNA level of LIMA1 in HNSC and paired normal samples (A). K-M (Kaplan-Meier) survival plot of LIMA1 methylation in HNSC (B). LIMA1 expression level negatively correlated with methylation level (C). Methylation data of CpG dinucleotides and its correlation with LIMA1 expression (D). The relationship between LIMA1 expression and tumor subgroups based on HPV status and TP53 mutation status in HNSC (E). The relative promoter methylation level of primary tumor and paired normal tissues (F). The relative promoter methylation level of tumor subgroups based on TP53 mutation status (G). HNSC indicates head and neck squamous cell carcinoma; HPV, human papillomavirus; TCGA, The Cancer Genome Atlas.

To further estimate the correlation between LIMA1 and the clinical characteristics of HNSC, a χ2 test or Wilcoxon test were employed to reveal LIMA1 expression correlation with multiple clinical parameters (Table 1). No significance was found for age, clinical stage, and tumor pathologic stage; however, sex was significant (P < .05). Additional analysis was performed in UALCAN. LIMA1 was significantly associated with TP53 mutation (P = .002 < .05) and human papillomavirus (HPV) status (P = .003 < .05). Higher expression of LIMA1 was discovered in HPV-negative and TP53 mutant tumors (Figure 2E). As a tumor suppressor gene, TP53 mutation has demonstrated its oncogenic effect in multiple cancer types. An increase in the TP53 mutation level might silence the TP53 function as a consequence of carcinogenesis. HPV is involved in up to 25% of HNSC cases, in which HPV positivity is closely related to a significantly longer survival. Collectively, the expression of LIMA1 was associated with oncogenesis and indicated poor prognosis.

Table 1.

Clinical characteristics and LIMA1 mRNA expression in HNSC.

	Low expression of LIMA1	High expression of LIMA1	P
n	250	250
T stage, n (%)			.791
T1	15 (3.1%)	18 (3.7%)
T2	72 (14.8%)	71 (14.6%)
T3	69 (14.2%)	61 (12.6%)
T4	86 (17.7%)	93 (19.2%)
N stage, n (%)			.237
N0	112 (23.4%)	127 (26.6%)
N1	37 (7.7%)	43 (9%)
N2	86 (18%)	66 (13.8%)
N3	4 (0.8%)	3 (0.6%)
M stage, n (%)			.372
M0	235 (49.5%)	235 (49.5%)
M1	1 (0.2%)	4 (0.8%)
Clinical stage, n (%)			.078
Stage I	8 (1.6%)	11 (2.3%)
Stage II	37 (7.6%)	58 (11.9%)
Stage III	56 (11.5%)	46 (9.5%)
Stage IV	142 (29.2%)	128 (26.3%)
Radiation therapy, n (%)			<.001
No	57 (13%)	96 (21.9%)
Yes	166 (37.8%)	120 (27.3%)
Primary therapy outcome, n (%)			.174
PD	15 (3.6%)	26 (6.2%)
SD	4 (1%)	2 (0.5%)
PR	2 (0.5%)	4 (1%)
CR	189 (45.4%)	174 (41.8%)
Race, n (%)			.454
Asian	3 (0.6%)	7 (1.4%)
Black or African American	23 (4.8%)	24 (5%)
White	216 (44.7%)	210 (43.5%)
Sex, n (%)			.001
Female	50 (10%)	83 (16.6%)
Male	200 (40%)	167 (33.4%)
Age, n (%)			.194
⩽60	130 (26.1%)	114 (22.8%)
>60	120 (24%)	135 (27.1%)
Histologic grade, n (%)			.257
G1	33 (6.9%)	28 (5.8%)
G2	138 (28.7%)	161 (33.5%)
G3	62 (12.9%)	57 (11.9%)
G4	2 (0.4%)	0 (0%)
Anatomic neoplasm subdivision, n (%)			<.001
Alveolar Ridge	10 (2%)	8 (1.6%)
Base of tongue	16 (3.2%)	7 (1.4%)
Buccal mucosa	6 (1.2%)	16 (3.2%)
Floor of mouth	34 (6.8%)	26 (5.2%)
Hard palate	0 (0%)	7 (1.4%)
Hypopharynx	7 (1.4%)	3 (0.6%)
Larynx	62 (12.4%)	49 (9.8%)
Lip	2 (0.4%)	1 (0.2%)
Oral cavity	26 (5.2%)	46 (9.2%)
Oral tongue	48 (9.6%)	78 (15.6%)
Oropharynx	7 (1.4%)	2 (0.4%)
Tonsil	32 (6.4%)	7 (1.4%)
Smoker, n (%)			.195
No	49 (10%)	62 (12.7%)
Yes	196 (40%)	183 (37.3%)
Alcohol history, n (%)			.385
No	74 (15.1%)	83 (17%)
Yes	172 (35.2%)	160 (32.7%)
Lymphovascular invasion, n (%)			.782
No	101 (29.8%)	118 (34.8%)
Yes	58 (17.1%)	62 (18.3%)
Lymphnode neck dissection, n (%)			.193
No	51 (10.3%)	39 (7.8%)
Yes	197 (39.6%)	210 (42.3%)
OS event, n (%)			<.001
Alive	168 (33.6%)	115 (23%)
Dead	82 (16.4%)	135 (27%)
DSS event, n (%)			<.001
Alive	192 (40.4%)	154 (32.4%)
Dead	48 (10.1%)	81 (17.1%)
PFI event, n (%)			.027
Alive	166 (33.2%)	141 (28.2%)
Dead	84 (16.8%)	109 (21.8%)
Age, median (IQR)	60 (53-67)	62 (54-72)	.026

Abbreviations: IQR, interquartile range; OS, overall survival; CR, complete remission/response; PR, partial remission/response; SD, stable disease; PD, progressive disease; OS, overall survival; DSS, disease specific survival; PFI, progression free interval.

DNA demethylation accounts for the increase in LIMA1 expression in HNSC

As common genomic modification, DNA methylation regulates gene expression via de novo DNA methylation and demethylation to adjust the binding of transcription factor(s) to DNA or by directly recruiting proteins involved in gene repression. Study of the methylation profile of LIMA1 might help to understand the underlying mechanism of LIMA1 overexpression in HNSC. We investigated the methylation and expression data of LIMA1, and mRNA expression was inhibited by DNA methylation (Spearman r = −0.28, Pearson r = −0.26, P < .001) (Figure 2C). Subsequently, the entire DNA methylation status of LIMA1 was explored in MEXPRESS. The relative methylation level (β-value) of single CpG parameters and its correlation with mRNA expression was calculated (*P < .05, **P < .01, ***P < .001), in which a close association of mRNA expression with DNA methylation was displayed (Figure 2D). The single CpG parameters in the promoter region were significantly correlated with the expression of LIMA1. The total relative methylation level (β-value) of DNA promoter region in different groups (tumor vs normal) was compared using UALCAN. The LIMA1 promoter was hypo-methylated in HNSC tumors compared with that in normal tissues (Figure 2F, P < .001). The demethylation of the LIMA1 promoter might account for the upregulated LIMA1 mRNA expression in patients with HNSC. Moreover, the promoter relative methylation level (β-value) of TP53 mutant versus TP53 non-mutant samples proved that a high promoter relative methylation level was related to TP53 non-mutant status. Additional prognosis analysis was carried out using SurvivalMeth. Survival analysis using a Kaplan-Meier plot was adopted to identify the prognostic value of LIMA1 methylation (Figure 2B). LIMA1 hyper-methylation was associated with longer survival, while hypomethylation was associated with a poor prognosis (P = .04686, 95% CI: 0.73712 [0.54571-0.99568]).

LIMA1 is an independent prognostic predictor in HNSC

Considering LIMA1 was associated with tumor malignant behavior in HNSC, we further validated the prognostic value of LIMA1. The prognostic significance was estimated in the TCGA cohort. LIMA1High indicated poorer survival (P = .0002, 95% CI, 1.265-2.186) (Figure 3A). Univariate and multivariate Cox regression analysis indicated that LIMA1 was an independent prognostic predictor in patients with HNSC. A modeled nomogram was constructed that could predict the 1-, 3-, and 5-year overall survival of patients with HNSC (Figure 3E).

Figure 3.

Prognostic analysis of LIMA1 in the TCGA set. Prognostic value of LIMA1 in HNSC. Left part (from top to bottom) shows the LIMA1 expression, survival time, and survival status of TCGA HNSC dataset. Rest part (from top to bottom) shows the Kaplan-Meier survival plot of LIMA1 and time-dependent ROC (receiver operating characteristic) curve of LIMA1 (A). Univariate and multivariate Cox regression of P value, risk coefficient (HR), and confidence interval (CI) of gene expression and clinical characteristics (B and C). Nomogram to predict the 1-, 2-, and 3-year overall survival of HNSC patients (D). Calibration curve for the overall survival (OS) nomogram model in the discovery group. A dashed diagonal line represents the ideal nomogram, and the blue, red, and orange lines represent the 1-, 2-, and 3-year observed nomograms (E). HNSC indicates head and neck squamous cell carcinoma; HR, hazard ratio; TCGA, The Cancer Genome Atlas pathological TNM system, the combinations of this system are grouped in five less-detailed stages: Stage 0, Stage I, Stage II, and Stage III,Stage IV.

Genomic CNV and mutation correlated with LIMA1 alteration in patients with HNSC

We identified genomic (CNV) and mutational events associated with LIMA1 alteration. The deletion and amplification (Figure 4A to C) data that correlated with LIMA1 alteration in TCGA HNSC were obtained and visualized using cBioportal. The Volcano map and Scatter plot illustrated the CNV events associated with the LIMA1 altered or unaltered groups (Figure 4A and B). The top remarkable 10 amplified or deleted genes were loaded into GeneMANIA to construct the PPI network. The deletion group was enriched in the cell cycle, whereas the amplification group was related to ribosome and protein targeting to the membrane (Figure 4C). These findings were consistent with corresponding tumor single functional status in the single-cell sequencing analysis.

Figure 4.

Genomic copy number variation (CNV) events and somatic mutation associated with LIMA1 in HNSC. The related deletion or amplification gene sets associated with LIMA1 with representative histogram of top 10 genes (A and B). Gene enrichment analysis of top gene sets conducted by GeneMANIA (C). Lollipop charts display the mutation distribution and protein domains for LIMA1 in HNSC with the labeled recurrent hotspots. Somatic mutation rate and transcript names are indicated by plot title and subtitle, respectively (D). Somatic landscape of TCGA HNSC cohort. The barplot above the legend exhibits the number of mutation burden. Mutation information of each gene in each sample (ordered by LIMA1 expression) is shown by the waterfall plot, in which genes are ordered by their mutation frequencies. A –log10-transformed q values estimated by MutSigCV are shown by side (E). Cohort summary plot displays the distribution of variants with corresponding variant classification, type, and SNV class. Bottom part (from left to right) indicates mutation load for each sample, variant classification type. A stacked barplot shows top 10 mutated genes (F). CNV indicates copy number variation; HNSC, head and neck squamous cell carcinoma; SNV, single nucleotide variation; TCGA, The Cancer Genome Atlas.

The mutational events that correlated with LIMA1 expression are summarized in the Lollipop plot shown as Figure 4D. Missense mutations were the dominant mutation types of LIMA1 and the somatic mutation rate was 0.79%. Figure 4E and F displays the somatic landscape of LIMA1 in the HNSC cohort: The waterfall plot summarizes the gene mutation information, including mutation types (annotated in different colors) and the mutation burden. A cohort summary plot displays the distribution of variants with annotations for variant classification, type, and single nucleotide variation (SNV) class. The mutation load of each sample is shown in the bottom part (from left to right). The somatic mutation profiles between LIMA1High and LIMA1Low were similar. The top 10 mutated genes associated with LIMA1 alteration are listed in a bar chart, including TP53, PIK3CA, CDKN2A, LRP1B, SYNE1, CSMD3, and FAT1, and have been demonstrated to be mutated in the oncogenesis process of tumors.[42 -47] Among these genes, CDKN2A exhibited additional Splice_Site mutation and total mutation rate in the LIMA1High group, which might drive oncogenesis.

Tumor-associated pathways enrichment analysis demonstrates that LIMA1-correlated gene sets are involved in PI3K-AKT and JAK-STAT signaling pathway in HNSC

The LinkFinder module of the LinkedOmics database was used to excavate the co-expression genes of LIMA1 in HNSC. A volcano plot divided the related genes into a positive correlation group (red) and a negative correlation group (green) (Figure 5A). The lateral heat-maps display the expression levels of the top 50 related genes for each sample in the positive correlation group and the negative correlation group, respectively. A bar chart showed the top gene sets that maximized gene coverage with their normalized enrichment score and significance. A P value and FDR < 0.05 was considered as statistically significant (solid bar) while a void bar indicated no statistical significance. The GSEA provided by LinkInterpreter indicated that the related genes were enriched in BPs, such as positive regulation of cell motility and extracellular structure organization; cellular components (CCs), such as cell-substrate junction, cell-cell junction, and extracellular matrix; and molecular functions (MFs), such as cell adhesion molecule binding. The KEGG pathway analysis reveals that the PI3K-AKT and Janus kinase (JAK)-STAT signaling pathway enriched in gene sets correlated with LIMA1 (Figure 5C). The results were sorted by the weighted set cover algorithm to display top gene sets while maximizing gene coverage.

Figure 5.

LIMA1 correlated gene set enrichment analysis (GSEA). The genes correlated with LIMA1 were analyzed by Pearson test in HNSC (A). The top 50 genes positively or negatively correlated with LIMA1 in HNSC (B). GO (BP, biological process; CC, cellular components; MF, molecular functions) and KEGG enrichment analysis, with P value, FDR < .05 was considered as significant (C). FDR indicates false discovery rate; GO, gene ontology; HNSC, head and neck squamous cell carcinoma; KEGG, Kyoto Encyclopedia of Genes and Genomes.

The PI3K/AKT and JAK/STAT signaling pathways with annotation are shown in Figure S1 in the Supplemental material.

LIMA1 knockdown inhibits the EMT progression and tumor-associated signaling pathways in HNSC

On the basis of bioinformatic analysis results, further verification was performed to gain an insight into the function of LIMA1. CAL27 and HSC4 cells were cultured and transfected with LIMA1 siRNA. Cell migration assays and wound healing assays illustrated that silencing LIMA1 suppressed cell migration. In addition, the protein level of tumor-associated pathway markers were measured using immunoblotting to further verify the tumor-associated signaling pathways predicted by GSEA. The protein levels of AKT, p-AKT, and p-STAT3 were significantly downregulated in the transfected CAL27 cells, while transfected HSC4 cells showed significantly decreased levels of p-STAT3 and AKT. LIMA1 knockdown’s effect on EMT was validated by detecting the protein levels of EMT-related proteins. The results revealed that N-cadherin and TWIST1 levels were decreased, whereas E-cadherin and Claudin 1 levels were increased after LIMA1 silencing in CAL27 cells, whereas silencing of LIMA1 in HSC4 cells significantly decreased the level of N-cadherin, and increased the level of E-cadherin. (Figure 6A to D). A LIMA1-inhibition treated group and a control group were designed for comparison. The relative N-cadherin level decreased and the E-cadherin level increased in the treated group compared with that in the control group (Figure 6E).

Figure 6.

Silencing LIMA1 suppresses EMT process and tumor-associated pathway in vitro. Representative images and analysis of transwell assay and wound healing assay (A and B). The expression level of proteins in the LIMA1-associated signaling pathways was determined by Western blotting assay (C). Western blot analysis of EMT-related proteins in CAL27 and HSC4 cells with LIMA1 knockdown (D). Representative images of immunofluorescent staining of N-cadherin or E-cadherin (E). (*P < .05; **P < .01; ***P < .001; ****P < .0001). EMT indicates epithelial-mesenchymal transition.

Discussion

In this study, multiple layers of evidence demonstrated that LIMA1 overexpression is associated with tumor aggressiveness in HNSC for tumor metastasis, hypoxia, angiogenesis, and EMT. Hypoxia is recognized as an inhibitor that modulates anti-tumor immunity. DNA demethylation of the promoter region of LIMA1 might explain the differential expression level of LIMA1 in HNSC. Human papilloma virus (HPV)-negative HNSC and TP53 mutation were found in the high LIMA1 expression group. Compared with HPV-negative HNSC, HPV-positive HNSC represents a distinct subset that responds better to concurrent chemo-radiation and is related to better prognosis, with distinct genomic features.[49,50] High LIMA1 expression might indicate an adverse outcome in patients with HNSC. A previous study demonstrated that TP53 mutation was capable of regulating LIMA1 at transcription level, and TP53 mutations account for the downregulation of LIMA1 associated with poor survival of patients with cancer. Genomic alterations, including CNVs and mutations, were closely associated with LIMA1, which, together with the cell cycle, ribosome, protein targeting to the membrane, and tumor suppressive gene mutations, might contribute to oncogenesis in patients with HNSC. Interestingly, the cell cycle is believed to modulate the drug resistance of tumors. LIMA1 might serve as a therapeutic biomarker of HNSC. Further functional analysis corroborated that LIMA1, as an oncogene, participates in the progression of HNSC. GO and KEGG pathway enrichment analyses implied that LIMA1 co-expressed gene sets were enriched in the PI3K-AKT and JAK-STAT signaling pathways. The PI3K/AKT downstream pathway was proven to participate in anti-apoptosis processes by targeting apoptosis-related proteins such as MCL1 and PIM1, and cell cycle progression or inhibition. A study has demonstrated that suppressing the PI3K/AKT signaling pathway induces cell cycle arrest (G2/M phase) and apoptosis in HNSC, which serves as a promising therapeutic target for patients with HNSCs. Besides, the cell cycle participates in the mechanism of drug resistance of tumors, which suggest a potential therapeutic value. The AKT activation might account for the cell invasion and metastasis of HNSC. Studies have demonstrated that PI3K inhibitors represent a therapeutic strategy against HNSC.[53 -55] Moreover, activation of JAK-STAT signaling contributes to cellular invasion and proliferation in HNSC.[56,57] This highlights the tumor aggression-associated role of LIMA1 in patients with HNSC. However, studies in other types of cancers showed that the loss of LIMA1 contributes to tumor cell migration and cultivate metastasis.[10,11,13,15] Despite its controversial expression pattern and distinct functions in HNSC compared with those in other types of cancer, LIMA1 could be a unique biomarker for HNSC. Interestingly, the 2 signaling pathways are closely associated with the process of EMT. The PI3K/AKT signaling pathway directly induces EMT. Activation of the JAK-STAT pathway enhances the EMT process, leading to strengthening of tumorigenic and metastatic behavior, cancer stem cell (CSC) transition, and chemoresistance in cancer. Supplementary in vitro assays further validated that LIMA1 participates in tumor-associated signaling pathways and markedly affects the EMT process.

The tumor microenvironment is widely implicated in tumorigenesis and progression.[60,61] Through data mining, LIMA1 was found to correlate negatively with B-cell infiltration. A lack of infiltrating B-cells might mediate the immune response failure in tumors. Deletion of the B-cell-associated immune response might account for the anti-tumor immune failure in HNSC as a consequence of increased LIMA1 expression.

A recent study illustrated the role of LIMA1 in cellular metabolism, in which LIMA1 depletion contributed to reduced mitochondrial adenosine triphosphate (ATP) production, which affected tumor growth. Functional mitochondria can regulate the intracellular reactive oxygen species level and promote EMT in tumor cells. These studies suggested that the functional state of mitochondria might be related to the promotion of EMT by LIMA1. Besides, EMT was demonstrated to affect the response to immunotherapy because of its immunomodulatory crosstalk with pan-cancer immune evasion. A recent study showed that escape of tumor cells from PD-1 blockade therapy is related to attenuated mitochondrial activity in T-cells. Attenuated mitochondrial activity results in a decrease in the c-ATP level, and the c-ATP level of tumor-infiltrating lymphocytes was negatively associated with PD1 expression on T-cells. Furthermore, enhanced mitochondrial function can improve the efficiency of anti-PD1 therapy. Hence, mitochondrial function might directly or indirectly (through EMT) affect anti-tumor immunity and the immunotherapy response, which indicates a potential complex molecular mechanism by which LIMA1 mediates EMT in HNSC.

A limitation of this study might be the lack of clinicopathologic data to support the conclusions.

Conclusions

LIMA1 is overexpressed in HNSC and contributes to the progression of EMT in HNSC, which further drives tumor invasion and metastasis. LIMA1 acts as a prognostic biomarker and indicates poor survival in HNSC. Silencing of LIMA1 suppressed the EMT process and the tumor invasion and metastasis in vitro. The expression of LIMA1 is affected by its methylation level and is associated with genomic alteration in HNSC, in which tumor-associated pathways are activated.

Click here for additional data file.

Supplemental material, sj-png-1-onc-10.1177_11795549221109493 for Overexpression of LIMA1 Indicates Poor Prognosis and Promotes Epithelial-Mesenchymal Transition in Head and Neck Squamous Cell Carcinoma by Wei Ma, Yiqun Liao, Ziwen Gao, Wenyan Zhu, Jianbing Liu and Wandong She in Clinical Medicine Insights: Oncology