Ferroptosis and WDFY4 as novel targets for immunotherapy of lung adenocarcinoma

FPI identification and sample grouping

In accordance with the findings of Liu et al. [18], a total of 24 ferroptosis regulator genes (FRGs) were identified and categorized into two groups based on their respective modes of action: positive and negative genes. The positive regulation of gene expression was observed through FPI’s influence on LPCAT3, NCOA4, ACSL4, GPX4, SLC3A2, ALOX15, SLC7A11, NOX5, NFE2L2, NOX3, NOX1, and NOX4. Conversely, the negative regulation of gene expression was evident in the case of FDFT1, COQ10A, HMGCR, and COQ10B, among others. The FPI value was determined by calculating the enrichment scores (ES) of positive regulatory genes minus the ES of negative regulatory genes by the single-sample gene-set enrichment analysis (ssGSEA) algorithm. Subsequently, a survival analysis was conducted on the FPI values in LUAD patients using the x-title algorithm. The results indicated that an FPI < 0.47 is associated with improved survival in LUAD patients, whereas patients with an FPI > 0.47 exhibit a poorer long-term prognosis. The samples were subsequently categorized into high- and low-score groups based on their PFI scores. Variance analysis was conducted using LIMMA, and the differential genes were determined based on the following criteria: log2 fold change ≥1.5 and Benjamin-Hochberg (B-H) adjusted P-value < 0.05.

Immune infiltration and correlation analysis

Studies have confirmed that FPI is closely related to immune regulation, although the precise regulatory mechanism in LUAD remains unknown. In order to investigate this potential relationship further, we utilized ssGSEA to calculate the relative degree of cellular infiltration for each subtype, based on the expression levels of 24 immune cell surface markers as outlined by Finotello et al. [19]. Additionally, we conducted bias correlation analysis to examine the relationship between the relative infiltration values and PFI values of each cell subtype of cells, following standard quantification. The resulting correlation statistics were adjusted by Pearson’s method.

Gene set enrichment analysis (GSEA)

The gene-set-based enrichment analysis method known as GSEA is employed in the analysis of gene expression data. Initially, the purpose is established, followed by the selection of gene sets from MSigDB are selected for systematic analysis. The gene expression data are then ranked based on their association with the utilizing normalized scores. Subsequently, the evaluation of the impact of synergistic changes in genes within the gene sets on phenotypic changes is conducted by determining if the genes within each gene set are enriched for the clinical phenotypic correlation. The calculated parameters for this analysis were as follows: minGSSize = 10, maxGSSize = 500, and B-H adjusted p-value Cutoff = 0.05.

Gene set variation analysis (GSVA) and survival analysis

Gene Set Variation Analysis (GSVA) is a non-parametric, unsupervised technique that distinguishes itself from Gene Set Enrichment Analysis (GSEA) by its ability to compute enrichment scores for individual gene sets in each sample without the need for pre-grouping of samples. In other words, GSVA transforms a gene sample matrix into a gene-set by sample matrix, facilitating the assessment of pathway enrichment for each sample. This approach offers a standardized quantification of gene enrichment results, thereby enabling statistical analysis of differences in subsequent pathway sets. Following that, by employing the LIMMA package for differential expression analysis, which relies on GSVA scores, it becomes feasible to discern genes that exhibit differentially expression across samples. The resultant sets of differentially expressed genes possess greater biological significance and are more amenable to interpretation in the context of disease progression. Moreover, leveraging GSVA scores, an additional prognostic analysis was performed using K-M survival analysis to elucidate the influence of pathways on the long-term survival of patients.

Single-cell RNA clustering

Based on the single-cell data, we first performed quality control on the downloaded raw data to filter out interfering impurities. The criteria employed for this purpose encompassed cellular gene count, mRNA count, ribosomal genes, etc. The preprocessing parameters utilized were as follows: min Gene = 500, max Gene = 4000, and the percentage of mitochondrial genes <5%. Subsequently, principal components analysis (PCA) was employed to downscale highly variable genes, and the resulting specific components analysis (PC) values were further downscaled to a two-dimensional space using linear downscaling techniques such as t-distributed stochastic neighbor embedding (tSNE) or uniform manifold approximation and projection (UMAP), following the selection of an appropriate PCA dimension. After the aforementioned procedure, the majority of feature values in the transcriptional profile were represented in a two-dimensional format. Prior to conducting PCA downscaling, the expression data underwent centering, and cell cycle regression analysis was performed using the Scale Data function.

