Association between Pregestational Vaginal Dysbiosis and Incident Hypertensive Disorders of Pregnancy Risk: a Nested Case-Control Study

ABSTRACT A balanced vaginal microbiome dominated by Lactobacillus can help promote women’s reproductive health, with Lactobacillus crispatus showing the most beneficial effect. However, the potential role of vaginal microbiomes in hypertensive disorders of pregnancy (HDP) development is not thoroughly explored. In this nested case-control study based on an assisted reproductive technology follow-up cohort, we prospectively assessed the association between pregestational vaginal microbiomes with HDP by collecting vaginal swabs from 75 HDP cases (HDP group) and 150 controls (NP group) and using 16S amplicon sequencing for bacterial identification. The vaginal microbial composition of the HDP group significantly differed from that of the NP group. The abundance of L. crispatus was significantly lower, and the abundances of Gardnerella vaginalis was significantly higher, in the HDP group than in the NP group. Of note, L. crispatus-dominated vaginal community state type was associated with a decreased risk for HDP (odds ratio = 0.436; 95% confidence interval, 0.229 to 0.831) compared with others. Additionally, network analysis revealed different bacterial interactions with 61 and 57 exclusive edges in the NP and HDP groups, respectively. Compared with the HDP group, the NP group showed a higher weighted degree and closeness centrality. Several taxa, including G. vaginalis, L. iners, and bacterial vaginosis-associated bacteria (Prevotella, Megasphaera, Finegoldia, and Porphyromonas), were identified as “drivers” for network rewiring. Notable alterations of predicted pathways involved in amino acid, cofactor, and vitamin metabolism; membrane transport; and bacterial toxins were observed in the HDP group. IMPORTANCE The etiology of HDP remains unclear to date. Effective methods for the individualized prediction and prevention are lacking. Pregestational vaginal dysbiosis precedes the diagnosis of HDP, providing a novel perspective on the etiology of HDP. Early pregnancy is the critical period of placental development, and abnormal placentation initiates HDP development. Thus, disease prevention should be considered before pregnancy. Vaginal microbiome characterization and probiotic interventions before pregnancy are preferred because of their safety and potential for early prevention. This study is the first to prospectively assess associations between pregestational vaginal microbiome and HDP. L. crispatus-dominated vaginal community state type is linked to a reduced risk for HDP. These findings suggest that vaginal microbiome characterization may help identify individuals at high risk for HDP and offer potential targets for the development of novel pregestational intervention methods.


