A systematic literature search of PubMed, Web of Science, and Embase was carried out for articles published between January 01, 2000, and February 28, 2022. The keywords included HPV [all fields], HPV infection [all fields], human papillomavirus [all fields], cervical intraepithelial neoplasia (CIN) [all fields], cervical lesion [all fields], cervical dysplasia [all fields], cervical neoplasm [all fields], cervical cancer [all fields], cervical carcinoma [all fields], microbiome [all fields], microbiota [all fields], and flora [all fields].

The inclusion criteria for the study, were as follows (1): all observational studies (including case-control, cross-sectional, prospective, and retrospective cohort studies or interventional studies) (2); women with HPV infection (HPV+), HPV persistence, low-grade squamous intraepithelial lesions (LSIL), HSIL, and Ca – categorized according to the results of HPV testing and histology of cervical biopsy, compared to healthy controls (HPV-negative and cytology-negative women) or HPV-negative women (HPV-), cervical cytology was used for classification if histology was not reported (3); studies reporting microbiota analyses of vaginal samples based on high-throughput sequencing (4); studies with information on vaginal microbial alterations at the community level (alpha and beta diversity) and/or taxonomic levels (relative abundance of different microbial species).

Records were excluded if (1) they did not use high-throughput sequencing approaches (2), sequenced the microbiota from cervical or cervicovaginal swabs [because of anatomy-associated discrepancies in the microbiome (21)] (3), analyzed the microbiome in intestinal samples (4), were duplicates (5), were not published in English, or were review articles or conference abstracts.

Prospective or interventional studies reporting altered microbial composition without relevant baseline measurements were also excluded.

Study characteristics and microbiome quantification methods were extracted from the included references. Study characteristics include the following details: publication information, subject demographics, and clinical features, whereas the latter consists of a sequencing platform, sequencing target, DNA extraction, data analysis pipelines, and reference databases. As major outcomes of the analysis, we extracted microbial alterations at the community level (alpha and beta diversity) as well as taxonomic levels (relative abundance of different microbial species). Within- and across-disease findings for the relative abundance of taxonomic assignments were qualitatively synthesized and categorized as increased, decreased, or inconsistent abundancies. No findings were consistent, with less than 75% agreement between studies reporting this taxon (22).

We utilized the National Institutes of Health (NIH) National Health, Lung, and Blood Institute Study Quality Assessment Tool for Observational and Cross-sectional Studies to evaluate the internal validity and potential bias of the included studies. Study quality was rated as “good”, “fair, or poor” by two authors (MW and HL), with differences addressed through discussion.

A total of 3322 relevant records were obtained in the initial electronic search, including 1042 from PubMed, 912 from Web of Science, and 1368 from Embase. Six hundred ninety-nine studies were imported after removing duplicates and screening the titles and/or abstracts. Of those, studies were excluded due to non-vaginal samples (n=15), non-high-throughput sequencing approaches (n=2), no information on any change in taxa compared with HPV- (n=1), and longitudinal studies in the absence of controls (n=2) after assessing the full text (n=42). Ultimately, 22 studies were selected for the final analysis. A flowchart of this process is shown in Figure 1 .

The 22 studies consisted of 17 case-control, three cross-sectional, and two longitudinal studies. Most studies (72.7%, 16/22) were conducted in Asia (China, Korea), five (22.7%) were performed in Western countries (USA, UK, Mexico, and Sweden), and one (4.5%) was performed in Africa (Nigeria). In the 17 case-control comparisons, six studies compared HPV+ or HPV persistence to the controls, ten studies compared women with cervical lesions and/or Ca to the controls, and 3 included women with HPV infection. Only one study compared ≥ HSIL with LSIL. Of these, one study grouped samples into LSIL, HSIL, and normal cytology groups based on observed cytology. In the two prospective studies, comparisons were performed between cases before treatment (local excisional treatment and neoadjuvant chemotherapy) and untreated controls. The detailed information is provided in Table 1 .