Cell identification of core subpopulations, pseudotime analysis, and variance analysis

To identify the significantly highly expressed genes within each cluster, we used the Find All Markers function from the Seurat package. This function allowed us to obtain differential genes by configuring the test use parameter and conducting Wilcox differential statistics analysis, with p-value < 0.05 to define differential genes. Following the completion of cell clustering and downscaled visualization we were able to demonstrate the heterogeneity of the cell population. To annotate the cell subpopulations within each cluster, we used the Single R and scCATCH algorithms on the differential genes obtained from each cluster. To elucidate the developmental status and potential regulatory nodes of each subpopulation, we conducted pseudotime analysis on the transcriptional profiles of cell subpopulations based on the Monole algorithm. Subsequently, we examined the pseudotime genes that are relevant to development.

Single cell GSVA and pathway enrichment analysis

We estimated the kernel density of the pathway sets for the transcriptional profile level in a single cell based on the unsupervised GSVA algorithm, which shares similarities with the nonparametric test method used in probability theory to estimate the unknown density functions. Gene sets were compiled based on specific gene sets of interest, such as signaling pathways and GO entries. The RNA-seq count data of genes in each cell were calculated iteratively on the gene expression matrix, allowing for the exploration of enrichment values for specific pathways in each cell. Based on the grouping, cell status and subpopulation information, variance analysis was performed on the GSVA values. The method was the same as described in Method 4.

Identification and functional analysis of core regulatory genes

To probe the potential role of ferroptosis genes in regulating core cell populations during the progression of LUAD, we conducted an analysis that involved the intersection of single-cell pseudotime genes and ferroptosis-related genes. Subsequently, we elucidated the functional changes of these genes using a pseudotime heat map. Additionally, for the core genes, we employed protein-protein interaction (PPI) network analysis to gain further insights into their protein function and interactions. To accomplish this, we initially utilized the String online database (https://www.string-db.org/) was used to analyze the network, followed by network contribution analysis using Cytoscape. The primary methods employed for the interactions included Textmining, Experiments, Databases, Co-expression, Neighborhood, Gene Fusion, and Co-occurrence, with an interaction score threshold of 0.4. Additionally, the biological functions and protein substructure binding thresholds were analyzed using the Toppgene database (https://toppgene.cchmc.org), and the drug targets of core genes by CMAP database. A significance level of B-H adjusted P-value < 0.05 was deemed indicative of a significant difference.

Cell lines and cell culture

Human lung adenocarcinoma A549 cells were obtained from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China) and were cultured in Roswell Park Memorial Institute 1640 (RPMI 1640) with 10% fetal bovine serum (FBS; Gibco Life Technologies, NY, USA), and 1% penicillin–streptomycin (Gibco Life Technologies, NY, US) at 37°C in a humidified incubator with 5% CO₂.

Establishment of WDFY4 overexpressing cell line

The lentivirus carrying the WDFY4 gene sequence was obtained from Hanbio Co. LTD (Shanghai, China). Lung adenocarcinoma A549 cells were seeded in 12-well plates and then infected with lentiviruses following the manufacturer’s recommended protocol. After 24 hours, serum-free medium was replaced with complete medium. To establish a stable cell line with overexpressed WDFY4, lentivirus-infected cells were chosen by incubation with 2 μg/ml puromycin.

Cell counting Kit-8 assay

The 2000 A549 cells transfected with lentivirus (ad-WDFY4) were inoculated in a 96-well plate and cultured in an incubator for a predetermined duration (24, 48, 72 h). After incubation with CCK-8 solution (Beyotime Biotechnology, Shanghai, China) for 2 h, the absorbance (450 nm) was measured by a microplate reader.

Cell apoptosis assay

Ad-WDFY4/A549 cells were cultured for 24 h, followed by digestion with 0.25% trypsin and subsequent harvesting. Subsequently, Annexin V-PE/propidium iodide (PI) apoptosis detection kit (Beyotime Biotechnology, Shanghai, China) was used to evaluate the degree of apoptosis by flow cytometry. The percentage of apoptotic cells was calculated as the sum of the early and late apoptotic cells located in the lower right quadrant and upper right quadrant, respectively.

