Vaginal Microbiome Signature Is Associated With Spontaneous Preterm Delivery

Background: Preterm delivery (PTD) represents an important public health and therapeutic challenge. Despite the reported link between the composition of vaginal microbiome and PTD, previous studies were inconsistent in their conclusions and utilized non-uniform designs. We performed an independent case-control study carried out on the Slovenian population, where we re-evaluated the role of the vaginal microbiome in PTD. Methods: Vaginal microbiomes of pregnant women who delivered preterm were compared to those delivered at term to examine differences in the microbial richness, diversity, and differential abundance of specific taxa. We obtained vaginal swab samples from 155 Caucasian women who were classified as either term (≥380/7 weeks, n = 107) or preterm (≤366/7 weeks, n = 48) in exclusion of any other medical or obstetric conditions. The vaginal microbiomes of these women were characterized by 16S ribosomal RNA (rRNA) gene sequencing of the V3-V4 region on the MiSeq platform. Results: Women who experienced PTD had a higher microbial richness (Chao1, P = 0.011) and alpha diversity (Shannon, P = 0.00059) than women with term deliveries. We report that overall vaginal microbial community composition (beta-diversity) was significantly different by delivery gestational age category (PWeightedUnifrac < 0.001). Women who delivered preterm had decreased Lactobacilli spp. abundance as well as increased abundance of Gardnerella and other bacterial vaginosis (BV) and aerobic vaginitis (AV) associated genera including Atopobium, Sneathia, Gemella, Megasphaera, Dorea, Streptococcus, and Escherichia/Shigella. Conclusions: In the present study, we provide further evidence that vaginal microbiome composition is associated with PTD.


