Firstly, API-related differentially expressed genes (DEGs) in squamous cell carcinoma-9 (SCC-9) cells were identified through RNA-seq and a preliminary analysis was conducted on the molecular mechanism of possible action of Apigenin. Next, combined with relevant individual information of TSCC patients from The Cancer Genome Atlas (TCGA) databases, 7 API-related genes were found to create a prognostic model and to survey the underlying mechanisms of TSCC patients [25]. The low- and high-risk groups were constructed to estimate patients’ outcomes by accurately predicting their overall survival (OS) from TSCC. Furthermore, this model was used to classify the tumor microenvironment (TME) and drug sensitivity of TSCC patients in different risk groups. A schematic of the analysis workflow is demonstrated (Fig. 1).

Fig. 1.

Analysis workflow. TSCC, tongue squamous cell carcinoma; API, apigenin; DEGs, differentially expressed genes; GSEA, gene set enrichment analysis; TME, tumor microenvironment.

The TSCC cell lines SCC-9 and CAL-27 were purchased from the American Type Culture Collection (ATCC). For culturing, SCC-9 cells were grown in DMEM-F12 maintained at 37 °C with 5% CO${}_2$. And CAL-27 cells were cultured in DMEM medium in the same environment. All cell lines were validated by STR profiling and tested negative for mycoplasma. For experimental studies, SCC-9 and CAL-27 cells were treated with or without API (IA0400, Solarbio, Beijing, China) for the indicated period with dimethyl sulfoxide (DMSO, D8371, Solarbio, Beijing, China) as the vehicle. None of the final concentrations of DMSO in the medium exceeded 0.1% in all treatments.

The semi-inhibitory concentration (IC50) of API was calculated by the Cell Counting Kit-8 (CCK-8) assay (BA00208, Bioss, Beijing, China). In brief, 100 µL of SCC-9 and CAL-27 cells growing exponentially at a concentration of 5 $\times$ 10${}^4$ cells/mL were seeded into 96-well plates and further incubated with different concentrations of API (0, 2, 4, 8, 16, 32, 64, 128 µmol/L). After 24 hours, the CCK-8 chemical agent (10 µL) was mixed in each well, and the mixtures were incubated at 37 °C for 120 minutes. The OD values were tested at a wavelength of 450 nm. The IC50 was determined by Prism 9 software (GraphPad Software, San Dieg, USA).

SCC-9 and CAL-27 cells (3 $\times$ 10${}^3$) in 100 µL of medium were inoculated into 96-well plates. Following the CCK-8 assay procedure, cell proliferation was analyzed at 1, 2, 3, 4, and 5 days. At the specified point in time, 10 µL of CCK-8 reagent was added to every well. After incubating for 2 h at 37 °C, a wavelength of 450 nm was used to analyze the OD values in each well, and growth curves were created to define the growth rates.

Using a 200 µL pipette tip, a scratch was made in the center of each well. Each well was washed 3 times with PBS to remove cell debris. Then, the residual cells were treated with or without 20 µM API for 24 h. The migrating cells were imaged at 0 or 24 h by using a Leica inverted microscope (DMI8, Leica, Solms, Germany; magnification, 50$\times$).

In the invasion assay, transwells were pretreated with matrigel (1:8). Following the suspension of SCC-9 and CAL-27 cells in serum-free medium (10 $\times$ 10${}^5$/mL), 200 µL of the cell suspension containing API (0, 5, 10 µM) were added to the upper compartment, and 800 µL of medium containing 20% fetal bovine serum (FBS) were used in the lower chamber. Then, SCC-9 cells were incubated for 24 hours in migration and 48 hours in invasion. As for CAL-27, cells were incubated for 48 hours in migration and 72 hours in invasion. Cells that had migrated were then fixed with 4% paraformaldehyde, followed by crystal violet staining. The number of migrating cells was determined after the cell chamber was photographed using a microscope (DMI8, Leica, Solms, Germany; magnification, 50$\times$).

Six-well plates were used to cultivate SCC-9 and CAL-27 cells (1000/well). 20 µM API was applied to the cells. And every three days, the culture medium was changed. After 10 days of incubation, cells were immobilized with 4% paraformaldehyde. Then, the cells were followed by 10 minutes of crystal violet staining. A digital camera was used to record colony images, and colonies with more than 50 cells each were tallied.

