Sialyltransferase-related genes as predictive factors for therapeutic response and prognosis in cervical cancer

Data collection and pre-processing

The data of somatic mutation and clinical phenotype in CC were obtained and downloaded from The Cancer Genome Atlas (TCGA) database. The samples without survival time information and survival status were removed while retaining those with the survival time >0 day. Relevant RNA-seq expression profile was also downloaded and converted to transcripts per million (TPM) format, followed by the log2 transformation. The final sample of 291 CC patients to be used for subsequent analyses and information on the clinical characteristics of all patients is shown in Table 1. GSE44001 dataset with survival information in CC was downloaded from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE44001) (Lee et al., 2013). The probes were converted into Symbol as per the annotation file, and the samples without overall survival status or clinical follow-up data were also removed. 300 CC samples were then sorted and applied in the analysis. In addition, the 20 SiaTs were all downloaded based on the existing study (Cao et al., 2023).

Weighted gene co-expression analysis (WGCNA)

The specific gene modules were analyzed and identified via “WGCNA” R package (Langfelder & Horvath, 2008), and all processes mentioned were based on a study published elsewhere (Jia et al., 2024). Specifically, the soft threshold for WGCNA was determined with the help of “pickSoftThreshold” function, and then the gene modules were identified via hierarchical clustering based on the criteria of at least 60 genes per module (minModuleSize = 60). “GSVA” R package to be used to calculate SiaTs scores. Next, genes in the module with relatively higher correlation were applied for subsequent analyses (Wang et al., 2024).

Enrichment analysis

The selected genes within the specific module were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses employing “clusterProfiler” R package at the sorting threshold of p < 0.05 (Yu et al., 2012; Xu et al., 2024).

Prognostic model construction and verification

Univariate Cox regression analysis was utilized to determine the association of the selected genes with the prognosis, and then the “glmnet” R package was applied for LASSO Cox regression analysis to narrow down the number of selected genes (Engebretsen & Bohlin, 2019). Multivariate Cox regression analysis was hereafter commenced for the generation of risk coefficient of each identified gene and the following formula was applied for the calculation of RiskScore:

$$\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\mathrm{\Sigma }\mathrm{\beta }\mathrm{i}\times \mathrm{E}\mathrm{x}\mathrm{p}\mathrm{i}$$

(Note: “β” refers to the corresponding regression coefficient of identified genes and “Expi” means the expression of identified genes).

The optimal threshold was then calculated via “survminer” R package to allocate patients into high- or low-risk group (Liu et al., 2021). Thereafter, the Kaplan-Meier (KM) survival analysis with logarithmic rank test was applied to analyze the survival of patients, and the efficacy of prognostic model was evaluated via plotting the receiver operator characteristic (ROC) curves and calculating the corresponding area under the curve (AUC) values.

Tumor immune microenvironment (TIME) analysis

