All patients who visited our hospital between October 2019 and November 2020 with RV were enrolled in the Department of Obstetrics and Gynecology Hospital. After obtaining informed consent, participants were screened to confirm their eligibility before enrollment. This study was conducted in accordance with the amended Declaration of Helsinki on the ethical conduct of research involving human subjects.

During the study period, all patients visited the obstetrician-gynecology hospital more than once with symptoms of RV. The enrolled patients ranged in age from 20 to 55 years old and were selected based on the following criteria: (a) patients who visited the hospital with vaginitis symptoms more than once in 1 year and (b) patients with apparent vaginitis symptoms (e.g., itching, vaginal discharge, and odor). The exclusion criteria were the use of chemotherapy or any antibiotics, vaginal suppositories, or consumption of hormone products 3 months before the hospital visit. Other exclusion criteria included the presence of systemic diseases (diabetes mellitus, cancer) or a history of total hysterectomy or vaginal surgery.

For comparison with healthy women (n = 100), the 16S rRNA gene sequencing data tagged as Project I from healthy, adult, North American females, collected from the posterior fornix and urogenital tracts, were obtained the Public Database of the Human Microbiome Project (BioProject PRJNA48333; Human Microbiome Project Consortium, 2012).

Samples from the enrolled subjects were collected by a physician using a flocked swab (NFS-2 swab, Noble Biosciences, Inc., Hwaseong-si, Korea) from the posterior fornix of the vaginal wall. The collected sample was suspended in 1 ml of phosphate-buffered saline and immediately transferred to the laboratory for storage at −80°C until further use.

For DNA extraction, a total of 40 vaginal swab samples were thawed on ice, vortexed for approximately 10 s, and centrifuged at 1,455 × g for 10 min. The pellets were resuspended in 100–180 μl of lysis buffer (20 mM Tris–HCl pH 8.0, 2 mM sodium EDTA, 1.2% Triton X-100, 20 mg/ml), containing 20 mg/ml lysozyme from chicken egg whites (Sigma, St. Louis, MO, United States), and incubated for 30 min at 37°C. After 30 min, proteinase K and buffer AL from the Qiagen DNeasy Blood and Tissue kit (QIAGEN, Venlo, Netherlands) were added and incubated for 30 min at 56°C according to the manufacturer’s protocol. DNA concentration and quality were measured using a spectrophotometer (Eppendorf, Hamburg, Germany). DNA concentration of more than or equal to 15 ng/μl with A₂₆₀/A₂₈₀ ratio of 1.8–2.0 was used for 16S rRNA sequencing.

The genomic DNA was amplified by polymerase chain reaction (PCR) using fusion primers 341F (5′-AATGATACGGCGACCACCGAGATCTACAC-XXXXXXXX-TCGTCGGCAGCGTC-AGATGTGTATAAGAGACAG-CCTACGGGNGGCWGCAG-3′) and 805R (5′-CAAGCAGAAGACGGCATACGAGAT-XXXXXXXX-GTCTCGTGGGCTCGG-AGATGTGTATAAGAGACAG-GACTACHVGGGTATCTAATCC-3′; the underlined sequence indicates the complementary region of the primer) binding to the V3–V4 regions of the 16S rRNA gene. The amplification conditions were as follows: initial denaturation at 95°C for 3 min, followed by 25 cycles of amplification (denaturation at 95°C for 30 s, primer annealing at 55°C for 30 s, and extension at 72°C for 30 s), with a final elongation at 72°C for 5 min. PCR amplification and sequencing were performed at ChunLab, Inc. (Seoul, Korea), using the Illumina MiSeq Sequencing System (Illumina, United States) according to the manufacturer’s instructions.

