Integrative Analysis of Shared Pathogenic Genes and Potential Mechanisms in Gardnerella vaginalis and Persistent HPV16 Infection

Abstract

Bacterial vaginosis, often accompanied by Gardnerella vaginalis (GV) overgrowth, is associated with persistent high-risk human papillomavirus (HR-HPV) infection, particularly HPV16. This study integrated transcriptomic data from in vitro GV infection experiments and a GEO dataset (GSE75132) of HPV16 persistence to elucidate shared pathogenic mechanisms. Differential expression analysis identified 4115 genes modulated by GV infection and 861 by HPV16 persistence, with 74 common differentially expressed genes (DEGs) displaying consistent trends. Enrichment analyses revealed that these DEGs participate in metabolic pathways, immune defense, host-pathogen interactions, and carcinogenesis. Protein-protein interaction networks and Random Forest (RF) feature selection pinpointed radical S-adenosyl methionine domain containing 2 (RSAD2) and Interferon-induced protein with tetratricopeptide repeats 1 (IFIT1) as central hub genes. Upstream transcription analysis identified the homer_AGTTTCAGTTTC_ISRE motif and established a ceRNA network involving hsa-miR-654-5p, IFIT1/RSAD2, and lncRNAs. Mendelian randomization (MR) and colocalization analyses linked RSAD2 downregulation to an increased risk of cervical carcinoma in situ (rs2595163, PPH4 = 0.62), while ROC analysis demonstrated strong diagnostic potential for the combined hub gene expression. Notably, single-cell transcriptomics revealed distinct RSAD2 and IFIT1 expression patterns in immune and epithelial cells during the progression from HPV infection to cervical cancer. Collectively, these findings support RSAD2 and IFIT1 as promising biomarkers and therapeutic targets for HPV-related cervical lesions.

Keywords: Gardnerella vaginalis; IFIT1; RSAD2; molecular mechanisms; persistent HPV16 infection.

Figure 1

DEGs co-regulated by GV infection and persistent HPV16 infection, and their enrichment analysis. (a) GV co-culture model diagram; (b) Grouping information for persistent HPV16 infected populations; (c) Volcano plot of DEGs regulated by GV infection; (d) Volcano plot of DEGs regulated by persistent HPV16 infection; (e) Venn diagram of common DEGs (co-regulated by GV infection and persistent HPV16 infection), showing the commonly upregulated genes; (f) Venn diagram of common DEGs, showing the commonly downregulated genes; (g) Bubble chart of GO enrichment analysis for common DEGs. The top 10 most significant terms (sorted by p-value) are shown; (h) Treemap of KEGG enrichment results for the common DEGs. DEGs, differentially expressed genes; GV, Gardnerella vaginalis; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 2

Research process flowchart. CIS, Carcinoma in situ; DEGs, differentially expressed genes; GO, gene ontology; GSEA, gene set enrichment analysis; GV, GV, Gardnerella vaginalis; HC, healthy control; HSIL, high grade squamous intraepithelial lesion; KEGG, Kyoto Encyclopedia of Genes and Genomes; LISL, low grade squamous intraepithelial lesion; MR, Mendelian randomization; SCC, Squamous cell carcinoma; scRNA-seq, single-cell RNA sequencing data; TCGA-CESC, data sourced from The Cancer Genome Atlas-Cervical squamous cell carcinoma and endocervical adenocarcinoma.

Figure 3

Identification of hub genes co-regulated by GV and persistent HPV16 infection (a) PPI network of the common DEGs; (b) Upset Venn diagram of the top 10 genes ranked by nine Cytohubba algorithms. The red bars represent genes consistently ranked within the top 10 across all nine algorithms (key genes); (c) Feature genes associated with GV infection; (d) Feature genes associated with persistent HPV16 infection; (e) AUC curve of the Random Forest model for predicting GV infection; (f) AUC curve of the Random Forest model for predicting persistent HPV16 infection; (g) Venn diagram of common feature genes shared by GV and persistent HPV16 infection; (h) Venn diagram showing the intersection between the common feature genes identified by random forest for GV and persistent HPV16 infections and the key genes ranked in the top 10 by all nine Cytohubba algorithms in the PPI network. AUC, area under the curve; DEGs, differentially expressed genes; GV, Gardnerella vaginalis; PPI, protein–protein interaction.

