This study was conducted in the Microbiology Institute at the Universidad San Francisco de Quito (USFQ) from June 2017 to November 2018. As previously reported (Salinas et al., 2020), 436 Ecuadorian women of Hispanic ethnicity between 18 and 56 years old volunteered to be part of the epidemiological study. Briefly, all women received a kit containing an informed consent approved by the Bioethics Committee of the USFQ, a standardized medical survey, and a vaginal transport swab system (Stuart’s transport media swabs; Copan Diagnostics Inc., Brescia, Italy). Volunteers were excluded if they reported having had sexual intercourse within the last 48 h, antimicrobial treatment in the last 3 months, or any evidence of bleeding. The study was supervised by a physician, a psychologist, and a full-time researcher from the USFQ. This investigation adopted a cross-sectional study design to determine the association between the presence of five well-known lactobacilli species and vaginal microbiota or opportunistic pathogens (such as Gardnerella spp., Fannyhessea vaginae previously known as Atopobium vaginae, Mobiluncus spp., Escherichia coli, and Candida albicans) previously diagnosed/detected in our last publication (Salinas et al., 2020), more exactly, healthy microbiota, intermediate microbiota, and vaginal dysbioses (AV and BV, and VC).

The present study was approved by the Ethics Committee of the USFQ (Protocol code: 2016-023IN by MSP-VGVS-2016-0244-O review board).

DNA extraction was realized through standard procedure following Peng and colleagues’ direct boiling point method (Peng et al., 2013). Briefly, the stored aliquots (0.9% NaCl) of 1 ml were incubated at 100°C in a water bath for 15 min and then immediately frozen at −20°C for 15 min. Next, the samples were centrifuged at 13,000 rpm for 15 min, and supernatants were aliquoted into 500-μl volumes. DNA quantification was performed with a NanoVue spectrophotometer (GE Healthcare Life Science, Marlborough, MA, USA), samples were eluted at 20 ng/µl with molecular grade water and stored at −20°C until the PCR analysis was performed. The quality of DNA was evaluated by measuring the concentration of phenolic compounds or the presence of salts (260/230) and protein contaminants (260/280). This procedure was adapted from Money (Money, 2005).

From our previous study (Salinas et al., 2020), 436 vaginal samples were selected for molecular characterization by PCR in a Bio-Rad Thermocycler (Bio-Rad, Hercules, CA, USA). DNA quantification was performed with a NanoVue spectrophotometer (GE Healthcare Life Science) to ensure the presence of amplifiable DNA. Concentrations of DNA in ng/μl were measured, as well as the phenolic contaminants (260/230) and the protein contaminants (260/280). Aliquots of DNA between 10 and 20 ng/µl were used for PCR analysis. Before lactobacilli detection was realized, all samples were analyzed for 16S conserved rRNA genes (fDD2-CCGGATCCGTCGACAGAGTTTGATCITGGCTCAG; rPP2-CCAAGCTTCTAGACGGITACCTTGTTACGACTT) by PCR, ensuring the absence of PCR inhibitors on samples, as previously described (Borja-Serrano et al., 2020). All samples were analyzed with a total of five primer pairs, targeting five Lactobacillus species (L. acidophilus, L. crispatus, L. gasseri, L. jensenii, and L. iners). Single-template PCR assays were performed for each primer set. The sequence, amplicon size, target species, and temperature of annealing for each primer pair are described in Supplementary Table 1 . A final volume of 20 µl was used according to the reference protocols (Galán et al., 2006; Fredricks et al., 2007; Sepehri et al., 2009; Henriques et al., 2012; DTU- National Food Institute, 2014), which included 0.5 U of Go Taq^® DNA Polymerase (Promega, Madison, WI, USA), 1× of Green GoTaq^® Flexi Buffer (Promega), 0.25 mM of MgCl₂ (Promega), 200 µM of dNTP mix (Promega), 0.5 µM of each primer and target template DNA concentration of approximately 4 ng/μl, and the remaining volume with molecular grade H₂O. The PCR thermal cycling consisted of initial denaturation at 94°C for 2 min, followed by 29 cycles of denaturation at 94°C for 30 s, annealing at each primer pair temperature for 30 s and extension at 72°C for 1 min, and final extension of 5 min at 72°C. The respective use of negative (without DNA sample and samples with other related bacteria) and positive (collection of identified strains of each species through DNA sequencing) controls were used in each PCR assay. These positive controls were provided by the Microbiology Institute at USFQ. All samples were randomly performed in triplicate with different negative and positive controls. After PCR amplification, a volume of 4 µl from each PCR product was visualized in 1.5% (w/w) agarose (Promega) gel electrophoresis using 0.1% ethidium bromide staining. The DNA analysis was performed under permit No. MAE-DNB-CM-2016-0046 (De Backer et al., 2007; Garg et al., 2009; Tsai et al., 2010; Zhang et al., 2012).