The raw sequence reads were pre-processed for quality checks, such as trimming of low-quality amplicons, MTP, diversity calculation, comparative MTP, and taxonomic biomarker discovery analyses, and analyzed using the 16 s rRNA gene sequencing-based MTP of an EzBioCloud App,¹ a ChunLab bioinformatics cloud platform (ChunLab, Inc., Seoul, Korea). Briefly, pre-processed raw sequences with low-quality (<Q25) reads were filtered using Trimmomatic ver. 0.32. After the quality check, paired-end sequence data were merged using VSEARCH version 2.13.4, followed by trimming of primers, extraction of non-specific amplicons that do not encode 16S rRNA, and detection of chimeric reads. Species-level identification was conducted using the EZBioCloud database. The Operative taxonomic unit (OTU) was picked using the open-reference method with a 97% cutoff similarity. For comparative MTP, copy number normalization (read count, gene copy number) and cutoff (%) of 1.0% for et cetera (ETC) were applied. MTP diversity was calculated using alpha diversity indices, such as ACE, Chao1, and Jackknife for species richness and Shannon, NPShannon, and Simpson for species diversity and evenness. Beta diversity distances were calculated to visualize the sample difference with the principal coordinate analysis (PCoA) by the Jensen–Shannon method. Beta set-significance analysis between two MTP sets was performed using the permutational multivariate analysis of variance (PERMANOVA) test. Additionally, taxonomic biomarkers were determined using the linear discriminant analysis (LDA) effect size (LEfSe) and Kruskal–Wallis H test.

Non-parametric comparison tests (Kruskal–Wallis test for more than two groups and Mann–Whitney U-test for two groups) were used to compare MTP results using MedCalc version 19.8 (MedCalc Software, Ostend, Belgium). Statistical significance was set at p < 0.05.

Forty women were enrolled in this longitudinal study. The study participants were all Korean women, with a mean age of 44.05 ± 7.3 years and a body mass index (BMI) of 23.47 ± 3.41 kg/m². The majority of women (n = 33, 82.5%) were premenopausal. Clinical characteristics of patients with RV are summarized and displayed using descriptive statistics in Table 1. Gynecological symptoms included vaginal discharge (n = 22, 55%), itching (n = 4, 10%), malodor (n = 1, 2.5%), with at least 13 patients reporting two or more symptoms (32.5%). Urological symptoms (increased frequency or urge of urination) were observed in 35 patients (87.5%). Among premenopausal women, 13 (32.5%) had irregular menstrual cycles, and eight had heavy menstrual bleeding (20%). Underlying gynecological diseases [uterine leiomyoma (n = 10), adenomyosis (n = 3), and endometrial polyps (n = 4)] were observed in 17 patients (42.5%).

The total valid reads from patient samples (n = 40) ranged from 3,277 to 79,775 (median = 33,017) after removing low-quality amplicons, non-target amplicons, and chimeras from the 16S rDNA gene sequencing results. The average read length was 420 ± 7 bp per sample. The number of OTUs ranged from 10 to 160 (median = 23), whereas the number of species ranged from 8 to 156 (median = 19) in samples obtained from 40 patients with RV.

The total number of valid reads obtained from the samples of healthy women (n = 100) ranged from 3,320 to 52,645 (median = 6,114). The number of OTUs found in the samples ranged from 8 to 99 (median = 26), whereas the number of species found in the samples ranged from 8 to 98 (median = 25). The overall percentage of valid reads identified at the species level in 40 patients (99.3 ± 2.7%) was similar to that of 100 healthy women (99.3 ± 0.8%).

We compared vaginal MTP and found that the abundance of Lactobacillus spp. was significantly lower in samples from patients with RV (n = 40) than those from healthy women (median 25.80 vs. 99.23, p < 0.001; Figure 1). When 24 representative microbial taxa were analyzed at the species level, the relative abundance of Lactobacillus crispatus, Lactobacillus gasseri, Lactobacillus jensenii, and Lactobacillus reuteri was significantly lower in patients with RV than in healthy women (p < 0.001, 0.044, <0.001, and <0.001, respectively). In contrast, the populations of Gardnerella vaginalis, Atopobium vaginae, Gemella asaccharolytica, Leptotrichia amnionii, Escherichia coli group, and Ureaplasma urealyticum were significantly higher in patients with RV than in healthy women (p < 0.001, p < 0.001, p = 0.005, p < 0.001, p = 0.020, and p < 0.001, respectively).