Control samples were defined as healthy women in most studies, whereas in two studies the controls were defined as having normal cytology. The controls were HPV positive in one study.

The different microbiome analysis methods are shown in Table 2 . Sequencing approaches included 16S rRNA (95.5%, 21/22) and metagenomics (4.5%, 1/22). Among the studies employing 16S sequencing, sequencing was predominantly performed using the Illumina MiSeq platform (52.4%, 11/21), followed by the HiSeq platform (33.3%, 7/21), the Ion Torrent PGM platform (9.5%, 2/21), and the Novaseq platform (4.8%, 1/21). Five different hypervariable regions were amplified, including V1-V2 (4.8%, 1/21), V1-V3 (4.8%, 1/21), V3 (9.5%, 2/21), V4 (42.9%, 9/21), and V3-V4 (33.3%, 7/21). However, one study did not provide a sequencing target. The bioinformatics pipelines used included QIIME (52.4%, 11/21), Mothur (19.0%, 4/21), QIIME and DATA2 (4.8%, 1/21), DATA2 (4.8%, 1/21), and USEARCH (4.8%, 1/21). The widely used reference databases for the taxonomic assignment were the SILVA (33.3%, 7/21), Greengenes (23.8%, 5/21), and Ribosomal Database Project (14.3%, 3/21) databases.

A total of 57 alpha-diversity assessments were performed across 21 studies ( Figure 2 ) Reported metrics included measurements of richness (observed species/operational taxonomic units, Chao1, abundance-based coverage estimator; 43.9%, 25/57), richness and evenness (Shannon, Simpson, inverse Simpson; 51.7%, 30/58), biodiversity (Faith’s phylogenetic diversity, 1.8%, 1/57) and sequencing depth (Good’s coverage, 1.8%, 1/57). Most of the reported metrics were based on the Shannon index (36.8%, 21/57).

When comparing HPV+ vs. controls or HPV−, five studies examined community richness using different metrics. Two studies indicated that there were no significant differences, two found significantly higher richness indices, whereas one found a non-significant decrease in richness indices. Regarding diversity, combining richness and evenness, significantly higher diversity was observed in three studies, significantly lower diversity was observed in one study, and no significant differences were observed in three studies. One study found a significant increase in Faith’s phylogenetic diversity.

When comparing HPV persistence vs. controls, lower richness indices were reported in two studies, and lower alpha diversity based on a combination of richness and evenness data was obtained in one study.

In comparisons between CIN and controls, two studies investigated higher richness and diversity by incorporating richness evenness indices in LSIL, four studies reported increases in both metrics in HSIL, and four studies reported a similar pattern in CIN. Among women with Ca, 4 studies reported differential richness data. Higher richness was reported in 2 studies, and lower richness was reported in the remaining studies. Differential diversity indices combining richness and evenness were reported in five studies, four of which reported higher diversity and two of which reported lower diversity. Four studies examined diversity based on a comparison between HSIL and LSIL, two of which reported no significant observations in terms of richness and/or evenness data, with the remaining studies providing higher indices. In addition, a significantly lower Good’s coverage was reported in one study comparing HSIL and controls.

Beta diversity comparisons were evaluated in 15 studies, with 13 studies comparing subgroup vs. controls, one study comparing combined subgroups vs. controls, and one study comparing within-subgroup results. Six of 13 studies identified significant differences between different samples, with 4 of 13 showing visually separated components with no significance, 2 finding no clear separation between groups, and 1 study without a report. In the across-group (controls excluded) comparisons, two studies reported no significant differences between subgroups (≥ HSIL vs. ≤ LSIL, HSIL vs. LSIL, and CIN vs. HPV+). One subgroup (HPV+ + LSIL + Ca) was significantly different from the controls. This study also showed a significant difference between the combined subgroup (LSIL + HSIL) and HPV infections. The beta diversity assessments are presented in Table 3 .

Twenty-one studies have reported the relative abundance of vaginal microorganisms at the phylum, family, genus, and species levels. Between-subgroup comparisons of CIN were conducted relative to controls or in the order of abundance in subgroups reported by previous studies.