Scratch wound healing assay

The cells were seeded in 6-well plates and cultured until reaching approximately 80% confluence. A scratch was made in the middle of the culture dish using a 10 μL pipette tip. Photographs were captured immediately after the scratch was made and again after 24 hours. Image J software was used to analyze the migration distance and calculate the migration rate.

Transwell migration assay

The transfected A549 cells were seeded in the upper chamber of a serum-free medium (filter membrane 8 μm). Complete medium containing 10% FBS was added to the lower chamber. After a 24-hour incubation period, the upper chamber was taken out and fixed in a 4% PFA solution at room temperature for 10 min. Subsequently, the cells were stained with a 0.1% crystal violet solution and imaged.

Western blotting

Cells were harvested and lysed in RIPA buffer (Beyotime Biotechnology, Shanghai, China), Subsequently, the protein concentration was measured using the BCA kit (Beyotime Biotechnology, Shanghai, China). An equal amount of cell lysate was separated by SDS-PAGE and transferred onto a PVDF membrane. The PVDF membranes were then blocked with 5% skimmed milk and incubated overnight at 4°C with primary antibodies against WDFY4 (1:1000, Abcam, Waltham, MA, USA). This was followed by incubation with polyclonal HRP-labeled secondary antibodies (1:5000, MultiSciences, Shanghai, China) at room temperature for 1 h. Finally, the membrane was developed using enhanced chemiluminescence.

Xenograft nude mice model

All animal experiments conducted in this study were granted approval by the Institutional Animal Care and Use Committee of Research Center for Drug Safety Evaluation of Hainan Province, Hainan Medical University (2023-0320 (NM01)) and were performed in strict adherence to the established guidelines. BALB/c nude male mice (4–5 weeks old) were provided by the Beijing Vital River Laboratory Animal Technology Co., Ltd. A total of 1 × 10⁷ A549 cells, transfected with control or WDFY4 lentivirus, were subcutaneously injected into the left dorsal region of nude mice. Tumor volume and weight were evaluated at 2-day intervals. Tumor volume was calculated according to the following formula: length × width2/2 (mm3). On the 15th day, the mice were euthanized via cervical dislocation following anesthesia with pentobarbital sodium. Tumor tissue blocks were collected for subsequent immunohistochemical staining analysis.

Immunohistochemical analysis

The IHC assays and IHC scores were performed according to a previously published protocol [20, 21]. High expression levels were defined as a cumulative score of 4 points, while low expression levels were determined by a cumulative score of less than 4 points. The LUAD tissue microarray (HLugA180Su08, XT19-019) was obtained from Shanghai Outdo Biotech, including 98 cases of LUAD tissues and 82 cases of adjacent tissues. The primary antibodies used, including WDFY4, Ki-67, and Caspase-3, were purchased from Abcam (1:200, Abcam, Waltham, MA, USA).

Statistical analysis

The statistical analysis was performed using GraphPad Prism software 8.0 (GraphPad Software, Inc., San Diego, CA, USA). The data were represented as means ± standard deviations. For the cell and animal experimental data, statistical significance was determined using the student’s t-test or one-way ANOVA. The expression level of WDFY4 in lung cancer tissues and noncancer tissues was analyzed by Student’s t-test. Chi-square test was used to analyze clinicopathological features, and Kaplan-Meier method and log-rank test were applied for survival analysis. P < 0.05 was considered statistically significant [22].

Data availability statement

The original contributions presented in this study are included in this article/Supplementary material, the research and analysis data are from the GEO database, and further inquiries can be directed to the corresponding authors.

FPI identification and sample grouping analysis

Each LUAD sample in the TCGA had its FPI quantified based on the GSAV analysis. As shown in Figure 1A, the FPI values in LUAD cancer tissues were lower than those in para-cancerous tissues The FPI values of all TCGA-LUAD cancer samples were then analyzed for survival. Utilizing the survival and survminer algorithms, the optimal cutoff value was determined to establish a correlation between FPI and survival in the samples. This cutoff value was then used to classify the samples into FPI high- and low-expression groups. Figure 1B illustrates the differential genetic profile of the two sample groups. The lower the FPI value, the shorter the long-term survival rate of LUAD patients (Figure 1C).

