A total of 90 women were recruited at the Gynecology Outpatient Department of Sichuan University West China Second Hospital in Chengdu between January 2017 and May 2021. Forty-five patients diagnosed with VLSC confirmed pathologically were included as the LSC group. Forty-five women with no recorded vulvar complication and vaginitis were included as the control group (H group). The excluding criteria were (1) the vaginal smear indicated vaginitis such as bacterial vaginosis, vulvovaginal candidiasis or trichomoniasis, etc; (2) patients who had sexual intercourse within 48 hours before the visit; (3) patients received vaginal medications or douches within 48 hours before the visit; (4) pregnant women; (5) patients had taken any oral antibiotics,oral antimycotics, oral contraceptive medication or hormone replacement treatment in the past one month; (6) patients had received topical corticosteroids, topical immunomodulators, antihistamines, antidepressants or transcutaneous electrical nerve stimulation in the past one month. After obtaining informed consent, the mid-vaginal secretion samples were dipped using sterilized cotton swabs. Two sterilized swabs were collected for each sample and the swabs were put in tubes on ice until genomic DNA extractions within 3 hours. This study was approved by the Institutional Ethics Committee of Sichuan University West China Second University Hospital (NO.1356) and the work has been carried out in accordance with The Code of Ethics of the World Medical Association (Declaration of Helsinki) for experiments involving humans.

The TIANamp Swab DNA kit (TIANGEN BIOTECH Beijing Co. Ltd, Beijing, China) was used to extract the whole genomic DNA of vaginal bacteria species. The swab suspensions were mixed with 20 $\mu$L proteinase K. The mixture was incubated at 56 °C for 60 min. Then, the samples were subjected to mechanical lysis by bead beating in a Fast Prep 24 machine (MPBio, Santa Ana, CA, USA) at 6 m/s for 40 s. The lysates were subjected to centrifugation at 12,000 rpm for 60 s to pellet the beads and filtered using the CR2 columns. DNA was purified, and the quality of total genomic DNA was evaluated by agarose gel electrophoresis (Invitrogen, Carlsbad, CA, USA).

The V3-V4 variable regions of the bacterial 16S rRNA gene were amplified with primers 338F (ACTCCTACGGGAGGCAGCA) and 806R (GGACTACNNGGGTWTCTAAT). The amplicon quality was visualized using gel electrophoresis. The PCR products were purified with AgencourtAMPure XP beads (Beckman Coulter Co., Bria, CA, USA) and quantified using Qubit dsDNA assay kit. The concentrations were then adjusted for sequencing. Sequencing was performed on an Illumina NovaSeq6000 with two paired-end read cycles of 250 bases each (Illumina Inc., San Diego, CA, USA; OE Biotech Company; Shanghai, China)