Figure 2 presents a summary of within- and across-disease comparisons for HPV-related diseases at the phylum, genus, and species levels. Taxa reported as increased or decreased in the subgroups of each study were combined. High inconsistencies were observed at different taxonomic levels in both comparisons. The following results were obtained for the taxa provided by at least two studies:

The differential numbers of some taxa were consistent between the subgroups. For example, higher Bacteroidetes was common to both HPV+ and HSIL; lower Firmicutes was common to HPV+, LSIL, and HSIL; higher Streptococcus was common to HPV+, HSIL, and Ca; lower Lactobacillus was common to all subgroups. The genera Gardnerella and Prevotella were inconsistently present in HPV+, HSIL, and Ca; Sneathia was inconsistent in HPV+ and LSIL; L. iners was inconsistent in HPV+ and HSIL ( Figure 3 ).

Most studies were rated ‘fair’ (63.7%, 14/22). Only four studies were rated as ‘Good’ (18.2%, 4/22), which not only included study features, microbiome methods, and rich reporting of major outcomes but also considered several potential confounders in the reports. Four studies were rated as ‘poor’ quality (18.2%, 4/22). The major sources of potential bias were the lack of reporting of tools or methods for measuring outcomes, lack of detailed information on demographic or clinical characteristics, and lack of detailed definition of cases. The internal validity of the included studies is shown in Table 4 .

Overall, our synthesis indicated a significant increase in richness and diversity in women with cervical lesions and cancer, but no difference was observed between HSIL and LSIL. Medium and large studies indicated a non-significant relationship between richness and/or evenness indices in women with HPV infection, suggesting that alpha diversity indices are of limited utility as a measurement of vaginal health, or for distinguishing between HPV infection and controls. A large number of intrinsic and extrinsic factors account for variations in diversity at the baseline, including race, age, menstrual cycles, sexual activities, contraception use, smoking, and diet, thereby complicating the measurement (33, 41–46).

However, reports with inconsistent data have been published. A few studies have reported a reduction in the richness and diversity of HPV infections. A similar pattern was observed in women with HPV persistence and Ca. These residual heterogeneities for certain diseases may be attributed to methodological effects and the study design. Alpha diversity is particularly affected by measurement errors in microbiome research, owing to the heavy bias of commonly used estimators (47). Although the number of sequences affects the alpha diversity estimates, several confounders in combination with data pre-processing may influence the data output (47, 48). Furthermore, richness metrics do not account for unobserved taxa or provide variance estimates (47). These multifactorial aspects may have contributed to the heterogeneous findings. It is worth noting that in one study investigating lower alpha diversity in Ca, a higher proportion of HPV+ women were a part of the control cohort, which may have resulted in less harmonization in the synthesis.

Regarding beta diversity, most studies have observed a significant cluster between the samples. In summary, two of the studies reported no significant differences between HSIL and LSIL, which is consistent with the pattern observed for alpha diversity.

Among the multiple demographic and clinical features that may account for the extensive divergences across studies, existing evidence indicates that regional and control selection are dominant. The geographical region and the correlative effect of ethnicity can influence the configuration of the vaginal microbiome (60). Our findings indicate that some of the disturbances may be specific to Asian cohorts (i.e., a differential abundance of L. iners). At the same time, divergence was observed in the categorization of subgroups and controls, which is problematic in disentangling the panel of vaginal microbial disturbances observed in HPV infections and cervical neoplasms. For example, Laniewski et al. (24) reported a significant increase in L. crispatus in HSIL versus HPV infection. However, a significant increase in L. iners was observed in another group compared to that in the control group (34).