Figure 1. Construction of an interaction network between the ferroptosis index (FPI) and the cancer immune microenvironment. (A) Variance analysis of standardized FPI between normal and cancer tissues in lung cancer. (B) Changes in the expression profile of ferroptosis-related genes. (C) Survival analysis and post-analysis of standardized FPI. (D) Analysis of immune cell infiltration in lung cancer tissues and relative expression abundance. (E) Correlation heat map of relative infiltration values of local immune cells and FPI. (F) Analysis of strongly associated immune cell subpopulations by Pearson correlation, with FPI as the core of regulation.

Immuno-infiltration analysis and correlation statistics

In order to quantify the relative infiltration of each cell within the LUAD cancer samples, an immuno-infiltration analysis was conducted using ssGSEA, which involved the assessment of 24 immune cell genes as markers (Figure 1D). Through partial correlation analysis, the correlation between FPI values and the standard relative infiltration values of immune cells was determined. The application of automated clustering techniques resulted in the identification of immune cell populations associated with FPI (Figure 1E). The analysis of the interaction network between FPI values and immune cells indicated a significant correlation between the infiltration levels of B cells, CD8 T cells, T helper cells, Th17 cells, NK CD56bright cells, and NK CD56dim cells with ferroptosis (Figure 1F). Among these, B cells exhibited the strongest association with ferroptosis (Figure 1F). Furthermore, the above-mentioned immune cell infiltration was closely related to the long-term survival prognosis (Figure 2A).

Figure 2. Identification of core immune subpopulations and variability analysis of GSVA pathway. (A) Correlation analysis between ferroptosis-related immune cell populations and survival prognosis of lung cancer patients. (B) Variability analysis of the Hallmark pathways in cancer and normal tissues.

GSVA and survival analysis

Subsequently, we conducted GSVA pathway quantification on genes belonging to the HALLMARKER gene set. According to the B-H difference calculation in the LIMMA algorithm, we observed significant upregulation in the expression levels of apical junction, epithelial mesenchymal transition, IL6/JAK/STAT3 signaling, apoptosis and Notch signaling pathways. Conversely, the expression levels of unfolded protein response, G2M checkpoint, oxidative phosphorylation, DNA repair, spermatogenesis and other pathways were significantly downregulated (Figure 2B). Furthermore, we conducted GSVA standard quantification and differential analysis on the immune gene sets. The results indicated that SIG_REGULATION_OF_THE_ACTIN_CYTOSKELETON_BY_RHO_GTPASES, SA_MMP_CYTOKINE_CONNECTION and SIG_BCR_SIGNALING_PATHWAY showed significant up-regulation. Moreover, SIG_REGULATION_OF_THE_ACTIN_CYTOSKELETON_BY_RHO_GTPASES and SIG_BCR_SIGNALING_PATHWAY was strongly associated with the long-term survival of LUAD patients (Figure 3A, 3B). In this case, patients with higher GSVA values of SIG_BCR_SIGNALING_PATHWAY had longer overall survival, whereas the opposite trend was observed for SIG_REGULATION_OF_THE_ACTIN_CYTOSKELETON_BY_RHO_GTPASES (Figure 3B).

Figure 3. Variability analysis of immunoregulation-related pathways. (A) Differential analysis of variability scores of immunoregulation-related pathways in GSEA. (B) Survival analysis of the SIG_BCR_SIGNALING_PATHWAY and SIG_REGULATION_OF_THE_ACTIN_CYTOSKELETON_BY_RHO_GTPASES in patients with LUAD. (C) Enrichment significance analysis of B Cell activation and immunoregulatory pathways.

GSEA analysis

The GSEA flow analysis revealed that the PFI low-score group exhibited elevated expression levels of pathways of B cell activation, B cell proliferation, B cell-mediated immunity and B cell receptor signaling pathway were high in the PFI low-score group, whereas the PFI high-score group displayed lower expression levels of these pathways (Figure 3C). These findings suggest a strong association between ferroptosis and B-cell activation, both of which are closely linked to the long-term survival of patients.

ScRNA cell clustering analysis