The scores of 10 cell populations were calculated via the “MCPcounter” R package (Zheng et al., 2021). The overall TIME infiltration score was evaluated using the ESTIMATE algorithm. In addition, the relative abundance of six cell types was determined via TIMER online tool (http://cistrome.org/TIMER).

SiaTs mutation landscape

The SiaTs mutation landscape and the gene mutation landscape in different risk groups were plotted based on the datasets processed by mutect2 software from TCGA.

Immuno/chemotherapy response and drug sensitivity analyses

TIDE website (http://tide.dfci.harvard.edu/) was employed to compute the TIDE score and the degree of anti-tumor immune evasion in patients receiving immunotherapy was then predicted. Higher TIDE score was indicative of higher possibility of immune evasion and lower benefit from immunotherapy. Also, relevant immune checkpoint genes were acquired from the HisgAtlas database (https://bio.tools/HisgAtlas) and their levels in different risk groups were compared.

The drug sensitivity in patients of different risk groups in TCGA-CESC was predicted and compared based on the estimated half maximal inhibitory concentration (IC50) value using the “pRRophetic” R package (Geeleher, Cox & Huang, 2014).

Cell culture and transfection

Both human cervical endometrial epithelial cells (HCerEpiC, CP-H058) and CC cell line Hela (CL-0101) commercially obtained from Procell (Wuhan, China) were cultured in minimal essential medium (PM150410; Procell, Wuhan, China) with the supplementation of 10% bovine calf serum (164210; Procell, Wuhan, China) and 1% penicillin-streptomycin (PB180120; Procell, Wuhan, China) at 37 °C with 5% CO₂.

The relevant small interfering RNAs against CXCL3 (target seq: 5′-AGCTTATCAGCGTATCATTGACA-3′) and the control small interfering RNAs were synthesized by GenePharma (Shanghai, China), which were then transfected into Hela cells with the use of Lipofectamine 2000 transfection reagent (11668-030; Thermo Fisher, Waltham, MA, USA) as per the manuals of the producer. All cells were harvested 48 h later for subsequent studies.

Reverse-transcription quantitative PCR

Each sample was run in three independent triplicates for the quantification assay. In detail, total cellular RNA from HCerEpiC and CC cells was extracted by RNAiso Plus total RNA extraction kit (9108; Takara Bio, Shiga, Japan) and subjected to the PCR analysis in the ABI7500 thermocycler (Thermo Fisher, Waltham, MA, USA). In detail, following the extraction of the RNA and the quantification on ND-2000 spectrophotometer (Thermo Fisher, Waltham, MA, USA), the relevant cDNA was prepared from 1 μg total RNA and synthesized with a commercial cDNA synthesis kit (6110A; Takara Bio, Kusatsu, Japan). The PCR was hereafter initiated using the PrimeScript^™ RT reagent kit (RR037Q; Takara Bio, Kusatsu, Japan) in the thermocycler at the following parameters: 95 °C for 30 s × 1 cycle and (95 °C for 5 s and 60 °C for 30 s) × 40 repeated cycles. 2^−ΔΔCT method was adopted to calculate the results, with GAPDH as a housekeeping control (Livak & Schmittgen, 2001; Sindhuja, Amuthalakshmi & Nalini, 2022). See Table 2 for the primers used.

Transwell cell migration/invasion assay

The transfected CC cells were transferred to the upper Transwell chamber (pore: eight-micron meter, 3422; Corning Inc., Corning, NY, USA) added with 200 micron liter serum-free media coated with/out the matrigel matrix (C0371; Beyotime, Shanghai, China). The lower Transwell chamber was added with 700 micron liter complete media and then the cells were cultured for 48 h. Next, all cells migrated/invaded to the lower Transwell chamber were serially fixed and dyed for 20 min, followed by the observation in three randomly picked areas in a light microscope (DP27, Olympus, Tokyo, Japan) (Zhang et al., 2023).

Statistics

The statistical analyses and the data visualization were realized in R software version 3.6.0 (R Core Team, 2019) and GraphPad Prism software version 8.0.2. Difference on the variables between two groups were compared using Wilcoxon test or student’s t test. The correlation was calculated based on Spearman’s correlation test. Logarithmic rank test was applied in the KM survival analysis. Statistically significant difference was defined when the p-value was less than 0.05 (In the figures, ns: p > 0.05; *: p < 0.05, **: p < 0.01, ***: p < 0.001 and ****: p < 0.0001).

Recognition of SiaTs-relevant biomarkers via WGCNA

A total of 20 SiaTs have been identified so far and their abnormal expression pattern in CC was observed. Specifically, it was noted in 34 CC sample that 14 SiaTs were mutated, with ST6GALNAC3, ST6GAL1 and ST8SIA5 ranking among the first three (Fig. 1A). The SiaTs score was hereafter calculated and relevant SiaTs-relevant gene modules were recognized and established. The samples were then clustered and screened for co-expression modules. The soft threshold β was set to 6 so as to ensure the scale-free network (Fig. 1B). Next, the specific gene modules were recognized via hierarchical clustering, and 17 modules were generated following the merging process (Fig. 1C). Among these modules, the grey module was the one which failed to be aggregated into other modules. The number of genes in the modules was listed in Fig. 1D. The correlation between the gene module and SiaTs score was determined in the end, where a strong correlation was seen in between the black module and the SiaTs score (Fig. 1E).

Enrichment analysis on genes in the black module

The ClusterProfiler R package was hereafter applied for the enrichment analysis on the genes in black module at the sorting threshold of p-value < 0.05. The relevant results were seen in Figs. 2A–2D. Based on the KEGG analysis, the genes were mainly enriched in metabolism-relevant pathways like Glycolysis/Gluconeogenesis, Glycosphingolipid biosynthesis-lacto and neolacto series and Alanine, Amino sugar and nucleotide sugar metabolism, aspartate and glutamate metabolism (Fig. 2A). The enrichment analysis on the relevant biological process has hinted the significant enrichment on organic anion transport, carboxylic acid catabolic process and O-glycan processing, small molecule catabolic process (Fig. 2B). Also, genes in the black module were observed in certain cellular components like ciliary part, apical part of cell as well as apical plasma membrane (Fig. 2C). In addition, the enrichment analysis on the molecular function has demonstrated the enrichment of these genes in oxidoreductase activity, acting on the CH-OH group of donors, NAD or NADP as acceptor, aspartic-type endopeptidase activity and aspartic-type peptidase activity (Fig. 2D).

Development of a prognostic model and verification

The samples in TCGA-CESC were grouped at the ratio of 7:3 into the train cohort and the test cohort. Then the genes in the black module were sorted at the parameters of p-value < 0.05 and correlation coefficient with SiaTs score > 0.3. A total of 445 genes were thereafter obtained and further narrowed down via Cox regression analysis based on the data in the training set at the threshold of p-value < 0.05, followed by the LASSO Cox regression analysis (Fig. 3A). Finally, eight genes were identified and applied for developing a prognostic model based on the following formula (Fig. 3B):

$$\begin{array}{rl}\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=& 0.455\ast KCNK\mathit{15}+0.26\ast CXCL\mathit{3}-0.715\ast PIH\mathit{1}D\mathit{2}\\ & -0.524\ast DTX\mathit{1}+0.295\ast LIF-0.443\ast TCN\mathit{2}-0.26\ast SERPINF\mathit{2}\\ & -0.286\ast GCNT\mathit{2}\end{array}$$

Patients were then allocated to different risk groups by the optimal cut-off value and the good prognosis was seen in patients of the low-risk group in the training and validation cohorts as well as the TCGA cohort (Figs. 3C–3E). The AUC value was calculated to reflect the prediction accuracy. High AUC value was seen in the training and validation cohorts as well as the TCGA cohort (Figs. 3F–3H). It was also noticed that the percentage of mortality in patients of high-risk group was evidently higher, based on the Chi-square test (Fig. 3I).

To better characterize the model robustness, similar validation was applied in the dataset GSE44001. Notably, the low-risk group had a better survival (Fig. 3J) and a high AUC value (Fig. 3K).

Immune cell infiltration in different risk groups

The immune cell infiltration was determined with ESTIMATE algorithm to compare the difference on the TIME in patients of different risk groups. A higher immune infiltration was seen in low-risk patients (Fig. 4A). With the purpose of exploring the correlation between the RiskScore and immune functions, the score of 6 types of immune cells was calculated by TIMER, discovering the higher score of B_cell, CD4_Tcells and Macrophages in the low-risk group (Fig. 4B). This means that low-risk patients had a stronger anti-tumor immune response. Furthermore, the association between RiskScore and the immune cell score was also calculated and determined by MCP-counter method. Visibly, the RiskScore was negatively correlated with NK cells, T cells, cytotoxic lymphocytes, monocytic lineage, CD8 T cells, myeloid dendritic cells, and B lineage (Fig. 4C), hinting the higher RiskScore was indicative of a lower immune infiltration status. In addition, the correlation in RiskScore and prognosis- and immune checkpoint-related genes were evaluated, and it was noted that RiskScore was negatively linked to the following genes including HAVCR2, CD28, LGALS9, PDCD1 and CTLA4 (Fig. 4D). In conclusion, these results indicate that patients in the low-risk group had higher levels of immune cell infiltration and stronger immune responses, suggesting that their immune microenvironment may be more active.

Abnormality in gene mutation and relevant pathways in different risk groups

The gene mutation status in different risk groups was investigated, and the top 10 genes with relatively higher mutation frequency in these two groups were listed in Figs. 5A,5B. It was seen that TTN and PIK3CA ranked the first and second in both risk groups, while KMT2C and MUC16 occupied the third position in the high- and low-risk groups, respectively.

The ssGSEA was then initiated using TCGA dataset to determine the relevant pathways in different risk groups. The ssGSEA score of each sample in HALLMARK Pathway was calculated and the correlation between the RiskScore and the calculated score was additionally analyzed. It was noted that the RiskScore was positively correlated with MYC Targets V2 and TGF Beta Signaling yet negatively correlated with Interferon Gamma Response (Fig. 5C), thus hinting that higher RiskScore may be associated with the promoted growth of tumor cells and the suppressed immunoreactivity.

Next, the differential activated pathways in different risk groups were compared and analyzed. The dataset h.all.v2023.1.Hs.entrez.gmt was applied for GSEA and the significant enrichment was defined by the threshold FDR < 0.05. Relevant results have demonstrated that Angiogenesis and MYC target V2 may be the differentially activated pathways in high-risk patients, while those in low-risk patients were KRAS Signaling DN and Interferon alpha response (Figs. 5D,5E).

Difference on the immuno/chemotherapy response in different risk groups

First, we compared the different immunotherapy response and applied TIDE software to assess the potential efficacy of immunotherapy in different risk groups. It was uncovered in Fig. 6A that insignificance on the TIDE score was seen in the two groups, whilst the lower Exclusion score and higher MDSC score were observed in the high-risk group than those in the low-risk group. Meanwhile, we also determined the difference on the expression levels of immune checkpoint genes in these groups, and it was unveiled that most genes had a high expression in the low-risk group, including CTLA4, HAVCR2, LAG3, TIGIT, and PDCD1 (Fig. 6B). Furthermore, the response to the chemotherapy was analyzed in different risk groups of TCGA-CESC cohort based on the estimated IC50 value, and the lower IC50 values of paclitaxel, cisplatin, crizotinib, erlotinib, sunitinib, and sorafenib were seen in high-risk group in comparison with those in low-risk group (Fig. 6C). These results indicate that low-risk patients with CC have higher levels of immune infiltration and more active levels of immune checkpoint expression, suggesting a starting point more amenable to immunotherapy, while high-risk patients are more likely to benefit from chemotherapy due to the presence of more immunosuppressive factors.

Quantification on the SiaTs-related genes in CC cells

For subsequent cellular validation, we firstly measured the expressions of all 8 SiaTs-related genes in CC cells Hela and HCerEpic cells. Relevant results have revealed the differential expressions of these genes in CC cells Hela, where the higher expression levels of KCNK15, LIF, TCN2, SERPINF2, and CXCL3 and the lower expression levels of PIH1D2, DTX1 and GCNT2 were observed in comparison with those in HCerEpic cells (Fig. 7A, p < 0.05).

CXCL3 has been reported to play a key role in tumor development and progression (Cheng et al., 2023). To further investigate its potential function in CC, we designed specific small interfering RNAs and transfected them into CC cells as needed. The in-vitro migration and invasion experiments of the effects of CXCL3 silencing on the CC cells showed that CXCL3 silencing could visibly reduce the number of migrated and invaded Hela cells at 48 h (Figs. 7B, 7C, p < 0.001).