Sections were cut from formalin-fixed, paraffin-embedded adjacent noncancerous tissues and cancer tissues. The sections were routinely deparaffinized, rehydrated, antigen-retrieved, and blocked with 3% hydrogen peroxide. Thereafter, primary antibodies for ITGA5 (1:200, Proteintech) were incubated on the sections overnight at 4 °C, followed by secondary antibodies that were HRP-labeled. The sections were stained with hematoxylin and then with diaminobenzidine. Unaware of any clinical information about our patients, two skilled pathologists independently evaluated and scored the stained sections based on the degree of staining and the ratio of tumor cells with a positive stain. The staining intensity has been separated into four groups: 0 for no staining, 1 for mild staining, 2 for moderate staining, and 3 for intense staining. The proportion of positively stained tumor cells was classified depending on the following criteria. 0 (5% positive cancer cells), 1 (5–25% positive cancer cells), 2 (26–50% positive cancer cells), 3 ($>$50% positive cancer cells). The total score is derived from the sum of the staining intensity and the stained cell score. Based on the total score, patients were separated into two groups: 0–2 negative and 3–6 positive.

Transcriptome RNA sequencing data from SCC-9 cells treated with or without 20 µM API for 24 hours were provided by Applied Protein Technology. RNA-seq data of TSCC patients obtained from the TCGA contained 148 TSCC samples and 15 normal tissue samples. Clinical data in 147 samples, including age, sex, tumor grade, and survival status, were also gathered from the TCGA. Besides, data from 58 additional samples were downloaded on the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/) (ID: GSE111395, GSE75540), which later was used a as test cohort.

The “limma” package, containing p as well as false discovery rate (FDR) 0.05 as the cutoff points, was employed to detect DEGs between tumor tissues and normal tissues acquired from TCGA datasets [26]. This “limma” package was also used to analyze DEGs in SCC-9 cells treated with or without 20 µM and prognosis-related genes of TSCC patients. Next, API-related DEGs were found in both tissue DEGs, and prognostic genes were identified. Furthermore, API-related genes with prognostic value were selected by univariate Cox analysis of overall survival (OS). Then, interacting networks based on these API-related genes were created via the STRING database (version 11.5, developed jointly by the Swiss Institute of Bioinformatics, the Novo Nordisk Foundation Protein Research Center and the European Molecular Biology Laboratory). Aiming to create a prognostic model, Lasso Cox regression analysis further revealed target genes with the “glmnet” package. In addition, multivariate Cox analysis was also carried out to evaluate the prognostic sign of the genes. While calculating the penalty parameter ($\lambda$) using the smallest possible criterion, normalized expression of the genes in TSCC patients and coefficients were gathered to analyze total risk scores. The algorithm was as follows: risk score = sum (each gene’s expression in each TSCC patient $\times$ corresponding coefficient). Then, TSCC patients in the TCGA cohort were allocated to the high- or low-risk group as well as the Gene Expression Omnibus (GEO) cohort [27]. Principal component analysis (PCA), from running the “stats” R package, was conducted to estimate this model’s ability to isolate samples into two groups [28]. In addition, the “Rtsne” R package was run to create t-distributed stochastic neighbor embedding (t-SNE) and then further determine the groups’ distribution. For the survival analysis of each gene, the optimal cut-off expression value was determined by the “surv_cutpoint” function of the “survminer” R package. Kaplan–Meier (K-M) analysis by the “survminer” package was built to predict the OS of TSCC patients between the high- and low-risk groups.

Receiver operating characteristic (ROC) curve analysis was performed using the “survival”, “survminer”, and “timeROC” packages. The 1-year, 3-year, and 5-year OS were determined by ROC curves to estimate the anticipation capacity of the model. The relationship between clinical data acquired from the TCGA and API-related gene expression was analyzed with Cox regression analyses. Moreover, the results were visualized through forest maps.

The input data were split into high- and low-risk groups obtained by risk scoring, and the underlying molecular pathways in API-related genes were further analyzed. Gene set enrichment analysis (GSEA) was performed by GSEA software (version 4.3.2, Broad Institute, UC San Diego, CA, USA). Simulated value = 1000, p-value 0.05, which was conducted as the criterion for statistical significance.