RESULTS
Participant characteristics. This study included 225 individuals with normal baseline blood pressure, including 75 women who developed HDP and 150 matched controls who delivered at term without any complications after embryo transfer (Fig. 1). The characteristics of the participants are detailed in Table 1. Risk factors for HDP, including maternal age, body mass index (BMI), baseline blood pressure, nullipara, indicators for glucose and lipid metabolism, number of fetuses, and incidence of polycystic ovary syndrome (PCOS), revealed no difference between the HDP and NP groups. Of the 75 women in the HDP group, one died of sudden elevated blood pressure at 35 weeks of gestation, and another had a late miscarriage at 26 weeks (Table S1 in the supplemental material).
Pregestational vaginal microbiome changes in women who developed HDP. A total of 986 amplicon sequence variants (ASVs) were obtained in the two groups. The rarefaction curve of the ASV number revealed that the sequencing depth was sufficient to capture most of the microbial characteristics (Fig. S1). Moreover, alpha diversities indicated by Shannon, Richness, Chao1, and Simpson diversity indexes were comparable between the two groups ( Fig. 2A, Fig. S2, conditional logistic regression, P . 0.05). The PCoA based on the weighted Unifrac distance indicated that the vaginal microbial composition of the HDP group significantly differed from that of the NP group (Fig. 2B, PERMANOVA, R 2 = 0.014, P = 0.037). This result was further validated by the Bray-Curtis and Jaccard distances ( Fig. S3A and B). Furthermore, the dissimilarities among individuals within the HDP group were significantly larger than those among individuals within the NP group (Fig. 2C, Fig. S3C and D, Wilcox's rank-sum test, P , 0.001) and were comparable to the intergroup distances, indicating that the vaginal microbiome was highly divergent among individuals in the HDP group.
To further illuminate the difference in microbial community compositions between the two groups, we compared the relative abundance of the dominant bacterial taxa with mean relative abundances over 0.1%. At the phylum level, the two groups showed similar compositional patterns, both dominated by Firmicutes, Actinobacteria, Bacteroidetes, Proteobacteria, Tenericutes, and Fusobacteria (Fig. S4A). The top five abundant genus in the two groups were Lactobacillus, Gardnerella, Atopobium, Prevotella, and Streptococcus (Fig. S4B). The abundance of Lactobacillus crispatus was significantly lower (q = 0.046), while the abundance of Gardnerella vaginalis was higher (q = 0.046) in the HDP group than in the NP group (Fig. 2D).
Associations of identified differential bacterial taxa with a gestational weeks at delivery and neonatal birth weight were examined using Pearson's correlation analysis (Fig. S5A). L. crispatus was positively associated with gestational duration, whereas women with a high abundance of G. vaginalis were vulnerable to an early delivery (P = 0.038, P = 0.005, respectively). Noticeably, G. vaginalis was significantly linked with a low neonatal birth weight (P = 0.012), and a similar correlation of G. vaginalis can be found in twin and singleton pregnancies (P = 0.031). Vaginal CSTs associated with HDP development. The vaginal microbial communities were heterogeneous and diverse in the population. To further explore whether the vaginal community structure diverged between the two groups and to stratify populations and simplify the analysis of the vaginal microbiome, we assigned vaginal microbial profiles to community state types (CSTs) according to the most abundant taxon (Fig. 2E,  Fig. S5B, Table S2). The L. crispatusand L. iners-dominated CSTs accounted for 46.2% and 41.3% of the total vaginal samples, respectively (Table S2). The L. crispatus-dominated CSTs constituted more than half (52.0%) of the samples in the NP group but only 34.7% of the samples in the HDP group, which coincided with the higher L. crispatus abundance in the NP group than in the HDP group (Fig. 2E, Fig. S5B Table S3). Nevertheless, the women with L. iners predominance were prone to HDP (OR = 1.732; 95% CI, 0.947 to 3.165, P = 0.074; Table S3).
Alterations in bacterial interactions of vaginal microbiota in the HDP group. We investigated whether bacterial interactions, a critical part of the ecosystem, were altered in the women who developed HDP. Cooccurrence networks were constructed to analyze the relationships between bacterial taxa in the two groups (Fig. 3A). To assess the concordance of microbial cooccurrence between the HDP and NP groups, we first counted the number of exclusive and shared nodes and edges. Although more than half the nodes (33/57) were shared, 61 and 57 unique edges were found in the NP and HDP groups, respectively, with only 42 shared edges ( Fig. 3B and C). Differences in the topological structure of networks were also found between the two groups. Compared with those in the HDP group, the bacterial interactions in the NP group were denser, with a higher weighted degree and closeness centrality (Fig. S6, Wilcox's rank-sum test, P , 0.001, P = 0.008, respectively). A similar pattern was also observed in the shared nodes (Fig. 3A, Fig. S6, Wilcox's rank-sum test, P , 0.001, P = 0.044). The comparable degree between the two networks indicated that the total connections of nodes in the two networks were similar (Fig. S6).
Finally, we inferred "driver" taxa for network rewiring from NP to HDP by analyzing the most common subnetwork. The following 10 taxa essentially contributing to the remodeling of the entire microbial community structure were identified as potential "driver" taxa (22): G. vaginalis, L. iners, Ureaplasma parvum, Porphyromonas uenonis, Finegoldia magna, Lactobacillus vaginalis, Prevotella disiens, Prevotella buccalis, Gordonibacter_unclassified, and Megasphaera lornae (Fig. 3A). Overall, the distinct components and topology of the networks indicated that the discrepancy in vaginal ecology between the two groups was driven not only by differential abundant taxa but also by the different interactions of shared microbes.
Functional characteristics of vaginal microbiota in normotensive women who developed HDP. The function profiles were inferred to determine the functional capacity of the vaginal microbiota (Fig. 3D, conditional logistic regression, all P , 0.05 with q , 0.2). A comparative analysis of the predicted Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways showed that the HDP group was enriched in pathways associated with environmental information processing and membrane transport, such as bacterial toxins, transporters, starch and sucrose metabolism, peptidases, methane metabolism, and taurine and hypotaurine metabolism. In contrast, the vaginal microbiota of the NP group was mainly enriched in pathways involved in the metabolism of amino acids, cofactors, and vitamins, including histidine, tyrosine, arginine, proline, and riboflavin, and biosynthesis of folate.

DISCUSSION
The etiology of HDP remains unclear to date. Consequently, effective methods for the individualized prediction, prevention, and treatment of this disease are lacking. The role of the microbiota in HDP development has attracted increasing interest (10 to 14). To the best of our knowledge, this study is the first to prospectively assess associations between pregestational vaginal microbiota and HDP. Study participants were from a large cohort with high follow-up rates and detailed medical information before and during pregnancy, which enabled us to exclude women with pregestational diabetes mellitus (PGDM) or gestational diabetes mellitus (GDM), considering they have increased risks of HDP and vaginal the NP and HDP groups based on Unifrac distance. The vaginal bacterial community composition is significantly different between the HDP and NP groups. Each dot represents a sample and is colored by group. The difference in beta diversity between the NP and HDP groups was determined using PERMANOVA. (C) Unifrac distance of the vaginal microbial community among individuals within each group and between different groups. Significance in comparison of Unifrac distance was determined using a two-tailed Wilcox rank-sum test. (D) Differentially abundant bacterial taxa between the two groups. Significance was determined using conditional logistic regression corrected for multiple testing with significance identified at q , 0.  (3,23). The good match between groups on confounders, including age, BMI, baseline blood pressure, and nulliparity, made our observed associations reliable. This study showed that pregestational vaginal dysbiosis precedes the diagnosis of HDP in normotensive women receiving in vitro fertilization or intracytoplasmic sperm injection and embryo transfer (IVF/ICSI-ET) treatment, who were considered to have a relatively high risk of HDP (3,4,6). Of note, L. crispatus-dominated vaginal microbiome was associated with a decreased likelihood of developing HDP by 0.436. As demonstrated in the present study, the abundance of L. crispatus was lower, while the abundance of G. vaginalis was higher, in the women who developed HDP than in the healthy controls. L. crispatus enhances the trophoblastic invasion by upregulating the expression of matrix metalloproteinases (MMPs) (24), which might explain its beneficial effects on placental implantation and pregnancy outcome (25 to 28). MMPs are important regulators in vasodilation, placentation, and uterine expansion during normal pregnancy (29). The protein expression of MMP-1 in extravillous trophoblast is decreased in preeclampsia (30). Reduced MMP-2 expression, decreased uterine vascularization, and arterial expansive remodeling have been observed in a previous model of hypertensive pregnancy, possibly contributing to uteroplacental ischemia (31). Depletion of L. crispatus may be related to the impaired invasion of the trophoblast, which may result in shallow placentation and inadequate remodeling of the spiral artery, leading to placental ischemia and preeclampsia (4). Furthermore, high L. crispatus abundance was positively associated with gestational age at delivery, whereas the opposite association was found for G. vaginalis. This result is reasonable because HDP was closely related to medically indicated preterm birth and fetal growth restriction, which also indirectly reflect the severity of the disease.
A decreased abundance of L. crispatus and a higher level of G. vaginalis in the female reproductive tract were also associated with spontaneous preterm birth (32 to 36), early missed abortion (37), recurrent pregnancy loss (38), and lower fecundability (39). The overlap of bacterial associations between adverse pregnancy outcomes and HDP suggests that maternal diseases might have a common vaginal dysbiosis that resembles the common gut dysbiosis in different diseases (40,41). Further studies are needed to depict detailed mechanisms underlying the overlap. The overlapping mechanisms could be partially mediated by the inflammation triggered by these pathogens because proinflammatory cytokines and complement proteins reportedly induce abnormal placentation (42 to 44).
Despite its lack of differential abundance, L. iners was identified as the "driver," essentially contributing to the remodeling of the entire microbial community from NP to HDP together with G. vaginalis. Vaginal microbiomes dominated by L. iners transform more often into a diverse community than L. crispatus-dominated CSTs (45,46). L. crispatus and G. vaginalis tend to be exclusive in the vagina, whereas L. iners often coexists with G. vaginalis (35). L. iners enhances the adhesion of Gardnerella to epithelial cells, and Gardnerella displaces L. crispatus but not L. iners from the epithelia (47). The decreased proportion of the L. crispatus-dominated CST and the increased proportion of the L. iners-dominated CST, and the latter's propensity to shift toward diverse communities, could partly explain the larger disparity among individuals in the HDP group than among those in the NP group.
Additionally, notable alterations in bacterial interactions, which are important in the microecosystem, were found in the women susceptible to HDP. Despite shared nodes, a large number of edges were exclusive to the women in the HDP or NP group, indicating the remodeling of the entire microbial structure. Hence, the discrepancy between the two groups was not only driven by differentially abundant bacteria but also by the distinct bac- FIG 3 Legend (Continued) are labeled. The size of nodes represents the mean relative abundance of the taxa. The edge width reflects the strength of correlations. The edge color denotes the direction of correlations, with red and blue representing positive and negative correlations, respectively. Only significant bacterial connections beyond the cutoff mean 6 three standard deviations are retained. The heatmaps shown at the top depict comparisons of weighted degree and closeness centrality of shared nodes in the two networks using a two-tailed Wilcox rank-sum test. The heatmaps are colored by rank of weighted degree or closeness centrality, with red denoting a high rank. *, P , 0.05; **, P , 0.01; ***, P , 0.001. (B) Exclusive and shared nodes (B) as well as edges (C) in the two networks. (D) Predicted KEGG pathways significantly different between the two groups. The left part depicts the mean relative abundance of differential predicted pathways (q , 0.2) between the NP and HDP groups. The right part depicts differences in abundance between the two groups with 95% CI. Red denotes pathways enriched in the HDP group, while blue denotes pathways enriched in the NP group. terial interactions of shared microbes. Several taxa were identified as potential "drivers" during the remodeling, of which the majority belong to bacterial vaginosis-associated bacteria, namely, G. vaginalis, Porphyromonas, Finegoldia, Prevotella, and Megasphaera (48,49). Gardnerella, also identified as a differential taxon, can initiate biofilm formation, which then serves as a scaffold for other species to adhere to during vaginal dysbiosis (50,51). Prevotella, Porphyromonas, Streptococcus, Ureaplasma, and L. iners were detected in preeclampsia amniotic fluid and placenta by cultivation and sequence-based methods, supporting the involvement of dysbiosis in HDP development (15,16).
In the present study, the NP group was enriched with pathways related to the metabolism of histidine, tyrosine, arginine, riboflavin, and folate, whereas the HDP group was enriched with pathways related to transporters and bacterial toxins as signs of bacterial virulence and predisposing factors of inflammation. High dietary histidine and tyrosine intake are associated with low blood pressure (52,53). Decreased serum L-arginine, which plays a key role in nitric oxide metabolism, is a risk factor for HDP (54). Riboflavin, a methylenetetrahydrofolate reductase (MTHFR) cofactor for folate metabolism, could lower blood pressure by restoring MTHFR activity in vascular cells, improving nitric oxide availability and endothelial function (55,56). Dietary deficiencies of folate and riboflavin are related to increased plasma homocysteine, reduced DNA methylation patterns, and genome instability events, which increase cardiovascular risk (57). These studies suggest that the microbiome participates in HDP development by modulating metabolism and inducing inflammation. A recent study has similarly demonstrated that reductions in short-chain fatty acid-producing bacteria accompanied by decreased propionate and butyrate contribute to HDP development by causing inflammation and impaired trophoblast invasion (14).
Several limitations of the present study need to be acknowledged. Although this study was based on a large-scale population, only 75 cases were included in the final analysis because of the low incidence of HDP (4.9%, 108/2224) and the exclusion of cases with other complications. Moreover, the study participants were infertile women of the Chinese Han population who received IVF/ICSI treatment. Caution should be taken when generalizing the results of this study to other populations in other geographical regions. Large-scale and longitudinal multicenter studies and animal experiments are therefore needed to verify these findings.
Conclusions. In conclusion, pregestational vaginal dysbiosis in normotensive women precedes the diagnosis of HDP. L. crispatus-dominated vaginal microbiomes are associated with a reduced risk for HDP. Our study extends the understanding of potential links between vaginal perturbations and pregnancy complications, provides novel insights into the etiology of HDP, and sheds light on a novel path for the prevention of this disease.

MATERIALS AND METHODS
Study participants. This study was derived from a comprehensive microbiome cohort established in 2019 with vaginal swabs collected from Chinese females aged 19 to 64 years at routine clinic visits at the Center for Reproductive Medicine, Shandong University.
The exclusion criteria were as follows: antibiotic treatment within 30 days, vaginal medications or vaginal douching within 7 days, engagement in sexual activity in the current menstrual cycle, vaginal bleeding, vaginal subjective symptoms (e.g., vaginal itching and abnormal discharge), any acute inflammation, and severe systemic disease. HDP is defined as new-onset hypertension during pregnancy (gestational hypertension or preeclampsia) after 20 weeks of gestation.
The workflow of the study is as follows (Fig. 1). A total of 5,849 normotensive women scheduled for IVF/ ICSI-ET were recruited. Among them, 4,345 individuals received fresh or frozen embryo transfer. Then, 2,224 achieved persistent pregnancy after 12 weeks of gestation, and one was lost to follow-up after 16 weeks of gestation. Finally, 108 women were diagnosed with HDP, of which 75 women who were not complicated with GDM or PGDM were included as cases (here denoted as "HDP group"). Meanwhile, 150 uncomplicated controls (here denoted as "NP group") with term delivery were matched by 1:2 nearest-neighbor propensity score matching (PSM) without replacement with a propensity score estimated using logistic regression models. For the model, the absence or presence of HDP was the outcome variable. The predictor variables included age, BMI, systolic blood pressure, diastolic blood pressure, fasting blood glucose, triglyceride, total cholesterol, high-density lipoprotein, low-density lipoprotein, transfer of fresh or frozen embryos, number of transferred embryos, and number of fetuses. Advanced maternal age, obesity, chronic hypertension, higher fasting blood glucose, multiple gestation, the use of assisted reproductive technologies, frozen embryo transfer, and certain cardiovascular risk factors such as elevated levels of triglyceride, total cholesterol, low-density lipoprotein cholesterol, and low levels of high-density lipoprotein cholesterol were in relation to a higher risk for HDP (3 to 9).

Pregestational Vaginal Dysbiosis and HDP mSphere
These clinical risk factors were selected to estimate the risk for HDP, and further PSM would balance the baseline characteristics between the two groups by reducing selection bias. The PSM method has been widely used in observational studies to reduce selection bias (58 to 60), and the procedure was conducted by the MatchIt package with default parameters (61). Each case, in descending order of the propensity score, was successively matched to the control with the closest propensity score. Two rounds of one-to-one matches were conducted, resulting in a match ratio of 1:2. This process yielded an adequate balance, as evidenced by the standardized mean differences of all variables being within a threshold of 0.15. Ethical approval and consent to participate. The study was conducted in accordance with the Declaration of Helsinki and approved by the Institutional Review Board of Reproductive Medicine, Shandong University (no. 2019LSZ14). All participants provided written informed consent.
DNA extraction from vaginal swabs. Vaginal swabs were collected from the posterior vaginal fornix of the women in the follicular phase of menstrual cycle at routine clinic visits. DNA extraction was performed using Magnetic Soil and Stool DNA kits (Tiangen Biotech, Beijing) following the manufacturer's protocol. The V1-V2 hypervariable regions were amplified. Further details about sample collection, 16S rRNA gene sequencing, and preprocessing of raw data are provided in the supplemental material.
Bioinformatics analysis and statistical analyses. An average of 59,635 reads per sample were obtained after quality filtering. The representative sequences were annotated against the Ribosomal Database Project (RDP) database (rdp_16s_v16_sp.fa accessible at http://www.drive5.com/sintax) (62) using the SINTAX algorithm with a confidence threshold of 0.6. Species allocation was performed by using the RDP and STIRRUPS databases (63).
Sequencing reads were subsampled to the lowest read count of 25,441, followed by rarefaction evaluation. Alpha diversity was quantified using Shannon, Richness, Chao1, and Simpson diversity indices. The weighted Unifrac distance, Bray-Curtis distance, and Jaccard distance were calculated to assess the dissimilarity of microbial composition among intragroup or intergroup samples. The dissimilarity was visualized using principal coordinate analysis (PCoA). Permutational multivariate analysis of variance (PERMANOVA) was conducted using the function adonis from the vegan package to evaluate community structure differences. Discriminatory bacteria were identified by conducting conditional logistic regression on high-abundance taxa with a mean relative abundance over 0.1%. Vaginal profiles were assigned to CSTs based on the most abundant taxon (32). Samples with the largest proportion of less than 30% were not assigned to any CST. Differences in the number of Lactobacillus crispatus-dominated and Lactobacillus iners-dominated CSTs between the two groups were evaluated using conditional logistic regression.
Microbial community networks were constructed as follows. Correlation coefficients were computed at the species level by using Fastspar v2.0 with 1,000 bootstrap replicates to estimate the P values. Bacterial taxa prevalent in at least 5% of the samples were used to filter rare taxa. Correlations with P values of ,0.05 and coefficients beyond the threshold of mean 6 three standard deviations were retained for further network analysis. Networks were visualized using the igraph package. Local property measures, such as degree, weighted degree, and closeness centrality, were calculated. Common and exclusive nodes and edges and the most common subnetwork between two networks were analyzed to decipher the network similarity and the extent of network rewiring. The taxon with an altered set of associations (identified by a high NESH score .1) and increasing importance in the case network (identified with a positive delta betweenness from control to case) was identified as a "driver" by using Netshift (22). Finally, the functions of vaginal microbiota were predicted using PICRUST, and the predicted KEGG pathways with mean relative abundances over 0.15% were compared between the two groups by using conditional logistic regression.
Conditional logistic regression was used to compare the difference between the case and control pairs unless otherwise specified. False discovery rate-corrected P values (q values) were used to define statistical significance at 0.2 for multiple testing.
Data availability. The sequencing data reported in this article have been deposited in the Genome Sequence Archive at the BIG Data Center, Chinese Academy of Science, under accession number CRA009956 (BioProject: PRJCA015148). The data are publicly accessible at http://bigd.big.ac.cn/gsa.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. TEXT S1, DOCX file, 0.02 MB.