Figure 4

Upstream transcriptional regulation analysis of hub genes. (a) Displays the top three motifs by highest AUC values. Red line shows each motif's average recovery curve, green line is the average plus standard deviation, and blue line represents the recovery curve for the current motif. The maximum distance between the current motif and green line marks the peak enrichment level; (b) Venn diagram of miRNAs targeting hub genes from the miRwalk database and HPV-related miRNAs from the HMDD database; (c) ceRNA network of hub genes, with orange representing miRNAs, green representing mRNAs, and purple representing lncRNAs. AUC, area under the curve; HPV, human papillomaviruses.

Figure 5

Single-gene GSEA enrichment analysis of hub genes (GO gene sets). GSEA enrichment analysis of IFIT1 in the GV infection dataset (a) and the persistent HPV16 infection dataset (b), with five GO terms co-enriched by IFIT1 and RSAD2 across both datasets highlighted in the plots; GSEA enrichment analysis of RSAD2 in the GV infection dataset (c) and the persistent HPV16 infection dataset (d), with the same five GO terms highlighted; (e) Venn diagram of GO terms commonly enriched by both IFIT1 and RSAD2 across the two datasets; (f) Detailed list of the remaining 12 GO terms consistently enriched by both genes in both datasets. All GO terms shown in the plots are upregulated entries. GO, gene ontology; GV, Gardnerella vaginalis; GSEA, gene set enrichment analysis.

Figure 6

Diagnostic potential of hub genes in samples with HPV16 and GV infection. (a–d) ROC curve of the two shared genes in test (a, b) and validation datasets (c, d) for HPV16 and GV infection; (e–h) ROC curve of the multi-marker diagnostic model in test and validation datasets for HPV16 and GV infection; (b, f) GV infection in Siha; (d, h) GV infection in ECT1/E6E7. GV, Gardnerella vaginalis; ROC, receiver operating characteristic.

Figure 7

Analysis of the relationship between hub genes and cervical cancer onset and prognosis. (a) Scatter plot illustrates the association between the effects of SNPs on RSAD2 and their effects on CIS; (b) Leave-one-out sensitivity analysis for RSAD2 on CIS; (c) Pleiotropic associations between RSAD2 and CIS. Top plot: shows −log₁₀ (p values) of SNPs from the CIS GWAS, with solid red rhombi indicating probes that passed the HEIDI test, middle plot: displays eQTL results, bottom plot: illustrates the locations of genes associated with the probes; (d) SMR indicating significant negative causal relationships between RSAD2 expressions and CIS onset (p_SMR < 0.05, p_HEIDI > 0.05). (e) Manhattan plot for the colocalization analysis of RSAD2 and CIS, where the diamond-shaped purple points represent the causal SNPs. (f) Comparison plot of loci for RSAD2 and CIS, with the diamond-shaped purple points indicating the lead SNP with the smallest p value. (g) Survival analysis of high-expression versus low-expression groups of RSAD2 in cervical cancer patients; (h) Survival analysis of high-expression versus low-expression groups of IFIT1 in cervical cancer patients. CIS, carcinoma in situ; RSAD2, radical S-adenosyl methionine domain containing 2; SNPs, single nucleotide polymorphisms.

Figure 8

Single-cell transcriptomic characterization of cervical tissues from HPV16-infected patients and cervical cancer patients. (a) UMAP visualization of 14,253 cells from the HPV16 infection scRNA-seq dataset (GSE171894); (b) UMAP visualization of 22,234 cells from the cervical cancer scRNA-seq dataset (GSE168652); (c, d) UMAP visualization of cell clusters corresponding to different sample groups in HPV16-infected and cervical cancer patients; (e, f) UMAP visualization depicting the expression levels of two hub genes in HPV16-infected and cervical cancer patients; (g, h) Dot plots illustrating the expression of the two hub genes across various cell populations in HPV16-infected and cervical cancer patients; (i, j) Violin plots comparing the expression of the two hub genes in HPV16-positive vs. HPV-negative cases, as well as in cervical cancer versus normal tissues. UMAP, Uniform Manifold Approximation and Projection.