Each sample in the high- or low-risk group was analyzed to compute the number of 22 types of immune cells infiltrating the tumors. Furthermore, the CIBERSORT method was used to assess the correlation between patient risk scores and immune cell density [25]. The risk score was determined to explore the connection between API-related genes and 22 levels of immune cell infiltration. Additionally, between the high- and low-risk groups, the immune status and stromal scores were evaluated by the ESTIMATE algorithm [29].

To explore drug sensitivity in the low- and high-risk groups, IC50 values of TSCC treated with multidrug treatment were determined through the “prophetic” software package.

Based on the Human Protein Atlas, each API-related DEG used to build a prognostic model was subjected to immunohistochemistry imaging to further confirm the protein expression of these genes.

Total RNA was obtained from SCC-9 cells treated with or without 20 µM API via TRIzol (Invitrogen, Carlsbad, CA, USA). Then, complementary DNA (cDNA) was generated using the Verso complementary DNA synthesis kit (Invitrogen, Carlsbad, CA, USA). Sequences of primers for detecting target genes are described in Table 1.

During data collection, gene expression was calculated through Student’s t-test. Survival analysis was determined using the R survival package. In addition, the log-rank test was used to test the survival rate of each group. The Kaplan‒Meier curve was utilized to predict the survival probability of the high- and low-risk groups in the data. Independent forecast factors of OS were employed through univariate and multivariate Cox regression. Moreover, the correlation of the two subgroups in the tumor immune microenvironment, including immune cells and immune pathways, was assessed employing the Mann–Whitney test. Additionally, the correlation coefficient was measured by Spearman analysis. R version 4.1.0. (The R Foundation for Statistical Computing, Boston, MA, USA) was used for all statistical analyses, and the p-value was 0.05.

To determine the API (Fig. 2A) IC50 value of the TSCC cell line, SCC-9, and CAL-27 cells were treated with API at progressively higher concentrations (0, 2, 4, 8, 16, 32, 64, 128 µM) for 24 hours (Fig. 2B). As shown in Fig. 2B, the IC50 value was approximately 35.16 µM in SCC-9 cells and 39.21 µM in CAL-27 cells based on the CCK8 assay. Cell proliferation was also evaluated using the CCK8 assay (Fig. 2C,D). The effect of API on the migratory capacity of SCC-9 and CAL-27 cells was then clarified by wound-healing assay. Our data showed that API could decrease gap closure at concentrations of 20 µM (Fig. 2E,F). Meanwhile, we also verified that API could effectively suppress cell migration and invasion by transwell assay (Fig. 2G,H). Furthermore, the colony-formation assays demonstrated that API could effectively inhibit clonogenic capacity in SCC-9 and CAL-27 cells (Fig. 2I,J). Our results show that API could markedly suppress cell proliferation, migration, and clonogenic ability of SCC-9 and CAL-27 cells.

Fig. 2.

Apigenin (API) reduced Tongue squamous cell carcinoma (TSCC) cell characteristics. (A) Chemistry architecture of API. (B) The IC${}_{50}$ value of API. The IC${}_{50}$ value was calculated by Prism 6 software. (C,D) API inhibited SCC-9 and CAL-27 cell proliferation. (E,F) API inhibited SCC-9 and CAL-27 cell migration. Typical images observed from the wound-healing assay of TSCC cells are displayed for 0 h and 24 h. (G,H) The migration and invasion ability of TSCC cells with or without API treatment in transwells. Representative photographs of transwell images were presented. (I,J) The indicated cells were seeded as single cells and maintained for 10 days. Representative photographs of colony formation were presented. The results displayed are representative of three stand-alone experiments. (F), (H), and (J) data are shown as means $\pm$ SD (standard deviation). ns p $>$ 0.05, ***p 0.001.