Following the exclusion of cells exhibiting abnormal gene numbers and ribosome ratio, cell populations were downscaled by PCA and tSNE. Subsequently, these populations were subjected to automated database annotation (Figure 4A). Notably, T cells emerged as the predominant cell population in LUAD. Intriguingly, Transitional B cells and B cells, despite belonging to the same B cell population, exhibited spatial separation within the cell population space, implying a potential difference in their biological functions. As for the overall cell population, the pseudotime analysis revealed a decrease in the expression levels of LTB, HLA-DRA, HLA-DRB1, CD52, HLA-DPB1 and HLA-DQB1 during the late stage of LUAD progression. In contrast, the expression levels of TXNDC5, CD38, GAS6, FKBP2 and SDC1 were increased in the late stage of LUAD progression, indicating potential distinct roles in the occurrence and progression of the cancer (Figure 4B).

Figure 4. Identification and pseudotime analysis of single-cell subpopulations in lung cancer. (A) Identification of immune cell subpopulations in lung cancer tissues, with T cells as the major subpopulation and B cell populations with two subpopulations spaced apart in the UMAP spatial map, suggesting distinct biological functions. (B) Single-cell pseudotime analysis heat map of B-cell subpopulations. (C–E) Spatial distribution of B-cell populations and Transitional B cells in the pseudotime analysis atlas.

Cell identification of core subpopulations, pseudotime analysis and variance analysis

Based on the above analysis, a subpopulation analysis of Transitional B Cells and B cells was conducted. The pseudotime trajectory revealed the presence of two cell populations with a core regulatory node in cell development (Figure 4C), with Transitional B Cells being classified within the late developmental population (Figure 4D, 4E).

Single cell GSVA and pathway enrichment analysis

We quantified the GSVA pathway in the expression profiles of cells from both populations, suggesting that HALLMARK DNA REPAIR HALLMARK UNFOLDED PROTEIN RESPONSE HALLMARK OXIDATIVE PHOSPHORYLATION and HALLMARK PROTEIN SECRETION were significantly elevated in the B cell population, while HALLMARK KRAS SIGNALING DN, HALLMARK EPITHELIAL MESENCHYMAL TRANSITION and HALLMARK HEDGEHOG SIGNALING demonstrated high expression in the transitional B cell population (Figure 5A).

Figure 5. Enrichment and variability analysis of B-cell subpopulation pathways. (A) Variability scores of B cell subpopulations and Transitional B cells in the Hallmark gene set pathway. (B, C) Results of enrichment significance of B cell subpopulations and Transitional B cells at the REACTOME gene set level. (D) Intersection analysis of pseudotime regulatory genes of B-cell subpopulations and ferroptosis-related differential genes in TCGA.

Results of GSEA revealed a significant association between B cells and the activation of various pathways, including MHC class II antigen presentation, TCR signaling and Generation of second messenger molecules (Figure 5B), Additionally, the activation of Transitional B cells was closely related to processes such as SRP-dependent cotranslational protein targeting to membrane, the stressed Protein Response (UPR) and Cargo concentration in the stressed ER et al. (Figure 5C).

Identification and functional analysis of core regulatory genes

By acquiring intersections, we have identified 14 ferroptosis-related differential genes that are closely associated with B cell activation in LUAD (Figure 5D). Pseudotime analysis has suggested that the expression levels of ACAP1, LINC00926, TLR10, MS4A1, WDFY4 and TRIM22 decrease during at the later stages of B cell differentiation in LUAD, while the expression levels of TMEM59, TP53INP1 and METTL7A increase (Figure 6A). The PPI network highlights the significance of WDFY4 in the biological interaction network and its role as a crucial regulatory target for B cell differentiation in LUAD (Figure 6B). Drug target prediction indicated that Acacetin, Quercetine dihydrate and Diperodon hydrochloride were essential target drugs for WDFY4-mediated B cell differentiation in LUAD (Figure 6C). Furthermore, the spermine catabolic process, polyamine catabolic process and spermine metabolic process appear to significant pathways of WDFY4 (Figure 6D). Additionally, Amino_oxidase, FAD/NAD-binding_dom and TP53INP1 were identified as structural thresholds for target regulation (Figure 6E).

Figure 6. Identification of ferroptosis-related target genes in B cell subpopulations of lung cancer. (A) Developmental changes of common core genes in gene intersection of B cell pseudotime core genes and lung cancer ferroptosis differential genes in the B cell pseudotime heat map. (B) Network of intersection genes at the level of protein-protein interactions. (C) Enrichment analysis of intersection genes in CMAP drug database. (D) Enrichment analysis of intersection genes at the levels of Biological Process, Cellular Component and Molecular Function. (E) Predictive analysis of intersection genes in the InterPro database at levels of substructural interactions.