The raw reads were deposited into the Genome Sequence Archive (GSA) database (Accession Number: CRA005753. Paired-end reads were preprocessed using Trimmomatic software (version 0.35, Usadel Lab, Aachen, Germany) to detect and cut off ambiguous bases (N). It also cut off low quality sequences with average quality score below 20 using sliding window trimming approach. After trimming, paired-end reads were assembled using FLASH software (version 1.2.11, Center for Computational Biology, Baltimore, MD, USA). Parameters of assembly were 10 bps of minimal overlapping, 200 bps of maximum overlapping and 20% of maximum mismatch rate. Sequences were performed further denoising as follows: reads with ambiguous, homologous sequences or below 200 bps were abandoned. Reads with 75% of bases above Q20 were retained and reads with chimera were detected and removed using QIIME software (version 1.8.0, Knight and Caporaso labs, Flagstaff, AZ, USA). Clean reads were subjected to primer sequences removal and clustering to generate operational taxonomic units (OTUs) using VSEARCH software (version 2.4.2, Institut Pasteur, Paris, France) with 97% similarity cutoff. The representative read of each OTU was selected using QIIME package. All representative reads were annotated and blasted against Silva database (Version 138) using RDP classifier (confidence threshold was 70%). The community structure of LSC and H groups were compared. Venn diagrams were drawn to analyze the amount of overlapped and unique OTUs in the two groups. In the alpha-diversity analysis, Shannon and inverse Simpson indices were estimated. In the beta-diversity analysis, principal co-ordinate analysis (PCoA) using binary_jaccard was performed. Linear discriminant analysis effect size (LEfSE) was performed to determine the dominant biomarkers of each group. Phylogenetic investigation of communities by reconstruction of unobserved states (PICRUSt) was performed to predict the function of different vaginal microbiota. Student’s t-test was used to compare quantitative data and p 0.05 was considered statistically significant.

The average age of the LSC and H groups was 44.89 $\pm$ 11.52 and 42.09 $\pm$ 10.3 years respectively (p $>$ 0.05). The mean BMI of the two groups was 23.95 $\pm$ 3.03 and 22.78 $\pm$ 3.21 kg/m${}^2$ (p $>$ 0.05). There was no significant difference in marital status within the two groups.

The data volume of valid tags ranges from 9359 to 74,949, and the length of the valid tags varied from 411.97 to 430.01 bps. 4475 OTUs were identified using a cutoff of 97% sequence similarity, including 31 phyla, 600 genera, and 905 species.

At the phylum level, the Firmicutes, Actinobacteria and Proteobacteria were the top three abundant bacterial phyla in H group. In the LSC group, the dominant three phyla were Firmicutes, Actinobacteria and Bacteroidota (Fig. 1A). At the genus level, the Lactobacillus and Gardnerella were the most abundant two bacterial genera in both of the two groups, with relative abundance of 64.32% and 13.64% respectively in the H group, and 64.82% and 9.36% respectively in the LSC group. Among the Lactobacillus, it was found that the Lactocillus iners. was the most abundant species. The relative abundance of Gardnerella decreased while Muribaculaceae and Streptococcus increased in the LSC group (Fig. 1B). The Venn diagram displayed that there were 2631 common OTUs of the two groups. H group had 402 unique OTUs, and LSC group had 1442 unique OTUs. The microbiota community of each group was different (Fig. 1C).

Fig. 1.

Microbial community structure bar plots of the LSC and H groups. (A) The microbial community bar plots of the LSC and H groups at the phyla level. (B) The microbial community bar plots of the LSC and H groups at the genus level. (C) OTU Venn diagram. The number in the overlapped area indicates the amount of common OTUs of the two groups.

The median of Shannon index of LSC group was 0.79 which was significantly less than the index of H group (1.43, p 0.01) (Fig. 2A). The median of inverse Simpson index of LSC group was 0.24 while the index of the H group was 0.42, and there was statistical difference between the two groups (p 0.05) (Fig. 2B). This indicated the alpha diversity of LSC group was less than the H group. The beta diversity of the two groups was compared using the principal co-ordinates analysis (PCoA) of binary_jaccard distances. The result was distinct between LSC and H groups (Fig. 3).

Fig. 2.

Alpha diversity comparison of LSC and H groups. Violin plots of the Shannon (A) and Simpson (B) index for the two groups. ‘**’ indicates p 0.01, ‘*’ indicates p 0.05.

Fig. 3.

PCoA analysis of the LSC and H groups. Some of the H group samples cluster separately from the LSC samples.

The analysis of similarity (anosim) of binary_jaccard revealed LSC group differed from H group significantly (p = 0.001). The top 5 discriminal genus with most relative abundance were Gardnerella, Atopobium, Tepidimonas, Sneathiaand and Massllia (Fig. 4). To identify the specific taxa, the LEfSE analysis was performed. Comparing with group H, family Leptotrichiaceae and genus Sneathia were discriminant taxa in group LSC (Fig. 5).

Fig. 4.

The box plot of the top 10 different genera in the LSC and H groups.

Fig. 5.

LEfSE analysis of the LSC and H groups. The rings represent the taxonomic levels in descending order, from phylum to genus. The circles in red or green represent more prevalent taxa of H or LSC groups. The diameters of the circles represent the relative abundances of the taxa.

To gain the function of vaginal microbiota in the etiology of VLSC, PICRUSt using Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways was used to predicate the functional profile of bacterial communities. It found that in addition to the difference detected in bacterial composition and diversity, microbial genes related to the signal transduction, metabolism of terpenoids and polyketides, transporters, excretory system, nervous system, energy metabolism and many other functions were extremely enriched among the LSC group comparing with H group, while genes related to immune system, cardiovascular diseases and digestive system, etc were reduced in LSC group (Fig. 6).

Fig. 6.

Heatmap summarizing the predicted functional factors based on KEGG databases in the level 2.