To explore the underlying molecular mechanism of how API restrains the malignancy of TSCC, we obtained the results of API-treated and untreated groups of DEGs through transcriptome RNA sequencing (Fig. 3A,B). Then, we also analyzed DEGs in RNA-seq data between TSCC samples and normal tissue samples from the TCGA (Fig. 3C,D). The individual information of these TSCC patients was generalized in Table 2. Integrated analyses were performed to demonstrate relationships between DEGs in SCC-9 cells and DEGs in the RNA-seq group or the prognosis-related group with $|$ log2 (fold change) $|$ $\ge$1 and FDR 0.05 (Fig. 3E). As shown in Fig. 3E, some API-related genes were differentially expressed between normal and tumor tissue samples. Nine of them (SLC7A14, ETNPPL, TMEM213, GPR158, GJB7, CASP14, HOXD-AS2, ITGA5, and MIR503HG) were correlated with OS by univariate Cox regression analysis (Fig. 3F). Unexpectedly, SLC7A14 and ETNPPL were downregulated in tumor samples. However, their HR expression was much higher, which can upregulate offset possibility. Consequently, these two genes were ruled out for further study. In total, 7 API-related prognostic DEGs were retained. (all FDR 0.05, Fig. 3G). The correlation between these 7 DEGs is presented in Fig. 3H.

Fig. 3.

Determination of prognostic API-related differentially expressed genes (DEGs). (A,B) Heatmap and volcano plot of API-related DEGs. (C,D) Heatmap and volcano plot of TSCC-related DEGs. (E) Venn diagram to determine API-related DEGs that were associated with overall survival (OS). (F) The expression of 9 overlaying key genes in tissues. (G) Through univariate Cox regression analysis, Forest plots display the relationship between OS and gene expression. (H) The network of prognostic API-related genes (red and blue lines represent positive and negative relativities, respectively). Darker colors represent stronger associations. FC, fold change.

In the Human Protein Atlas database, only 4 genes (CASP14, TMEM213, GPR158, and ITGA5) of this prognostic model have been analyzed immunohistochemically. As displayed in Fig. 4A, the protein expression of these four genes was shown to be different in normal tissues than in tumor tissues (Fig. 4A).

Fig. 4.

Validation of genes in this prognostic model. (A) Immunohistochemistry results of tissue samples from the Human Protein Atlas database (image available from v15.proteinatlas.org). (B) mRNA levels of caspase 14 (CASP14), transmembrane protein 213 (TMEM213), G protein-coupled receptor 158 (GPR158), integrin subunit alpha 5 (ITGA5) in SCC-9 cells. (*p 0.05, **p 0.01 and ***p 0.001, compared with the control group). (C) Representative immunohistochemical images of ITGA5 expression in primary TSCC cancer tissues and normal aracancerous tissue. Scale bar, 50 µM.

To validate these 4 genes in SCC-9 cells, RT-qPCR was utilized to evaluate their mRNA expression. These SCC9 cells were treated with or without 20 µM API for 24 hours. According to the results, TMEM213 and GPR158 mRNA levels were increased in the API group. In comparison with the CTRL group, the mRNA expression of CASP14 and ITGA5 was downregulated (Fig. 4B). In these 4 genes, ITGA5 and CASP14 are oncogenes, and API inhibits their expression. Since HR in ITGA5 is much higher (Fig. 3G), we further confirmed that ITGA5 was enhanced in tumor tissues than in paracancerous tissues through IHC (Fig. 4C).

Then, a prognostic model was constructed from the 7 genes noted above using LASSO Cox regression analysis. The risk score formula is listed below [26]: e (0.070 $\times$ expression level of CASP14 + 0.216 $\times$ expression level of ITGA5 + 0.212 $\times$ expression level of GPR158 + 0.319 $\times$ expression level of MIR503HG + 1.99 $\times$ expression level of ETNPPL + 0.214 $\times$ expression level of GJB7 + 0.309 $\times$ expression level of TMEM213). In accordance with the cut-off value, the TSCC patients in the TCGA cohort were divided into two groups: a low-risk group (n = 74) and a high-risk group (n = 73) (Fig. 5A). In addition, high-risk patients had a higher death rate than low-risk patients (Fig. 5B). Moreover, patients in the high-risk or low-risk group were determined to be located independently via PCA and t-SNE analysis (Fig. 5C,D). The Kaplan‒Meier curve steadfastly demonstrated that patients with TSCC in the low-risk category had notably better OS than those in the high-risk group (Fig. 5E, p 0.001).

Fig. 5.