Correlation between downregulation of WDFY4 and poor prognosis

To conduct a more comprehensive examination of the expression of WDFY4 in LUAD tissues and its relationship with the prognosis of lung cancer patients, we analyzed the expression of WDFY4 in LUAD tissue microarray through immunohistochemistry (IHC). The findings revealed a lower expression level of WDFY4 in LUAD tissues compared to noncancerous tissues (Figure 7A, 7B), and a higher incidence of low WDFY4 expression in LUAD tissues compared to noncancerous tissues (Figure 7C). However, there was no significant correlation observed between the expression of WDFY4 and various clinical pathological features, including age, gender, tumor size. Descriptive grade, lymph node positivity, and TNM stage (Figure 7D, 7E, Table 1). Notably, Kaplan-Meier analysis demonstrated a consistent association between down-regulation of WDFY4 and a poorer prognosis (Figure 7D), indicating the potential oncogenic role and prognostic significance of WDFY4 in LUAD.

Figure 7. The expression and clinical significance of WDFY4 in LUAD. (A, B) Tissue microarray (TMA) analysis by IHC staining showed a low expression of WDFY4 in LUAD tissues. (C) High and low expression rates of WDFY4 in LUAD tissue. (D) Kaplan-Meier analysis showed that LUAD patients with downregulation of WDFY4 were positively correlated with worse prognosis and shorter overall survival. (E) Correlations of WDFY4 expression levels in LUAD tissues and clinicopathological features.

Overexpression of WDFY4 inhibits the proliferation and migration of A549 Cells and promotes apoptosis in vitro and in vivo

To investigate the biological effect of WDFY4 on lung adenocarcinoma cells, we employed the lung adenocarcinoma A549 cell line and used lentivirus to transfect A549 cells to establish a WDFY4 overexpression cell model. Our findings demonstrate that the transfection of ad-WDFY4 significantly increased the expression of WDFY4 in A549 cells (Figure 8A). Moreover, the overexpression of WDFY4 inhibited the proliferation of A549 cells (Figure 8B), as evidenced by a notable increase in the apoptosis rate within the ad-WDFY4 group compared to the ad-Ctrl group (Figure 8C). Furthermore, we employed the Scratch wound healing assay and Transwell migration assay to investigate the effect of WDFY4 on the metastatic potential of the A549 cells. The results depicted in Figure 8D–8F demonstrate that the upregulation of WDFY4 substantially inhibits the migratory capacity of the A549 lung adenocarcinoma cell line. To investigate the potential influence of WDFY4 on the progression of lung cancer in vivo, we subcutaneously injected A549 cells that were stably expressing WDFY4 following lentivirus infection into nude mice. The analysis of subcutaneous tumor formation revealed a significant decrease in tumor volume and weight of nude mice in the ad-WDFY4 group compared to those in the ad-Ctrl group. Furthermore, immunohistochemical analysis demonstrated a lower expression level of the proliferation-related antigen ki-67 in the ad-WDFY4 group compared to the ad-Ctrl group, while the expression level of the apoptosis-related protein Caspase-3 was higher in the ad-Ctrl group (Figure 8G–8I). The data presented herein indicate that the upregulation of WDFY4 exerts a substantial inhibitory effect on the proliferation of lung adenocarcinoma cells both in vitro and in vivo, thereby highlighting its potential as a viable therapeutic target for the treatment of lung adenocarcinoma.

Figure 8. The effect of WDFY4 on the biological characteristics of lung adenocarcinoma A549 cells in vitro and in vivo. (A) The expression level of WDFY4 in A549 lung adenocarcinoma cells transfected with ad-Ctrl and ad-WDFY4. (B) Cell viability was measured by CCK-8 assays. (C) Cells were collected and stained with PI and annexin V-FITC for apoptosis analysis by flow cytometry. (D–F) Scratch wound healing and Transwell migration assay analyzed the metastasis of A549 cells after transfection of ad-WDFY4. (G) A549 cells transfected with ad-Ctrl and ad-WDFY4 were transplanted subcutaneously into nude mice, and the cancer growth curve and cancer weight were monitored (n = 5). (H, I) The expression levels of Ki-67 and Caspase-3 in cancer tissues were analyzed by immunohistochemistry (n = 5).