Expression and clinical data of 304 samples from The Cancer Genome Atlas (TCGA)-CESC were acquired from TCGA Genomic Data Commons portal (https://portal.gdc.cancer.gov/). Metagenomic profiles of the same cohort were sourced from TCGA metagenomic microbiome study (https://ftp.microbio.me/pub/cancer_microbiome_analysis/TCGA). The detailed pipeline for obtaining metagenomic profiles was outlined in a previous study (26). Relative abundance normalization was applied to the metagenomic data, and Streptococcus absence was defined as a relative abundance equal to 0.

DESeq2 was employed to investigate differentially expressed genes (DEGs) between Streptococcus-present and -absent groups. Significance criteria included an adjusted P<0.05 and an absolute log₂-fold change ≥0.5 (27).

Weighted gene coexpression network analysis (WGCNA) is a computational method that organizes genes into clusters or modules based on their coexpression patterns across different samples. This approach reveals the complex relationships between genes, providing insights into their roles in biological processes and their connections to disease phenotypes (28,29). In the present study, WGCNA was applied to identify significant genes associated with Streptococcus in samples from patients with AC. A soft threshold of β=7 and a scale-free topology fitting index (R²) of 0.95 were used for matrix transformation according to the scale-free topological criteria. The modules exhibiting the highest correlation with Streptococcus presence were selected. Key genes related to Streptococcus in the AC samples were then identified if they met the following criteria: Gene significance (GS) >0.5, module membership (MM) >0.8 and P<0.05.

Pathway enrichment analysis was conducted using two web-based tools, including g:Profiler (https://biit.cs.ut.ee/gprofiler/gost; version 2023) (30) and Enrichr (https://maayanlab.cloud/Enrichr/; version June 8, 2023) (31), employing the Kyoto Encyclopedia of Genes and Genomes pathway database (32) and Molecular Signatures Database Hallmark 2020 (33). Statistical significance was established with an adjusted P<0.1 and a minimum of four genes per pathway. Significant enriched pathways were visualized using ggplot2 version 3.4.2 (34).

HeLa cells, also known as human AC cell lines (KCLB number 10002; lot no. 59726; passage no. 98; Korean Cell Line Bank; Korean Cell Line Research Foundation), were thawed, cultured and passaged within T75 cell culture flasks using Minimum Essential Medium (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% fetal bovine serum (Gibco; Thermo Fisher Scientific, Inc.) in a humidified atmosphere with 5% CO₂ at 37°C.

Two subculturing instances were conducted to ensure cellular stability for experimentation. Initially, cells were seeded at 5–6×10⁵ cells per flask. Media exchange occurred every 2 days or as needed based on cellular conditions, with harvesting conducted once cells achieved 90–95% confluence.

GBS or Streptococcus agalactiae strain NCTC 818 [G19] (13813; American Type Culture Collection) was aerobically grown at 37°C in Tryptic Soy Broth (Becton, Dickinson and Company) according to the manufacturer's protocol.

HeLa cells were infected with GBS for 0, 2, 6 and 24 h at a multiplicity of infection equal to 100 (GBS:cell=100:1). Right before the experiment, GBS in the bacterial medium was precipitated and resuspended in human medium (Minimum Essential Medium supplemented with 10% fetal bovine serum), ensuring the concentration remained consistent with the 100:1 ratio. For the control group, human medium without bacteria was administered. Following this, the flasks were incubated for 0 (baseline), 2, 6 and 24 h in a humidified atmosphere with 5% CO₂ at 37°C.

All experiments were conducted in triplicate.

Total RNA was isolated from HeLa cells cocultured with GBS using the RNeasy Mini Kit (74104; Qiagen GmbH) according to the manufacturer's protocol. Ribosomal RNA depletion was performed using the MGIEasy rRNA Depletion Kit (1000005953; MGI Tech Co., Ltd) and library preparation was performed using the MGIEasy RNA Directional Library Prep Set (1000006386; MGI Tech Co., Ltd.). The final library concentration was 113 nM. Sample quantification was performed using the Qubit ssDNA Assay Kit (Q10212; Invitrogen; Thermo Fisher Scientific, Inc.) and the Invitrogen Qubit Fluorometer (Q33216; Thermo Fisher Scientific, Inc.). Sample quality was assessed using an Agilent 2100 Bioanalyzer (G2939AA; Agilent Technologies, Inc.). Sequencing was conducted on the DNBSEQ-G400 sequencer (MGI Tech Co., Ltd.) using the DNBSEQ-G400RS-High throughput sequencing FCL PE100 kit (1000016949; MGI Tech Co., Ltd.), generating 100-bp paired-end reads.

FASTQ files from the sequencer underwent quality assessment using FASTQC v0.12.1 and Multiqc version 1.24.1 (35,36). Cutadapt version 4.9 was employed to eliminate low-quality reads and sequencing adapters (37). The alignment and quantification of RNA-seq data were performed using STAR version 2.7.11b (38). Detailed information on sequenced and mapped samples can be found in Table SI.

TimeSeriesAnalysis (TiSA) is a tool including analysis and visualization packages for RNA-seq and microarrays. It facilitates the extraction of significant genes from time series transcriptomic data by assessing differential gene expression along both condition and temporal axes (39).

In the present study, TiSA was used to identify significant genes meeting the criteria of an absolute log₂-fold change ≥1 and an adjusted P<0.05 from longitudinal transcriptomic data, elucidating the interaction between HeLa cells and GBS.

A principal component analysis (PCA) plot was generated using the first and second principal components to visualize the differences between control and GBS-exposed HeLa cell groups at various time points. This was conducted using the plot_PCA_TS() function from the TiSA package. The partitioning algorithm based on recursive thresholding (PART) from the clusterGenomics R package, version 1.0 (40), was applied for clustering. DEGs and clusters were visualized using heatmaps, created with the ComplexHeatmap package version 2.8.0 (http://www.bioconductor.org/packages/devel/bioc/html/ComplexHeatmap.html). To elucidate the biological meaning of the clusters, Enrichr and g:Profiler were used, with the criteria as aforementioned.

Total RNA was extracted from GBS-exposed HeLa cells after 24 h to validate the mRNA expression of the genes belonging to the MAPK pathway [fibroblast growth factor 21 (FGF21), nerve growth factor (NGF), IL1A and IL1B] and mTORC1 pathway [DNA damage inducible transcript 3 (DDIT3), cystathionine γ-lyase (CTH), asparagine synthetase (glutamine-hydrolyzing) (ASNS) and nuclear protein 1 (NUPR1) using the RNeasy Mini Kit (74104; Qiagen GmbH). RT-qPCR was performed using a one-step TOPreal™ SYBR Green RT-qPCR Kit (cat. no. RT432S; Enzynomics Co., Ltd.) on a Thermal Cycler Dice Real Time System TP800 (Takara Bio, Inc.). This one-step kit combines high-yield TOPscript reverse transcriptase (cat. no. RT002; Enzynomics Co., Ltd.) with chemically modified Taq polymerases. According to the manufacturer's instructions, reverse transcription was performed at 50°C for 30 min, and the thermocycling conditions were as follows: Initial denaturation at 95°C for 10 min, followed by 45 cycles of denaturation at 95°C for 5 sec, and annealing and elongation at 60°C for 30 sec. The primer sequences utilized are presented in Table SII.

All reactions were performed in triplicate. The relative expression levels of the target genes were normalized to the reference gene (GAPDH) using the ΔCq method. Fold changes in gene expression were calculated using the 2^−ΔΔCq method, with the control sample serving as the calibrator (41).

To validate the activation of significantly enriched pathways in HeLa cells after 24 h of exposure to GBS, western blotting was performed. HeLa cells were co-cultured with GBS for 24 h, washed with cold PBS and lysed using cold M-PER™ Mammalian Protein Extraction Reagent containing 0.2% Halt™ Protease Inhibitor Cocktail (Thermo Fisher Scientific, Inc.). Tumor cell lysates (20 µg/lane; protein concentration quantified using a Bradford assay) were loaded onto a 4–20% Mini-PROTEAN TGX precast gel (Bio-Rad Laboratories, Inc.) and electrophoresed for 10 min at 80 V, followed by 60–90 min at 120 V. The proteins were transferred to PVDF membranes for 45 min at 12 V using the semi-dry blotting system (ATTO Corporation). After blocking with 5% skim milk for 1 h at room temperature, the membranes were incubated overnight at 4°C with the following primary antibodies: ERK1/2, phosphorylated ERK1/2 (p-ERK1/2), phosphorylated mTOR at Ser2448 (p-mTOR), p38, phosphorylated p38 (p-p38), panAKT (each at a dilution of 1:1,000), mTOR (dilution, 1:2,000), phosphorylated AKT at Ser473 (p-AKT) and Ki67 (each at a dilution of 1:5,000). Detailed information about the antibodies is provided in Table SIII.

Following incubation, the membranes were treated for 1 h at room temperature with anti-rabbit IgG HRP-linked antibody and anti-mouse IgG HRP for β-actin antibody at a dilution of 1:5,000. The membranes were washed three times with 1X TBS with 0.1% Tween-20 for 10 min. The PVDF membranes were then processed with a PicoEPD Western Blot Detection kit (EBP1073; ELPIS-Biotech, Inc.). The chemiluminescence images were captured using the Amersham Imager 600 (GE Healthcare), and the resulting bands were semi-quantified to determine relative protein levels using ImageJ software version 1.54 (National Institutes of Health). The experiments were conducted in triplicate for both control and treatment groups.

Candidate therapeutic agents were selected through a drug repurposing approach. One valuable tool in drug repurposing is the Connectivity Map (CMap) tool (42). CMap is a powerful platform for unraveling connections among drugs, genes and diseases. The fundamental principle involves comparing the gene expression profiles induced by a specific drug with those associated with particular diseases (42).

Genes from TiSA clusters of interest were uploaded to the CMap web-based tool version 1.1.1.43 (https://clue.io/). Potential anticancer drugs were identified based on their tau score, which is converted from normalized connectivity scores, comparing them with values of that disease for all the drugs in the reference database (43,44). Tau scores exceeding 90 or below-90 between two signatures indicate strong connectivity. A negative connectivity signifies the reversal of the disease signature by the drug (43).

Binimetinib (MEK162) was purchased from Selleck Chemicals (cat. no. S7007). A 10 mM stock solution was generated in DMSO and stored at −20°C.

Ridaforolimus (deforolimus; cat. no. HY-50908; MedChemExpress) was prepared as 10 mM stock solutions with DMSO and stored at −20°C.

Specially tailored working solutions of these drugs were prepared for immediate use. In the case of monotherapy, binimetinib and ridaforolimus were serially diluted to concentrations of 0.625, 1.25, 2.5, 5, 10, 15 and 20 µM. For combination therapy of binimetinib and ridaforolimus, a 1:1 ratio of the drugs was applied, with serial dilutions of 1.25, 2.5, 5, 10, 15 and 20 µM used. For all treatments the cells were incubated at 37°C with 5% CO₂. Each experiment was conducted in triplicate.

A trypan blue assay was utilized to evaluate the viability of HeLa cells following exposure to GBS or therapeutic agents. After washing with PBS, adherent cells in a 90-mm plate were dissociated using 2 ml of 0.25% trypsin-EDTA for 2 min at room temperature. The cells were then resuspended in 8 ml cell culture medium. A portion of this cell suspension was mixed with an equal volume of 0.4% trypan blue (Gibco; Thermo Fisher Scientific, Inc.) for 2 min at room temperature. The number of viable and non-viable cells was counted in triplicate using a hemocytometer.

After reaching confluence in a T75 flask, HeLa cells were subcultured into 96-well plates at a concentration of 2×10⁴ cells/100 µl. Following a 24-h incubation period, binimetinib and/or ridaforolimus were added into each well using serial dilutions of 1.25, 2.5, 5, 10, 15 and 20 µM.

Subsequently, the cells were incubated at 37°C with 5% CO₂ for either 24 or 48 h. Afterwards, 20 µl CellTiter 96^® AQueous One Solution (G3582; Promega Corporation) reagent was added to each well. The plate underwent an additional 4-h incubation at 37°C with 5% CO₂, and the absorbance at 490 nm was measured using a microplate reader (Multiskan™ FC; Thermo Fisher Scientific, Inc.).

The absorbance intensity at 490 nm for both the control and treatment groups was normalized to the absorbance of the background. Subsequently, cell viability percentages for different treatment groups were calculated using the following formula: % viability=adjusted optical density (OD)_sample/mean (adjusted OD_control), where adjusted OD_sample=OD_sample-mean (OD_background), and adjusted OD_control=OD_control-mean (OD_background).

To measure the induction of apoptosis, 5×10⁵ HeLa cells were plated in 90-mm plates and allowed to adhere and proliferate for 48 h. Cells were then treated with either 15 µM binimetinib, 20 µM ridaforolimus, a combination of both (5 µM binimetinib + 5 µM ridaforolimus) or DMSO (vehicle treatment) as a control at 37°C with 5% CO₂.

After 48 h, cells were harvested and transferred into tubes at a concentration of 5×10⁵ cells/ml. Cell death was assessed using an annexin V-FITC apoptosis detection kit (ab14085; Abcam) with 5 µl annexin V-FITC and 5 µl PI, followed by a 5-min incubation at room temperature in the dark. The apoptosis assay data were obtained using a BD FACSCalibur flow cytometer (BD Biosciences). The graph showing the percentages of apoptotic and necrotic cells was created using the Floreada tool (https://floreada.io; WASM version SIMD).

To estimate the efficacy of anticancer drugs, the dose-response curve R package version 3.0.1 was used to determine the IC₅₀ (45). The IC₅₀ values represent the concentrations of each drug required to achieve a 50% inhibition of cancer cell viability.

To predict the combination (binimetinib and ridaforolimus) effect, the synergy score was calculated using the SynergyFinder R package version 3.8.2 (46). If the synergy score is between −10 and 10, the effect of the combination of the two drugs is additive, if the score is >10, the effect is synergistic, and if the score is <-10, the combination is considered antagonistic (46,47).

To evaluate the effect of therapeutic agents on GBS-exposed HeLa cells and to compare these effects with those on HeLa cells not exposed to GBS, GBS and HeLa cells were co-cultured at a ratio of 100:1 in a humidified atmosphere with 5% CO₂ at 37°C for 6 h. Following this exposure, the GBS-exposed HeLa cells were treated with 15 µM binimetinib, 20 µM ridaforolimus, a combination of both (5 µM binimetinib + 5 µM ridaforolimus) or DMSO as a control at 37°C with 5% CO₂. The number of dead and viable cells after 24 h of treatment was measured using a trypan blue assay in triplicate as aforementioned.

To compare two groups, an unpaired Student's t-test was performed. For analyses involving multiple groups, one-way ANOVA followed by post hoc testing using the Holm-Bonferroni correction method was performed (48). All statistical analyses were performed using the stats package version 4.3.3 (https://search.r-project.org/R/refmans/stats/html/00Index.html) in R version 4.3.3 (https://www.R-project.org/). The results were visualized using the ggplot2 package version 3.4.2 (34) and ggpubr package version 0.6.0 (https://github.com/kassambara/ggpubr). Analyses were based on three experimental repeats. In the plots, error bars were included to represent the dispersion of the variable (standard deviation). P<0.05 was considered to indicate a statistically significant difference.

To investigate the relationship between Streptococcus and CC, a cohort of 304 TCGA CC cases, including transcriptomic data and Streptococcus abundance profiles, were analyzed (Table SIV). As shown in Fig. 1A, 31 cases (10.2%) exhibited the presence of Streptococcus.

DEG profiles (Fig. 1B) revealed 128 upregulated and 209 downregulated DEGs identified using DESeq2 (Table SV). Pathway enrichment analysis, conducted using the web-based tools g:Profiler and Enrichr, of the list of 128 upregulated DEGs, highlighted statistically and biologically significant pathways associated with cancer, including ‘MAPK signaling pathway’, ‘mTORC1 signaling’, ‘Ras signaling pathway’ and ‘Wnt signaling pathway’, as well as metabolism pathways such as ‘metabolism of xenobiotics by cytochrome P450’, ‘fatty acid metabolism’, ‘glutathione metabolism’ and ‘arachidonic acid metabolism’ (Fig. 1C; Table SVI).

To further explore the specificity of the relationship between Streptococcus and AC, 31 patients with AC from a cohort of 304 patients with CC in TCGA were analyzed. Among them, 3 patients with AC exhibited the presence of Streptococcus (Table SIV). Using the WGCNA tool, two significant key modules with high correlation values, r=0.77 (P=5×10⁻⁷) for the black module and r=0.59 (P=4×10⁻⁴) for the brown module, were identified (Fig. 2A-C; Tables SVII and SVIII).

A total of 500 key genes related to Streptococcus were selected from the two key modules based on the following criteria: GS >0.5, MM >0.8 and P<0.05. Pathway enrichment analysis showed that these genes were involved multiple oncogenic pathways, including ‘MAPK signaling pathway’, ‘mTOR signaling pathway’, ‘PI3K-AKT signaling pathway’, ‘TGF-beta signaling pathway’ and ‘Ras signaling pathway’ (Fig. 2D; Table SIX).

To further enhance the insights gained from TCGA-CESC data, a time-series experiment was performed to investigate HeLa cell-GBS interactions. Briefly, HeLa cells were infected with GBS for 0, 2, 6 and 24 h, with corresponding control groups treated solely with bacteria-free human medium. Subsequently, the transcriptomic data generated from this experiment were analyzed using the TiSA package.

The time series PCA plot demonstrated a clear separation between the GBS-exposed HeLa cell groups at 6 and 24 h and the other groups. PC1, representing variability over time, and PC2, representing differences between groups, were utilized for this analysis. As shown in Fig. 3, the dominance of PC1, accounting for 45% of the variation, while PC2 accounted for 14%, indicated that the variation over time was more pronounced than the variation between groups, which was especially evident after 6 h. Notably, the groups associated with GBS stimulation at 6 and 24 h exhibited a distinct separation from the other groups.

The conditional differential gene expression analysis for four-time points is presented in Fig. 4. The number of DEGs markedly increased after 6 h of exposure. Specifically, there were 16, 11, 460 and 1,494 DEGs after 0, 2, 6 and 24 h of exposure, respectively.

All DEGs identified from the conditional differential gene expression analysis, totaling 1,981 genes, were utilized for PART clustering, resulting in six clusters, as shown in Fig. 5A and Table SX. Genes in clusters 1, 2, 3 and 4 showed an upward trajectory between 6 and 24 h, while genes in cluster 6 decreased during the same period. Cluster 5 exhibited mixed expression patterns (Fig. 5B).

Functional analysis of each cluster was performed using Enrichr and g:Profiler. Cluster 2 exhibited the most significant pathways, encompassing both oncogenic pathways (‘MAPK signaling pathway’, ‘mTORC1 signaling’ and ‘p53 pathway’) and immune response pathways (‘TNF signaling pathway’, ‘NF-kappa B signaling pathway’, ‘IL-17 signaling pathway’ and ‘NOD-like receptor signaling pathway’) (Table SXI).

The pathway enrichment analyses revealed the upregulation of ‘MAPK signaling pathway’ and ‘mTORC1 signaling’ in the HeLa cells-GBS exposure experiment (Fig. 6) and the TCGA-CESC public data (Fig. 1C). Additionally, activation of the ‘MAPK signaling pathway’ and ‘mTOR signaling pathway’ was observed in the group of 3 patients with AC with Streptococcus presence (Fig. 2D). These findings suggested that the MAPK and mTORC1 pathways may serve significant roles in AC associated with GBS.

To validate the activation of the MAPK and mTORC1 pathways, the upregulated expression of genes associated with the MAPK pathway (FGF21, NGF, IL1A and IL1B) and the mTORC1 pathway (DDIT3, CTH, ASNS and NUPR1) from cluster 2 was first confirmed using RT-qPCR. As shown in Fig. 7C, most of these genes exhibited higher expression levels in HeLa cells exposed to GBS for 24 h compared with control cells, with the exception of NGF. This discrepancy may arise from differences in sensitivity or specificity between RNA-seq and RT-qPCR methods. Overall, the RT-qPCR results support the activation of the MAPK and mTORC1 pathways in response to GBS exposure.

Western blot analysis revealed a significant increase in the ratios of p-ERK1/2 to total (t-)ERK1/2, p-AKT to t-AKT and p-mTOR to t-mTOR in HeLa cells exposed to GBS for 24 h compared with control cells. Although the ratio of p-p38 to t-p38 was higher in the treated HeLa cells than in the control group, the difference was not statistically significant (Fig. 7D). The elevated phosphorylation of ERK1/2 indicated enhanced MAPK pathway activity, while the upregulation of p-AKT and p-mTOR suggested activation of the mTORC1 pathway, further supporting the involvement of these signaling cascades in the cellular response to GBS exposure.

Additionally, the trypan blue assay demonstrated a decrease in viable HeLa cell numbers at 6 h compared with control cells, followed by a significant increase in HeLa cell numbers, reaching levels similar to control cells by 24 h (Fig. 7A). This cell proliferation phenomenon was confirmed by examining the expression of Ki67, a well-known cell proliferation marker. The highest Ki67 expression was observed after 6 h of exposure and the expression was still increased at 24 h compared with that in the control group (Fig. 7B). This pattern may reflect cancer cell proliferation influenced by MAPK and mTORC1 pathway activity after 6 h of exposure to GBS.

Overall, transcriptomic analysis using public data, time-series co-culture experiments, and validation through RT-qPCR and western blotting assays underscored the pivotal role of the MAPK and mTORC1 signaling pathways. These promising results led to the hypothesis that inhibiting these pathways could be effective in treating AC.

To further investigate candidate drugs targeting the MAPK and mTORC1 signaling pathways, the list of genes from cluster 2 was uploaded to the CMap web-based tool. Ridaforolimus (or deforolimus) was chosen as the mTOR inhibitor due to its highest negative range score (Table SXII). Notably, to the best of our knowledge, it has not been extensively studied in the context of AC or HeLa cells.

While CMap suggested several anticancer drugs to inhibit the MAPK pathway, these drugs were previously investigated in HeLa cells (49–51). Therefore, binimetinib, a MEK inhibitor that has shown efficacy in other cancer types (52,53) but has not yet been explored in AC, was selected for further evaluation in the present study.

The MTS assay revealed minimal inhibitory effects of ridaforolimus on HeLa cell proliferation. No significant reduction in cell viability was observed after treatment with ridaforolimus for 24 and 48 h (Figs. 8B and S3). The IC₅₀ of 59.61±10.01 µM further supported this observation (Fig. 8B).

Based on this finding, the highest concentration of ridaforolimus (20 µM) was used for the apoptosis assay. This dose exhibited a weak cytotoxic effect (Figs. 9 and S2C).

A 1:1 ratio of binimetinib and ridaforolimus was administered to assess the impact of combination therapy. SynergyFinder was used to predict potential synergy, and the results suggested that this combination may exhibit an additive effect, as indicated by a Loewe synergy score of ~4.31, with statistical significance at P<0.05 (Fig. 10). Compared with ridaforolimus alone, the combination of binimetinib and ridaforolimus showed significantly greater inhibition of cell viability in the MTS assay. However, significant differences between the combination and binimetinib alone were observed only at concentrations of 10 and 20 µM (Fig. 11).

For the apoptosis assay, 5 µM binimetinib and 5 µM ridaforolimus were used, as this combination effectively reduced cell viability by 50% as indicated by the MTS assay. Although the predictive tool suggested an additive effect, the combination therapy induced notable cytotoxic effects, as evidenced by the highest percentage of annexin V observed 48 h after treatment in the apoptosis assay (Figs. 9 and S2D).

Furthermore, the IC₅₀ of this combination, calculated using the dose-response curve package, was determined to be 6.67±0.45 µM for each drug (Fig. S4). In practical terms, this translates to a dose reduction of 2.16 times for binimetinib and 8.94 times for ridaforolimus to reach the IC₅₀ when combined, emphasizing the potent impact of this combination.

Based on the promising therapeutic effects observed, we hypothesized that these anticancer agents (binimetinib, ridaforolimus and their combination) might affect HeLa cells exposed to GBS more than the group of cells without exposure to GBS. Experiments were conducted to confirm this hypothesis.

After 6 h of exposure of HeLa cells to GBS, the cells were treated with therapeutic agents (monotherapy and combination therapy). As shown in Fig. S5A, both monotherapies and combination therapies exhibited enhanced effects on HeLa cells exposed to GBS compared with the control group after 24 h of treatment, with the combination therapy demonstrating a more pronounced effect than monotherapy.

In addition, the data specifically indicated that combination therapy was more effective in HeLa cells exposed to GBS than in those not exposed to GBS (Fig. S5B).