Establishment and Validation of an API-related prognostic model. (A) The distribution of TSCC patients with risk scores in the The Cancer Genome Atlas (TCGA) cohort and Gene Expression Omnibus (GEO) cohort. (B) The locations and OS status of TSCC patients with different risk scores in the TCGA cohort and GEO cohort. (C) Principal component analysis (PCA) plot of the prognostic model. (D) t-SNE analysis of the prognostic model. (E) Kaplan‒Meier survival curve of TSCC patients. (F) ROC curves demonstrate the predictive validity of risk scores. (G) Forest plots for univariate and multivariate Cox regression analysis of patients in TCGA. (H) Nomogram predicting 1-, 2-, and 3-year OS in TSCC patients. (I) Calibration curves for consistency testing between forecasted and real OS.

Subsequently, to verify the reliability of the model established in the TCGA cohort, the patients from the GEO cohort were also placed in either a high- or low-risk group by the median value obtained as the same formula as that in the TCGA cohort (Fig. 5A). And from Fig. 4B, Patients in the low-risk group were more probably to own a longer OS (Fig. 5B). Similar to the outcome received from the TCGA cohort, PCA and t-SNE analysis revealed that the distribution of patients in two subgroups was in discrete directions (Fig. 5C,D). In addition, patients in the high-risk group had worse survival rates compared to the low-risk group (Fig. 5E).

After eliminating samples with inconclusive clinical information, ROC analysis was utilized to calculate AUCs. As revealed in Fig. 5F, the AUCs of the 7-gene model for 1-, 2-, and 3-year survival were 0.736, 0.718, and 0.709, respectively. The risk score was proven to be a stand-alone forecaster of OS in TSCC patients through univariate and multivariate Cox regression analyses (Fig. 5G). Then, all associated clinical factors for patients in the TCGA were used to build a predictive nomogram for OS (Fig. 5H). The red dot represents the total clinical characteristics score in a TSCC patient from the TCGA, which is 263 points. This red dot also indicates that the survival probabilities in this patient for 1, 3, and 5 years were 85.4, 69.5, and 53.2%, respectively. The 1-, 3-, and 5-year OS of TSCC patients were also evaluated using calibration curves (Fig. 5I).

Aiming to determine the function of API-related prognostic genes, GSEA was utilized to diagnose relative underlying pathways between the low- and high-risk groups. Based on the results, the enriched pathways in the high-risk group were mostly specialized in focal adhesion, ECM receptor interaction, viral myocarditis, regulation of actin cytoskeleton, leishmania infection, and coagulation cascades (Fig. 6A). The enriched mechanisms of the low-risk group were Parkinson’s disease, oxidative phosphorylation, Huntington’s disease, Alzheimer’s disease, seven amino acid metabolism, and one carbon pool by folate (Fig. 6B).

Fig. 6.

Functional analysis based on the prognostic model. (A,B) Mechanisms dramatically enriched in high- (A) and low-risk groups (B) through gene set enrichment analysis (GSEA).

As shown in Fig. 7A, the risk score was linked to regulatory T cells (Tregs), follicular helper T cells, CD8 T cells, resting mast cells, memory B cells, and M1 macrophages (Fig. 7A). In addition, a correlation between the 7 genes and the number of immune cells was also observed. The results showed that some immune cells were easily affected by these 7 genes (Fig. 7B). Additionally, T cells can be inhibited by stromal activation of the TME. The tumor stroma and ESTIMATE scores were much lower in the low-risk group than in the high-risk group (Fig. 7C).

Fig. 7.

Relationship between risk score and tumor microenvironment (TME) or Drug susceptibility. (A) Relations between risk score and immune cell types. (B) Seven genes in the prognostic model were associated with the number of immune cells. (C) Relationships between risk score and TME scores. (D) Drug susceptibility. Data are shown as the mean $\pm$ SD. (*p 0.05, **p 0.01 and ***p 0.001).

A drug susceptibility test was performed to analyze the effect of anticancer drugs in different risk groups. In the high-risk group, the results indicated that the IC50 values of bexarotene, bicalutamide, pazopanib, and WO2009093972 were much lower. The chemotherapeutic agents in the high-risk group, including A443654 and vinorelbine, had higher IC50 values than those in the low-risk group (Fig. 7D).