INTRODUCTION
Preterm delivery (PTD), defined by the World Health Organization (WHO) as the delivery of an infant before completed 37 weeks of gestation (1), represents a major global public health burden. It is the leading determinant of neonatal morbidity and mortality in developed countries (2) and the second most common cause of death in children under the age of five worldwide (3). Annually, 5-18% of pregnancies end preterm, resulting in 15 million premature babies (3), and prevalence is still escalating. Moreover, survivor premature babies have an increased risk of neonatal complications associated with the immaturity of organ systems, such as respiratory distress syndrome, cerebral palsy, retinopathy, hearing impairment, and bronchopulmonary dysplasia as well as increased long-term risk for a spectrum of cognitive and neurologic disorders (4)(5)(6).
Underlying causes of PTD are not completely understood, with more than 65-75% spontaneous PTDs with idiopathic onset (2,7). Thus, it is still challenging to identify women at higher risk. Epidemiologic data suggest that complex genetic, environmental and sociodemographic factors are involved. Spontaneous PTD is a syndrome of numerous pathological processes (3), including certain maternal or fetal conditions, such as preeclampsia, intrauterine infection, cervical insufficiency, uterine anomalies, preterm premature rupture of membranes (PPROM), polyhydramnios, and fetal malformations (8). Previously observed risk factors also include, prior PTD, multifetal gestation, young or advanced maternal age, assisted reproductive technology (ART), black race, smoking, extremes of body-mass index (BMI), and low socioeconomic status (2).
Multiple lines of evidence support the role of maternal vaginal microbiota as a potential risk factor for PTD. Most notably, as much as 40-50% of premature deliveries are assumed to be due to maternal infection and inflammation (9)(10)(11). Prior research has shown that the most frequent cause of bacterial infection represents ascending pathway from the vagina through the cervix and into the uterus (3), where bacteria provoke deciduitis and chorioamnionitis or enter into the amniotic cavity and the fetus. Immediately after passing amniotic membranes, bacteria and their products are detected by pattern recognition receptors, such as toll-like receptors (TLRs), which causes the production of pro-inflammatory cytokines, chemokines, prostaglandins, and proteases leading to myometrial contractions. Additionally, Abbreviations: AA, African American; ART, assisted reproductive technology; ASV, amplicon sequence variant; AV, aerobic vaginitis; BLAST, Basic Local Alignment Search Tool; BMI, body-mass index; BV, bacterial vaginosis; CST, community state type; DADA2, Divisive Amplicon Denoising Algorithm 2; EPTD, early preterm delivery; FDR, false-discovery rate; GDM, gestational diabetes mellitus; GLM, generalized linear model; HIV, human immunodeficiency virus; IUGR, intrauterine growth restriction; LPTD, late preterm delivery; MDS, multidimensional scaling, MMP-8, matrix metalloproteinase-8; OTU, operational taxonomic unit; PBS, phosphate-buffered saline; PCoA, Principal Coordinates Analysis; PCR, polymerase chain reaction; PERIS, Perinatal Information System of Slovenia; PERMANOVA, permutational multivariate analyses of variance; PPROM, premature rupture of membranes; PROM, premature rupture of membranes; PTD, preterm delivery; RDP, Ribosomal Database Project; rRNA, ribosomal RNA; spp, species; TLRs, toll-like receptors; WHO, World Health Organization.
inflammatory cascade leads to remodeling and disruption of fetal membranes (3,11,12). Hematogenous dissemination has also been proposed as an alternative route of intrauterine infection, as bacteria involved in periodontal disease have been found in the amniotic fluid (13). Furthermore, in utero microbial colonization of fetus has been suggested to have a role in the early establishment of the infant's microbiome (14,15).
During pregnancy, a rise in the level of progesterone and estrogens, as well as immunological changes, significantly remodel composition of the human vaginal microbiome toward stable, low richness and diversity community structures dominated by Lactobacillus spp., which provide greater resistance and a protective role against genital tract infections (16,17). Although dysbiosis (disruption of the vaginal microbiota) is asymptomatic in about 40% of affected women, it may predispose to a higher risk of urogenital conditions including aerobic vaginitis (AV) and bacterial vaginosis (BV) (18). The latest was reported to increase the risk of PTD ∼5-fold (19).
The link between bacterial diversity and PTD has previously been demonstrated (20). Additionally, recent advances in the culture independent-techniques, using next-generation sequencing technologies proposed the causative role of overall vaginal microbiota composition in the development of PTD, although findings across studies have been inconsistent and largely ethnically dependent (7,(21)(22)(23)(24)(25)(26)(27)(28)(29). Therefore, the role of the vaginal microbiome in the etiology of PTD remains to be delineated. However, studies (7,17,25,30) have unequivocally agreed that racial differences in the composition of vaginal microbiome play an important role in the detection of associations. Additionally, different conclusions about the role of vaginal microbiota in the development of PTD may have been reached due to the following reasons: (1) various definitions of PTD, (2) inclusion criteria varied widely (from low-risk PTDs to medically indicated and high-risk PTDs), (3) small power, (4) diverse designs, including different sequencing technologies and different 16S regions sequenced.
Although most of the previous 16S sequencing studies have identified a degree of association between the vaginal microbiome and PTD, there are inconsistencies in results that cannot currently be attributed to either different methodologies, heterogeneous definitions of PTD, or actual variations between study populations. Therefore, we used high-throughput 16S rRNA sequencing of the V3-V4 region in an independent Caucasian population to characterize whether the composition of maternal vaginal microbiota prior delivery differs between women who subsequently delivered preterm and those who delivered term.

Ethics Statement
All the women provided written informed consent prior to recruitment. This study was approved by the Slovenia National Medical Ethics Committee (approval number 90/02/15) and was conducted according to the principles expressed in the Declaration of Helsinki.

Sample Collection and Clinical Assessment
We included 48 pregnant, Caucasian women, with spontaneous onset of premature labor with or without PPROM (≤36 6/7 weeks) (cases) and 107 pregnant women with a spontaneous onset of labor with or without premature rupture of membranes (PROM) or planned cesarean section at term (≥38 0/7 ) (controls). The women gave birth at the Department of Perinatology of the University Medical Center Ljubljana in the period from October 2015 to January 2019. We assessed the gestational age at delivery by last menstrual period and confirmed it through ultrasound examination in the first trimester by the measurement of the fetal crown-rump length (31). Vaginal swabs of the lateral wall of the vagina were collected by experienced obstetricians of the women who exhibited signs of labor based on the presentation of uterine contractions and cervical dilation at the admission to the delivery department.
National Perinatal Information System of Slovenia (NPIS) was reviewed to collect additional data on the pregnancy, delivery and postpartum. Women participating to this study were also asked to complete a brief in-person questionnaire to collect additional demographic and medical information (age, working status, cigarette smoking, pre-pregnancy BMI, antibiotics use during pregnancy, gestational diabetes, previous PTD, and vaginitis.

Bacterial DNA Isolation, 16S rRNA Gene Amplification and Sequencing
The vaginal swabs for genomic bacterial DNA isolation were transferred to the laboratory in the Amies transport medium and stored at −80 • C until isolation. The enzymatic and mechanical cell lysis was performed as described by Ravel et al. (30). The lysate was processed using the QIAamp DNA Mini Kit (QIAGEN, Hilden, Germany), according to manufacturer's guidelines for DNA purification from buccal swabs (recommendations for cotton or Dacron swabs) with minor modifications, including incubation at 70 • C for 10 min instead of 56 • C and 5-min incubation step at room temperature after the addition of 100 µl of heated (56 • C) elution buffer. The DNA concentrations in the samples were quantified using the Qubit dsDNA high-sensitivity Assay Kit (Thermo Fischer Scientific). After isolation, the selected region (V3-V4) of the species-specific 16S rRNA gene was amplified by a polymerase chain reaction (PCR), using a universal set of Illumina adapter containing primers: 341F (5 ′ -TCGTCGGCAGCGTCAGAT GTGTATAAGAGACAGCCTACGGGNGGCWGCAG-3 ′ ) and 805R (5 ′ -GTCTCGTGGGCTCGGAGATGTGTATAAGAGA CAGGACTACHVGGGTATCTAATCC-3 ′ ). V3-V4 region was chosen to also assess differences in Gardernella abundance between groups (32). For the preparation of libraries, we used the standardized 16S Metagenomic Sequencing Library Preparation Protocol (Illumina R ), using KAPA HiFi HotStart ReadyMix 2X for amplification. The amplicons of each sample were labeled with Nextera XT Indexes. The PCR products were quantified and quality checked by using Bioanalyzer DNA 1000 chip (Agilent Technologies, Santa Clara, CA, USA). Libraries of pooled samples were sequenced on the Illumina MiSeq sequencer according to the manufacturer's specifications in the 2 × 300 bp pair-end runs (MiSeq Reagent kit v3).

Bioinformatics Analysis
The sequence data were imported into the R statistical environment (version 3.5.1) from demultiplexed FASTQ files, containing sequences with an average Phred score >30. Microbiome analysis was performed using Bioconductor software 1 following Callahan et al. (33) workflow for Big Data. Trimming and filtering were carried out with The Divisive Amplicon Denoising Algorithm 2 (DADA2) (34) on paired reads jointly, with forward reads truncated at position 260 nt, and the reverse reads at position 240 nt based on the quality scores. First 17 nucleotides of forward reads and 21 nucleotides of reverse reads were trimmed to remove primers with ambiguous sequences. Filtering parameters were as follows; zero ambiguous base calls; enforcement of a maximum of 2/3 expected errors per forward/reverse reads, respectively. The relationship between quality scores and error rates was assessed separately for each run in order to minimize batch effects resulting from run-to-run variability. DADA2 was used to produce a sequence table by independently inferring amplicon sequence variants (ASVs) from the forward and reverse sequences of each sample using the run-specific error rates and to further merge read pairs. ASV represents a higher-resolution analog of the common operational taxonomic unit (OTU) table (34), such that it clusters the sequence reads into a finer taxonomic resolution without the need for a particular similarity cut-off (34,35). Sequences identified as chimeric were removed. Taxonomy was assigned to ribosomal sequence variants using the RDP (Ribosomal Database Project) naive Bayesian classifier algorithm in the DADA2 R package and Silva v128 database (36). Only annotations that had more than 0.6 in their bootstrap confidence estimation values were accepted for further analyses. Clinical, taxonomic, AVS and phylogenetic data were combined into a single object using the phyloseq package (37) (version 1.24.2) for R. Due to the limitations of species-level classification with 16S rRNA sequencing, an additional aligning of representative sequences was performed by the Basic Local Alignment Search Tool (BLAST) (38), where necessary.

Statistical Analysis
All statistical analyses were conducted using R programing language (version 3.5.1). Alpha-diversity metrics were calculated at ASV-level, after rarefying at the minimal sequencing depth of 39,998 reads by performing 1,000 iterations of random subsampling. Rarefaction analysis was performed to evaluate whether further sequencing would likely detect additional AVSs, using the vegan 2 package for R. The microbial richness (number of unique ASVs) was evaluated by the Chao1 index, and alpha diversity by Shannon diversity index, which takes into account richness as well as the abundance of ASVs for each sample, using the phyloseq package (version 1.24.2) (37). Continuous variables were tested for the normality with the Shapiro-Wilk normality test and visually inspected by histograms. Differences in richness and diversity between the different sample groups were calculated by the nonparametric Wilcoxon rank-sum test and the Kruskal-Wallis test when comparing three groups (39).
To estimate differences in beta-diversity, ASVs occurring in fewer than 3 samples and with a count of <10 across all samples were excluded from raw data. The ASV table normalized to relative abundance values and multidimensional scaling (MDS) was performed on phylogenetic distance matrices; weighted UniFrac (40). Variation in the overall microbial community structure between groups was assessed with permutational multivariate analyses of variance (PERMANOVA), using adonis function implemented in the R package vegan with 999 permutations for significance testing. Differences between the outcome groups were visually presented by generating Principal Coordinates Analysis (PCoA) plots and confidence ellipses in the R package ggplot.
The differential relative abundance of taxa between groups was determined at ASVs level using the Bioconductor package DESeq2 (version 1.20.0) in R (41). After variancestabilizing transformations, we normalized samples using the nonparametric Wald negative binomial regression. Results were expressed as log2 fold change and taxa were considered significantly differentially abundant between groups if their adjusted P-value was <0.05. ASVs with prevalence >3% and >50 reads study-wide were included in the DESeq2 analysis.
Vaginal community state types (CST) were determined at ASV level, following the script by DiGiulio et al. (24). The Bray-Curtis dissimilarity metric, based on relative abundance was used to calculate the distance between all samples. Medoid clustering approach (pam in R) was applied to PCoA distances. The appropriate number of clusters was determined using the gap statistic (k = 5) (42). For hypothesis testing of specific-taxa associations, we performed one-sided Wilcox rank-sum tests on the mean relative abundance of a taxon across samples from the same group. P-values were adjusted for multiple comparisons using the Benjamini-Hochberg false-discovery rate (FDR) correction method (43).

Participants Characteristics
Vaginal microbiome sequencing was performed in 158 Caucasian pregnant women who contributed vaginal swabs at the admission to the perinatal delivery department. Three cases were excluded after the medical records were reviewed because the PTD was indicated due to IUGR. The majority of the participants had a vaginal delivery (94.8%), except eight controls, who delivered with a planned C-section (5.2%). Table 1 reports demographic and clinical features of participants. Of the included 155 women, 48 (31.0%) subsequently experienced PTD and 34 (70.8%) of these PTDs were classified as EPTD (<34 0/7 weeks.). Women who delivered at term did not differ significantly from those who delivered preterm in terms of age, working status, parity, pre-pregnancy BMI, the presence of thyroid disease, gestational diabetes mellitus (GDM), presence of asymptomatic bacteriuria, presence of vaginitis and vaginal treatment during pregnancy, urogenital infection, and smoking status ( Table 1). All of the included participants were Caucasians of Slovenian ethnicity.
As reported previously, the women with prior PTD history were more likely to experience subsequent PTD (P = 0.0476). Of the PTD group, 18 women (37.5%) experienced PPROM. As expected, the PTD group had a greater number of women, receiving progesterone therapy (P = 0.00177). None of the women reported vaginal intercourse during the 3 days before sampling.

16S Sequencing Results
A total of 155 swab samples were analyzed, comprising of 48 PTD samples and 107 controls. After initial quality control steps and chimera removal, Illumina MiSeq sequencing of 16S rRNA gene resulted in a total of 12,233,991 processed paired-end V3-V4 reads, with an average read count of 78,928.97 reads per sample (min = 39,998; max = 223,503) (Supplementary Figure S1A).
Next, rarefaction curves were generated by subsampling evenly to the minimal depth of 39,998 reads. Rarefaction analysis showed that curves for each sample reached the plateau and thus no additional ASVs were expected to be found by increasing sequencing depth (Supplementary Figure S1B).

Alpha and Beta Diversity
We demonstrated that the microbiome of participants across EPTD, LPTD, and term gestational delivery groups differed significantly in richness ( Figure 1A) (Chao1, Distributions of categorical variables are compared using Pearson chi-square test at a 5% level of significance and distribution of continuous variables with non-normal distribution with Wilcoxon rank-sum test. Characteristic with significant differences (at 5% level) between preterm and term cohorts are in bold.
Kruskal-Wallis, P = 0.034). Individual comparison using Wilcoxon rank sum test showed that Chao1 index decreased in the term group compared to the all preterm cases (Chao1, P = 0.011) (Supplementary Figure S2) or EPTD group (Chao1, P = 0.012), indicating the reduced richness in women who delivered term babies.
Comparison of bacterial diversity between EPTD, LPTD, and term, revealed significant differences by Shannon index in the Kruskal-Wallis test (P = 0.0021) (Figure 1A). By Wilcoxon rank sum test we demonstrated a significant decrease in diversity in the term group compared to all preterm cases (Shannon, P = 0.00059) (Supplementary Figure S2A) or EPTD group (Shannon, P = 0.00079).
Next, we independently tested PPROM for alpha diversity, comparing cases who experienced PPROM (n = 18) to cases who delivered spontaneously without PPROM (n = 30). Among mothers with PTD, those who experienced PPROM showed a trend toward greater richness when compared to spontaneous PTD cases without PPROM, although the results were not statistically significant (Wilcoxon test, Chao1, P= 0.097) (Supplementary Figure S3A).
With respect to beta diversity, weighted UniFrac analysis demonstrated that overall microbial community composition of

Composition of the Vaginal Microbiome
In women, who delivered at term, the Lactobacillus genus constituted 88.3% and Gardnerella 4.63% of total bacteria. In that group 91/107 (85%) samples had Lactobacillus abundance above 70%, and 55/107 (51%) samples harbored Lactobacillus genera above 98%. In the EPTD group, Lactobacillus comprised 50.3% and Gardnerella 17.6%, and in LPTD 75.5 and 5.22%, respectively. Atopobium was present at higher percentages in the EPTD (5.05%) and LPTD (2.13%) groups when compared to the term group (1.59%). Relative abundances of genera within each sample according to a week of delivery are presented in Supplementary Figure S4.

Relative Abundance Comparisons of Individual Taxa
At the genus level, we analyzed and confirmed the reduced mean relative abundance of Lactobacillus in the PTD group when compared to term (P adj = 0.000177), using one sided-Wilcoxon rank-sum test on taxonomically grouped data.
Furthermore, we tested for some of the additional previously reported associations between PTD and genera observed to have increased abundance during BV or AV (24,25) (Table 3). At the genus level, significant associations for following taxa, previously associated with BV were confirmed in our dataset: Gardnerella, Atopobium, Sneathia, Gemella, Megasphaera, and Dorea. We did not detect significant associations with PTD in mean relative abundances of the following genera: Ureaplasma, Mycoplasma, Mobiluncus, and Peptoniphilus. Among genera associated with AV, we detected an increased abundance of Streptococcus and Escherichia/Shigella ( Table 3).
At the species level, we individually tested for L. crispatus and L. iners associations with PTD. The mean relative abundance of ASV belonging to Lactobacillus acidophilus/casei/crispatus/gallinarum (classified as 100% L. crispatus/L. gallinarum according to BLAST) was enriched in the term group compared to PTD (one-sided Wilcox-test, P adj = 0.0201), with 38.1% in the term and 26.9% in the PTD group. Association of L. iners with gestational age delivery group was not significant.

Vaginal Microbiome Community State Types (CSTs)
Hierarchical clustering analysis demonstrated that vaginal microbiota clustered into 5 distinct community profiles, consisting of L. crispatus (CST I), L. gasseri (CST II), L. iners (CST III), and L. jensenii/L. iners -dominated clusters (CST V) and a cluster depleted of Lactobacilli species and enriched for anaerobic bacteria (CST IV) (Figure 3). Frequencies of specific CSTs in the entire sample were as follows: CST I (31.6%), CST II (11.6%), CST III (32.3%), CST IV (16.8%), and CST V (7.74%). In CST IV cluster were 41.2% of EPTD participants, 21.4% of LPTD participants, and only 8.41% of the term participants ( Figure 3A). CST I (dominated by L. crispaus) was present in 38.3% of term delivery group, 17.6% in EPTD and 14.3% in the LPTD group ( Figure 3A). We observed a significant difference in the overall frequency of five CSTs between three delivery groups (term, LPTD, and EPTD), P Fisher exact test = 0.0025, which was even more profound when only term and the EPTD groups were considered (P Fisher exact test = 0.0005).
When we tested previously reported association of Lactobacillus-depleted CST IV with PTD, the CST IV was significantly overrepresented in the EPTD group compared to term (P Fisher exact test = 0.0008).

DISCUSSION
Our findings showed that pregnant women who delivered preterm relative to term have higher richness and diversity of the vaginal microbiome, along with reduced Lactobacillus species, especially reduced L. crispatus and the increased abundance of the genus Gardnerella. Additionally, we found that overall frequencies of CSTs were significantly different between EPTD, LPTD, and term delivery groups, with CST IV (enriched for BVassociated taxa and higher diversity) observed in the EPTD group more often than in the term group. These observations are in accordance with the majority of the previously published 16S sequencing studies investigating the similar research question in Caucasian populations. The higher alpha diversity and lower abundances of Lactobacillus spp. (CST IV) in the early stages of pregnancy were previously correlated with gestational age at delivery (24,25). CSTs were previously shown to differ between term and preterm deliveries also in the study of Stafford et al. (27), where PTD has been associated with CST V (L. jesenii) dominance, while CST I (L. crispatus) and CST II (L. gasseri) were beneficial for term delivery. The association between PTD and vaginal dysbiosis (CST IV) was not confirmed in the study of Kindinger et al. (22). This discrepancy may be explained by the inclusion of only women at high-risk for PTD, possibly making the microbialcause for PTD less prominent in that cohort. The study of Freitas et al. (44), conducted in Canadian cohort also failed to detect specific vaginal microbiome CSTs early in pregnancy that could predict PTD, however, the findings of increased richness and diversity in PTD are in agreement with our study, despite different methodologies used (pyrosequencing of cpn60).
While the protective role of L. crispatus was repeatedly demonstrated in most previous 16S sequencing studies (21,22,25,27,29), association of L. iners with PTD yielded conflicting results. The largest 16S sequencing study of white European population investigating the role of vaginal microbiota in PTD (21), concluded that L. iners/Ralstonia solanacearum oligotype was associated with decreased risk of early (<34 weeks), but not late (34-36 weeks) PTD, as well as L. gasseri/L. johnsonii, L. crispatus/L. acidophilus, and Bifidobacterium longum/Bifidobacterium breve oligotypes. In contrast, a recent 16S sequencing study of V1-V3 region by Kindinger et al. (22) showed that dominance of L. iners at 16 weeks of gestation is a potential risk factor for early PTD (<34 weeks). We did not manage to confirm the association with L. iners with either term or preterm delivery, however, our results indicated that its mean relative abundance was the highest in the LPTD group (51.8%), while it was more prevalent in term group (32.6%) compared to EPTD (21. 9%).
Consistent with the previous reports of an association between BV and AV with PTD (24,25,29), we detected increased abundance of BV-associated taxa, including Gardnerella, Atopobium, Sneathia (Sneathia sanguinegens), Gemella, Dorea, and AV-associated bacteria, Streptococcus (Streptococcus agalactiae/pyogenes), and Escherichia/Shigella genera in the PTD group, when compared to term. Contrary to the study of Callahan et al. (25), we did not confirm the association of genera Mobiluncus, Peptoniphilus, and Mycoplasma with PTD; however, we did detect the positive association of Megasphaera with PTD, which was not confirmed in a before-mentioned study. The association of Ureaplasma with PTD (24) was not replicated in our cohort, nor was it in the studies of Callahan et al. (25) and Freitas et al. (44). Similarly to our study, Elovitz et al. (7) reported the significant association of Atopobium vaginae with increased rates of PTD in non-Afro-American (AA) women, while Sneathia, Atopobium, and Megasphaera were associated with PTD when all women (AA and non-AA) were considered.
Our results differ from the results of the studies by Stout et al. (26) and Romero et al. (23), where no specific bacterial taxa or community was correlated with pregnancy outcome. Nevertheless, indicated PTDs (preeclampsia, fetal indications) were not excluded from the first study and only taxa that were present in at least 25% of samples were considered in the differential abundance analysis in the study of Romero et al. (23). We suppose that geographical and cultural differences (predominantly AA population), as well as a low number of PTD cases included (n = 23 and n = 18, respectively), might also contribute to negative results of those two studies.
Decreased Lactobacilli abundance manifests in the reduced lactic acid, hydrogen peroxide, and bacteriocins production, which enhances the proliferation of pathogenic bacteria known to cause BV (45,46). Additionally, increased bacterial diversity has previously been associated with mucus and cytoskeleton alterations, cell death, increased expression of pro-inflammatory cytokines (29) and metalloproteinases and decreased IgG1/2 (47), which potentially leads to the inflammatory activation and stimulation of PTD. Furthermore, in the recent study, women with a non-Lactobacillus spp.-dominant microbiome were more likely to develop premature cervical dilation (48). It was also demonstrated (28,49) that PPROM, which precedes 30% of all PTDs, is associated with high bacterial diversity and reduced Lactobacillus abundance prior to rupture. In the present study, we detected a tendency toward increased richness when preterm cases with PPROM were compared with preterm cases who delivered spontaneously without PPROM, but the result did not reach statistical significance.
It is already known that 16S sequencing experimental design impacts sensitivity to different taxa, and here we used the V3-V4 primer set, which is also sensitive to Gardnerella. Consequently, in our dataset Gardnerella represented 7.7% of the genus relative abundance, which is similar to study of the Callahan et al. (25) (5% of reads). The higher relative abundance of L. iners in the PTD group in the study of Kindinger et al. (22) may, therefore, be due to undetected Gardernella (V1-V3 primer set), which was shown to often coexist at high frequencies with L. iners (25), rather than a direct association of L. inners with PTD. Another explanation is that increased L. iners abundance during the 16th week of gestation might be associated with increased bacterial instability later in pregnancy, as L. iners is often associated with BV-related vaginal bacteria, due to its possible lack of protection against pathogens (50).
In contrast, L. crispatus has the largest genome of four Lactobacilli spp. (45), including genes encoding for bacteriocin, adhesin (anti-adhesive molecules against Gardnerella vaginalis), lactate, and succinate (51) and was reported to have protection against inflammation induced PTD (51). It has also been demonstrated that D-lactic acid levels were significantly higher when L. crispatus was the dominant vaginal bacterium compared to L. iners or Gardernella dominated vaginal microbiota. It was demonstrated that D-lactic acid enhances activation of anti-microbial innate and acquired immune system and also minimizes matrix metalloproteinase-8 (MMP-8) levels that would otherwise promote proteolytic activity and thus enhance bacterial migration past endocervix into the uterus, leading to infection-related PTD (52).
The strength of the present study is the inclusion of a large number of early PTDs (<34 weeks) compared to previous reports. Furthermore, we utilized strict inclusion criteria, meaning that PTDs indicated due to maternal or fetal complications, such as IUGR, preeclampsia, or fetal anomalies were excluded from the study. With an aim to exclude potential pre-existing risk factors for PTD and factors that could independently influence on microbiota composition, we also excluded individuals with chronic diseases, such as hypertension, diabetes mellitus, or women with the uterine anomalies. Another strength is the homogeneity of the study population, while it is known that microbiota composition is largely influenced by ethnicity (7,17,30). All participants were Caucasians and belonged to the same Slovenian ethnicity. Furthermore, compared to the most studies, we utilized bioinformatic approach with improved resolution and comprehensiveness that resolves exact ASVs, instead of imposing the arbitrary dissimilarity thresholds that define molecular OTUs (35).
Our study also has some limitations that need to be addressed. Firstly, sampling was conducted at a single time point; longitudinal sampling through pregnancy could add information about microbiome stability. Although previous studies showed that vaginal microbiome in the pregnancy remains relatively stable (24,53), the increased Lactobacillus (23) and significant decrease in abnormal vaginal colonization rate were reported as the pregnancy progressed (54). Another limitation of our study is that 16S-sequencing technology does not assess the total abundance load of the taxa in the sample (gene copy measurement) that could be detected using a broad-range 16S rRNA quantitative PCR. Immune factors, that were suggested to lower the risk of PTD associated with vaginal microbiota, such as beta-defensin-2 (7) or proinflammatory cytokines (29), were not measured.
With a more comprehensive understanding of the vaginal microbiota in pregnancy in different populations, and its effect on PTD, detection of mothers at high-risk and development of pharmacological interventions for PTD prevention will become plausible.
In conclusion, our results suggest that maternal vaginal microbiome has an important role in the PTD etiology. In our independent analysis, we confirmed that vaginal microbiome of women who delivered preterm has increased richness and diversity prior delivery, decreased Lactobacillus abundance and the increased abundance of BV-associated taxa. Nevertheless, the role of vaginal microbiota in the etiology of PTD remains to be elucidated in further research considering different ethnicities and ideally with the integration of hosts and fetus genetic data.

DATA AVAILABILITY
Generated Statement: The datasets generated for this study can be found in the Sequence Read Archive (NCBI) repository, BioProject PRJNA544732, http://www.ncbi.nlm.nih.gov/ bioproject/544732.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Slovenia National Medical Ethics Committee, Ministry of Health, Štefanova 5, SI-1000 Ljubljana, Slovenia (approval number 90/02/15). The patients/participants provided their written informed consent to participate in this study.