Furthermore, the species richness, diversity, and complexity (alpha and beta diversity indices) were significantly different between the microbiota of the two groups. The species richness indices, including ACE, Chao1, and Jackknife, were significantly lower in patients with RV than in healthy women (median 24.50% vs. 40.60%, p = 0.011; 17.00% vs. 29.30%, p = 0.003; 21.00% vs. 32.00%, p = 0.086, respectively). In contrast, the NPShannon and Shannon diversity indices were significantly higher in patients with RV than in healthy women (median 0.68 vs. 0.33, p = 0.035; 0.68 vs. 0.31, p = 0.029, respectively), whereas the Simpson index was lower (median 0.55 vs. 0.89, p = 0.012) in patients with RV than in healthy women. These results indicate higher species diversity and evenness in patients with RV than in healthy women. Beta diversity distance was also significantly different between patients with RV and healthy women (p = 0.001, pair-wise PREMANOVA with 999 permutations). A comparative analysis of PCoA plots between RV and healthy women is outlined in Figure 2, based on the Jensen–Shannon beta diversity distances at the species level.

Based on the Lactobacillus spp. dominance, we taxonomically classified 40 vaginal microbiome compositions into six CST groups: L. crispatus-dominant CST group I (n = 5, 12.5%); L. gasseri-dominant CST group II (n = 2, 5.0%); Lactobacillus iners-dominant CST group III (n = 8, 20.0%); CST group IV with low abundance of Lactobacillus spp. and increased diversity and abundance of anaerobic species (n = 21, 52.5%); L. jensenii-dominant CST group V (n = 1, 2.5%); mixed CST group with co-dominance of more than two Lactobacillus spp. (n = 3, 7.5%; De Seta et al., 2019). The detailed microbiome taxonomic composition is shown in Supplementary Table 1.

The relative abundance (%) of 24 representative microbial taxa at the species level was presented on a heatmap. Hierarchical clustering (based on unweighted pair group method with arithmetic mean, UPGMA) showed accurate clustering of the relative abundance according to the CST group (Figure 3).

On comparing the relative abundance of different bacterial taxa among the CST groups, we observed that the relative abundance of Lactobacillus spp. was significantly lower in the CST group IV than in the other CST groups (p < 0.01; Figure 4). The populations of L. crispatus, L. iners, L. jensenii, L. reuteri, Bifidobacterium scardovii, and G. vaginalis were also significantly different among the CST groups (p < 0.001, <0.001, <0.001, 0.004, <0.001, and 0.003, respectively). The statistical significance of the relative abundance of bacterial taxa between the two groups was displayed as ^*p < 0.05 and ^**p < 0.01 in Figure 4.

The LEfSe analysis used to determine the relative abundance of bacterial taxa (LDA score > 2.0; p < 0.05) showed significant statistical enrichment. The phylum Bacteroidetes and lower taxa, including the order Bacteroidales, phylum Actinobacteria, and genera KQ959671_g (Coriobacteriaceae) and Gardnerella, and the phylum Firmicutes and lower taxa including family Peptoniphilaceae were significantly enriched in CST group IV. In contrast, the phylum Firmicutes and lower taxa, including genus Lactobacillus spp., were significantly enriched in CST groups I and III (Figure 5).

In our study, alpha diversity indices associated with species richness and diversity showed variations depending on the presence or absence of uterine diseases (Table 2; Figure 6). Patients with uterine diseases (n = 17) showed higher species richness indices of Chao and Jackknife (p = 0.006 and 0.003, respectively) and species diversity indices of Shannon, NPShannon, and Simpson (p = 0.035, 0.034, and 0.043, respectively) than those without diseases (n = 23).

We further found that the Simpson index was also lower in patients going through menopause (n = 7; p = 0.039), indicating increased species diversity in menopause status.

However, no significant association was found between the CST groups and other variables, such as uterine diseases, menopause, urological symptoms, gynecologic symptoms, menstrual cycle, menorrhagia, and BMI (data not shown).