Despite advancements in bioinformatics techniques in this field, vaginal microbiota studies continue to face methodological challenges. Another source of bias is methodological variations across a range of laboratory processing methods (i.e., DNA extraction, sequencing target, and platform) and data preprocessing methods (i.e., reference database and quality filtration criteria for sequences). Our analyses may also suffer from the use of different sequencing platforms as well as differences in hypervariable regions since heterogeneity in evaluating microbial diversity has been reported (60). For example, five studies reported enriched Prevotella in women with HPV infection using 16S rRNA gene sequencing, although a decrease was observed in one study using 16S rRNA amplicon pyrosequencing. This observation was also similar to the inconsistent changes in the taxa Sneathia reported by 1 study. Technical and clinical factors across studies may increase the difficulty of comparing effect sizes, indicating the need for consistency among methodologies and encouraging data sharing with sufficient metadata. The VAginaL community state typE Nearest CentroId Assifier (VALENCIA) clustering tool (12) is a novel established resolution classification program for the assignment of vaginal microbiome patterns that may be a promising step toward methodologically reinforcing consistency and reproducibility.

Other covariates, including hormonal fluctuations, contraception use, temporal activity, and HPV subtypes, are also intimately correlated with perturbations in the vaginal microbiome (23, 33, 41, 61–63). There was evidence of an association between HPV 16/18 infection and specifically increased genera, viz., Gardnerella, Prevotella, and Atopobium (62); however, insufficient evidence was presented in our analysis. Thus, a rigorous collection of these factors and their careful consideration in assessments and interpretation should be employed by future research groups.

To our knowledge, this is the first systematic review analyzing vaginal microbiome disturbances across a spectrum of HPV-associated diseases globally, thus forming the basis for the identification of reproducible and generalizable potential biomarkers. This enabled us to further explore its clinical use in predicting the severity of HPV-related cervical diseases as a non-invasive diagnostic tool. Functional analyses have provided new insights into the roles of specific bacteria in cervical carcinogenesis. In particular, examining the role of Lactobacillus species may provide better rational targets for the advancement of novel probiotic-based prevention and treatment agents. It has been shown that maintenance of the vaginal microbiota and improvement of cervical epithelialization favors regression of cervical lesions through a prebiotic effect (64). However, the present study has some limitations. First, none of the selected studies provided microbial information for cervical or cervicovaginal samples, which prevented adequate comparison of the relationship between HPV-induced cervical disease and vaginal microbiota. The decision to remove records in these samples was dictated by knowledge of the anatomic potential and substantial site-associated discrepancies in the microbiome (21, 62). Second, the geographic distribution of the included studies indicates a high proportion of studies performed in China. This circles back to the link between regional variances and vaginal bacterial clusters, an imbalance that may have affected data synthesis in our analysis. Third, most studies had relatively small or moderate sample sizes, indicating that our results may still be underpowered. Fourth, given that the selected studies varied widely in terms of study design and bioinformatic analyses, meta-analysis of the selected studies seem complicated. Therefore, we performed qualitative systematic analyses rather than a meta-analysis. Further meta-analysis based on the harmonization of methodologies may provide robust and reproducible taxonomic changes. Finally, this systematic review aimed to integrate current analyses that commonly use 16S rRNA amplicon sequencing at the horizontal taxonomic level, rather than functional profiling. The generated evidence implies that local metabolic patterns correlated with BV include amino acid, dipeptide, polyamine, and ketone body pathways (55). Recent integrative work highlighted the close association of 3-hydroxybutyrate, macrophage migration inhibitory factors, pathobionts, and dysbiotic microbiota with cervical carcinogenesis (65). Relatively little is known about the sophisticated interactive mechanisms among vaginal microbiota, metabolites, and the host. Given the recognized functional redundancy (11), harnessing multi-omics techniques will provide a window on functional candidates or metabolites to clarify the contribution of host-microbe interplay in cervical tumorigenesis.

Going forward, definitive causality for the role of specific taxa in HPV clearance and cervical diseases will stem from a more molecular-related epidemiological study at a larger scale, as well as identification of the mechanisms involved in both clinical research and experimental animal research. It is important to consider individual- and disease-related covariates using multi-omics techniques. Moreover, elucidating the potential for therapeutic manipulations of the vaginal microbiota holds promise for improving outcomes in HPV clearance and cervical lesions.