A multivariate logistic regression model was used to calculate the odds ratios (OR) of the clinical outcomes that included demographic variables (age, sex, city, and marital status), socioeconomic variables (occupation and level of education), personal habits (sex relationships, hygiene, and other habits), and the type or number of vaginal Lactobacillus species associated with the presence or absence of vaginal infection using logistic regression. These data were also considered categorical variables. Firstly, the variable of vaginal infection in the samples was categorized as the presence and absence, so a comparison of the different risk factors of both groups can be performed. After further statistical analysis, the study was defined by the type of vaginal dysbiosis (BV and AV) for testing differences in the previously analyzed factors and vaginal microbiota. The chi-square test was used to evaluate associations between the prevalence of vaginitis with the other risk factors. A value of p < 0.05 and 95% CIs were considered significant for the test. Logistic regression was also performed to calculate crude ORs for each variable mentioned; adjusted ORs were produced for variables with statistical significance in both tests applied for association (Ozaydin et al., 2013; Porras et al., 2014; Syam et al., 2015). Therefore, the chi-square test was used as a test of association, while the OR was then used as a measure of association (Kim, 2017). The statistical analysis of association with risk factors was performed for each type of vaginal infection but negative for the remaining types of vaginal infection to observe a significant difference between those populations. Each type of vaginal infection, normal or healthy microbiota, and intermediate microbiota were classified as dependent variables against sociodemographic and behavioral variables or the presence of Lactobacillus species as independent variables. All initial values of p < 0.05 obtained by univariable logistic regression, chi-square, and multivariable logistic regression analyses were then evaluated through the BH adjustment to detect false discovery rate (FDR) for conducting multiple comparisons. All statistical analyses were performed using SPSS version 22.0 (SPSS Statistics for Windows Version 22.0, IBM Corp, Armonk, NY, USA), except for the BH adjustment. The BH adjustment was realized using Seed-based d Mapping software (SDM, version 6.21, https://www.sdmproject.com, formerly “Signed Differential Mapping”) (Radua and Mataix-Cols, 2009; Radua et al., 2012). A clustering model was realized through Ward’s Minimum Variance Clustering Method with Euclidean squared distance to perform hierarchical clustering fed by a dimensionality reduction algorithm Principal Component Analysis (PCA) implemented in RStudio software (version 1.3.1073; https://rstudio.com/), using the option method = “ward” of the hclust function from the stats base R package (Package stats version 4.1.0) (Murtagh and Legendre, 2014).

A total of 436 women volunteered in our last study (Salinas et al., 2020), with their vaginal samples and epidemiologic data selected for lactobacilli characterization in the present study. The stored samples were chosen for the molecular analysis by PCR. As shown in Table 1 , our population set was constituted by Ecuadorian women between 18 and 56 years old, with 76.3% between 18 and 28 years old. Briefly, 66.1% of the women have healthy vaginal microbiota, 10.8% have an intermediate microbiota, and finally, 23.1% showed vaginal dysbiosis or infections. Among women with vaginal dysbiosis or infection, AV was the main vaginal dysbiosis being diagnosed in 52.5% (53/101), followed by BV (23.8%; 24/101) and VC (6.9%; 7/101). Eighty-four women were diagnosed with a single type of dysbiosis (83.2%), and the remaining 17 had vaginal coinfections (16.8%). The most common coinfections found in women were BV and AV (12/17), followed by BV and VC (3/17), AV and VC (1/17), and lastly, all studied vaginal infections (1/17).

Approximately 87.2% of the population set was constituted by undergraduate students or young professionals (380/436). The categories of professionals included the following: health professionals (23.0%), administrative clerks (20.3%), educators (14.9%), and general employees with college degrees (18.0%). The remaining professions without specialization or need for college degrees were classified as unprofessional careers. Most of the volunteers were single women (78.7%) and of Hispanic ethnicity (88.3%). Among the participants, 56.4% had a steady sexual partner, and 38.5% reported not having any sexual partner. Concerning birth control methods, 17.0% of participants reported using a condom, 40.8% reported the use of other birth control methods, and the remaining women did not use any birth control method (34.9%) or merely did not answer (7.3%). Alternative birth control methods included hormone treatment or other forms of protection (e.g., spermicides, diaphragm, cervical cap, and sterilization), intrauterine device (IUD), and natural methods (abstinence, fertility awareness method (FAM), and withdrawal). In our study, the most used alternative contraceptive method was hormonal, through oral contraceptives (46.7%) and local implants (6.6%).

When the results of the chi-square test in the prevalence of each Lactobacillus sp. between healthy microbiota, intermediate microbiota, and vaginal dysbiosis were analyzed, each group showed statistically significant differences in the presence of L. acidophilus, L. jensenii, and L. gasseri, as shown in Figure 1A . However, only L. acidophilus showed simultaneous statistical differences between healthy and intermediate microbiota (p = 0.026) and between healthy microbiota and vaginal dysbiosis (p = 0.015). The prevalence of both L. jensenii and L. gasseri was statistically different between healthy microbiota and vaginal dysbiosis (p = 0.017 and p = 0.020, respectively). No statistically significant differences were observed in these lactobacilli presence between intermediate microbiota and vaginal dysbiosis (see Figure 1A ). Finally, L. iners and L. crispatus did not demonstrate statistical differences among these groups of the vaginal microbiota.

Further statistical analysis was then realized between healthy microbiota and specific types of vaginal dysbiosis (more exactly AV and BV), as well as between vaginal dysbioses. As shown in Figure 1B , some statistically significant differences were found on certain Lactobacillus species when comparing healthy microbiota against BV and AV, but no statistical differences were found between BV and AV. L. acidophilus showed statistically significant differences in its prevalence on healthy microbiota against both dysbioses (BV, p = 0.041; and AV, p = 0.045), while L. jensenii only showed statistically significant differences between healthy microbiota and AV cases (p = 0.012).

To evaluate if the statistically significant differences found in the prevalence of lactobacilli could have a protective effect against the development of vaginal dysbiosis, univariable logistic regression analyses were then performed. As shown in Figure 2A , each Lactobacillus species was normalized according to the importance of their presence against the vaginal dysbiosis establishment, showing L. acidophilus, L. jensenii, and L. crispatus importance of 100%, 79.3%, and 74.8%, respectively. However, only L. acidophilus and L. jensenii exhibited statistically significant differences (p = 0.035 and 0.050, respectively), suggesting a potential protective effect against the development of vaginal dysbiosis.

In addition, little is still known about epidemiological factors and lactobacilli colonization among women (Vaneechoutte, 2017a; Auriemma et al., 2021). Therefore, multiple chi-square analysis was also performed to evaluate possible correlations with each individual Lactobacillus species. As shown in Figure 2B , L. iners (p = 0.017) and L. gasseri (p = 0.013) were more prevalent in women with previous antimicrobial treatment in their clinical background. L. crispatus (p = 0.048) was associated with the presence of vaginal secretion among women. However, L. iners was also related to other epidemiological factors, such as occupation (p = 0.009) and education level (p = 0.013), showing statistical differences in its distribution among women in these categories. More exactly, a higher prevalence of L. iners was found in women with a secondary level of education (see Table 1 ). Finally, other factors, such as having a sexual partner and contraceptive use, demonstrated statistically significant values in relation to multiple Lactobacillus species (L. acidophilus, L. gasseri, L. crispatus, and L. iners) differing only in the absence of L. iners in contraceptive use. Although these results evaluated the species’ individual probiotic role in the vaginal microbiota, it is well known that a probiotic microbiota is characterized by a multi-microbial effect character and is not caused merely by an individual effect (Vaneechoutte, 2017b; Wieërs et al., 2020). Therefore, a multi-microbial analysis was performed to evaluate a potential symbiotic or antagonistic relationship between these Lactobacillus species against both cases of vaginal dysbiosis.

Nowadays, it is well known that the probiotic activity provided by a certain microbiota is caused not just by the effect of an individual Lactobacillus species but also by its multi-microbial interaction. So further analysis was also done through the clustering model of these Lactobacillus species against vaginal dysbiosis (by itself and then AV and BV) and epidemiological factors. A clustering model was realized through Ward’s Minimum Variance Clustering Method evidencing multiple clusters. The clustering of samples was developed according to the presence of different Lactobacillus species in vaginal samples. As shown in Figure 3A , six clusters were selected for multiple chi-square analysis to evaluate statistically significant differences (p ≤ 0.05). Cluster 1 was characterized by the presence of L. iners and L. jensenii, while Cluster 2 was only formed by L. iners. Cluster 3 was constituted by L. iners, L. jensenii, and L. gasseri. Cluster 4 assembled four Lactobacillus species, more exactly, L. iners, L. jensenii, L. gasseri, and L. acidophilus. Finally, Cluster 5 gathered all studied Lactobacillus species, and Cluster 6 evidenced no lactobacilli presence.

As shown in Figure 3B , only Clusters 3 and 4 demonstrated statistically significant differences against Clusters 1 and 6 in the establishment of vaginal dysbiosis. Clusters 3 and 4 shared L. iners, L. jensenii, and L. gasseri, but Cluster 4 also comprised L acidophilus. Both clusters evidenced a multi-species effect, being more notorious in Cluster 4 due to the obtained statistical values (p < 0.030). Interestingly, Cluster 1 also gathered L. iners and L. jensenii as Clusters 3 and 4; however, the absence of L. gasseri and L. acidophilus led to the lack of probiotic protection in vaginal dysbiosis establishment. The absence of Lactobacillus species in Cluster 6 was expected to relate to vaginal dysbiosis. It is also important to highlight that no statistical differences were found between Cluster 3 and 4 or even between Cluster 1 and 6, suggesting a potential probiotic connection among them. When individually evaluating each vaginal dysbiosis, greater statistically significant differences were found in the presence of AV (three p-values ≤0.05) than BV (one p-value ≤0.05). In AV cases, clustering analysis showed statistically significant differences in Clusters 3 and 4 when compared to Cluster 6 (p = 0.022 and p = 0.020, respectively), but only Cluster 3 showed statistical difference against Cluster 1 (p = 0.048). However, in BV cases, only Cluster 4 showed a statistically significant difference against Cluster 6 (p = 0.009), which is characterized by the absence of Lactobacillus species. In AV and BV cases, the combination of L. iners, L. jensenii, L. gasseri, and L. acidophilus from Cluster 4 reflected a multi-microbial consortium with statistical differences in the establishment of both dysbioses.

Finally, multiple chi-square analysis was also performed to evaluate the epidemiological factors related to the presence of each Cluster. As shown in Figure 3C , only Clusters 2 and 6 showed statistically significant differences among epidemiological factors. Both clusters shared statistically significant differences in women with a sexual partner and no contraceptive use (see Figure 3C ), evidencing an association between these epidemiological behaviors and the lack of lactobacilli apart from L. iners. However, only Cluster 2 was associated with vaginal secretion in women (p = 0.048), more exactly, the presence of L. iners and the absence of the remaining analyzed lactobacilli, while Cluster 6 was correlated to women with a previous clinical history of antibiotic treatment for vaginal dysbiosis (p = 0.023), suggesting a potential correlation between treatments and eradication of vaginal lactobacilli.

To determine the relationship between the presence of any opportunistic pathogens and each Lactobacillus species analyzed in this study or the clusters formed by them, multiple chi-square tests were realized (see Supplementary Tables 2 , 3 ). The opportunistic pathogens previously identified in our last study were Gardnerella spp., F. vaginae (previously known as A. vaginae) and Mobiluncus spp. related with BV, E. coli related with AV, and C. albicans related with candidiasis (Salinas et al., 2020). As shown in Supplementary Table 2 , every opportunistic pathogen showed a statistical significance for at least one Lactobacillus species, except for Gardnerella genus.

Interestingly, the presence of L. iners was statistically correlated with the absence of Mobiluncus species (p = 0.033). On the contrary, the absence of L. jensenii was statistically associated with the absence of C. albicans (p = 0.034), while the absence of L. acidophilus evidenced the same association with F. vaginae (p < 0.001) and E. coli (p = 0.015). Meanwhile, L. crispatus showed multiple statistical associations with F. vaginae, C. albicans (both p < 0.001), and E. coli (p = 0.005), where its absence was correlated with the absence of these opportunistic microorganisms. Finally, no statistical correlation was found in the presence or absence of L. gasseri in the vaginal epithelium among Ecuadorian women.

To better understand the multispecies probiotic activity of lactobacilli, multiple chi-square tests were further studied in the lactobacilli clusters. The results showed significant values among clusters 1, 2, 4, and 6, as shown in Supplementary Table 3 . As previously stated, these clusters represent the presence of L. iners and L. jensenii (Cluster 1), only L. iners (Cluster 2), all lactobacilli except for L. crispatus (Cluster 4), and the absence of Lactobacillus sp. (Cluster 6). As shown in Supplementary Table 3 , Cluster 1 evidenced an inhibition of the presence of Gardnerella genus (p = 0.04). Cluster 2 is particularly interesting showing multiple statistical associations, illustrating a proliferation of F. vaginae (p = 0.001) and an inhibition of C. albicans (p = 0.001) and E. coli (p = 0.006). Mobiluncus spp. showed an opposite effect with the presence of Cluster 4, while the absence of Lactobacillus sp. in Cluster 6 could be inhibiting the proliferation of F. vaginae (p = 0.05).