Vaginal microbes alter epithelial transcriptome and induce epigenomic modifications providing insight into mechanisms for susceptibility to adverse reproductive outcomes

The cervicovaginal microbiome is highly associated with women’s health, with microbial communities dominated by Lactobacillus species considered optimal. Conversely, a lack of lactobacilli and a high abundance of strict and facultative anaerobes, including Gardnerella vaginalis, have been associated with adverse reproductive outcomes. However, how host-microbial interactions alter specific molecular pathways and impact cervical and vaginal epithelial function remains unclear. Using RNA-sequencing, we characterized the in vitro cervicovaginal epithelial transcriptional response to different vaginal bacteria and their culture supernatants. We showed that G. vaginalis upregulates genes associated with an activated innate immune response. Unexpectedly, G. vaginalis specifically induced inflammasome pathways through activation of NLRP3-mediated increases in caspase-1, IL-1β and cell death, while live L. crispatus had minimal transcriptomic changes on epithelial cells. L. crispatus culture supernatants resulted in a shift in the epigenomic landscape of cervical epithelial cells that was confirmed by ATAC-sequencing showing reduced chromatin accessibility. This study reveals new insights into host-microbe interactions in the lower reproductive tract and suggests potential therapeutic strategies leveraging the vaginal microbiome to improve reproductive health.


Introduction
The female lower reproductive tract is comprised of a complex ecosystem made up of host epithelial and immune cells, microorganisms including the microbiome (bacteria) 1 , mycobiome (fungi) 2 , virome (viruses) 3,4 and their metabolites.The cervicovaginal microbiome has become the focus of many studies due to its highly integrated and complex role in reproductive health and disease.Using high throughput 16S rRNA gene sequencing, the vaginal microbiota composition has been well characterized in both non-pregnant and pregnant individuals [5][6][7] .The vaginal microbiome has traditionally been de ned by the presence or absence of Lactobacillus species 5,8,9 .Cervicovaginal microbial communities dominated by Lactobacillus species are generally considered optimal and associated with positive reproductive health outcomes, while those lacking lactobacilli and comprising a wide array of strict and facultative anaerobes are associated with a diverse spectrum of adverse gynecological and reproductive outcomes including infertility 10,11 , sexually transmitted infections (STIs) (human papilloma virus (HPV) 12 and human immunode ciency virus (HIV)) 13 , and pregnancy complications (preterm birth) 14,15 .
Although the cervicovaginal microbiome is less taxonomically diverse than those found at other body sites, the identi cation of con-speci c genotypes cohabitating in the vaginal microbiome adds to the complexity of these communities 16 .Despite these advances in vaginal microbiome characterization and its association with clinical outcomes, there is a paucity of data regarding the precise mechanisms by which vaginal microbes modulate the function of host epithelia and promote adverse health outcomes.
The vaginal microbiome interacts with the epithelial barriers of the cervicovaginal space.The cervicovaginal epithelial barrier is unique as the epithelial cells that line this space have different embryological origins, resulting in distinct cell-speci c functions that include mucus production [17][18][19][20] .
The cervicovaginal epithelial barrier acts as the primary site of entry for invading pathogens and barrier disruption is associated with adverse health outcomes such as increased risk of STI (e.g., chlamydia, gonorrhea.HIV) acquisition 21,22 .Undoubtedly, revealing the molecular mechanisms by which vaginal microbes modulate host responses in the lower reproductive tract is necessary to understand how vaginal microbial communities drive reproductive health and disease.Therefore, the objectives of this study were to 1) perform discovery-based RNA-sequencing to unbiasedly identify genes and functional pathways in cervical and vaginal epithelial cells that are altered after exposure to G. vaginalis or L. crispatus and their culture supernatants, 2) elucidate the immune pathways activated by G. vaginalis and 3) identify molecular mechanisms by which L. crispatus optimizes cervical and vaginal epithelial function.

Results
Cervicovaginal epithelial cell gene transcription and associated functional pathways are differentially modulated by L. crispatus, G. vaginalis and their culture supernatants Exposure of cervical and vaginal epithelial cells to live L. crispatus or G. vaginalis or their bacteria-free culture supernatants resulted in signi cant differences in gene expression pro les.In ectocervical, endocervical and vaginal epithelial cells, PCA plots show distinctive clustering by bacterial exposure (Fig. 1A-C).Exposure to L. crispatus culture supernatants and live G. vaginalis were noted to be most distant from their respective controls (NYCIII for culture supernatants and NTC for live bacteria) (Fig. 1A-C).Gene expression pro les appeared to be epithelial cell-speci c as exposure to G. vaginalis culture supernatant resulted in clustering that was different from the NYCIII control in endocervical and vaginal cells that was not seen in ectocervical cells (Fig. 1A-C).To avoid false discovery, we then identi ed differentially expressed genes with an adjusted p-value < 0.05 and a Log2 fold change of greater than 1 or less than − 1 (Supplemental Table 1).Using these criteria, in all three cell types, the number of differentially expressed genes is highest after exposure to L. crispatus culture supernatants followed by live G. vaginalis and G. vaginalis supernatants (Fig. 1D-F).In contrast, very few genes were differentially expressed in all three epithelial cell types after exposure to live L. crispatus (Fig. 1G).A minority of the differentially expressed genes overlap between cell types for each bacterial exposure (Fig. 1G-J).For G. vaginalis exposures (live or supernatant), endocervical and vaginal epithelial cells have the highest number of overlapping differentially expressed genes.By contrast, exposure to L. crispatus culture supernatant elicits the most overlapping differentially expressed genes in ectocervical and endocervical epithelial cells (Supplemental Table 2).Additional comparisons were performed within each cervicovaginal epithelial cell type to identify commonly modulated genes by these different bacterial exposures (Fig. 1K-M, Supplemental Table 3).From these analyses, we were able to identify unique genes where expression was altered by both live G. vaginalis and L. crispatus culture supernatant but were differentially regulated; thus, demonstrating speci c molecular effects of these bacterial exposures on regulating cervicovaginal epithelial cell transcription (Supplemental Table 4, 5, 6).
Gene ontology analysis accounting for both upregulated and downregulated genes was performed on differentially expressed genes in each cell line/exposure combination to uncover the overlapping or different responses between cell types to bacterial exposures (Supplemental Tables 7a-i, notably without live L. crispatus as too few genes changed for analysis) 23,24 .Then, an aggregate score per sample was generated for all the dysregulated pathways and averaged for each cell line to re ect the diversity of gene expression differences between bacterial exposures within individual cell lines by heatmap (Fig. 2A) 23 .Unsupervised clustering of the Gene Ontology terms was used to identify if there were speci c trends.Then, a word cloud of the Gene Ontology terms was generated to describe top terms associated with each cluster 25,26 .Speci c themes of in ammatory and transcriptional dysregulation emerged.To further understand the most important dysregulated pathways in each cell type, we performed a clustering analysis of the Gene Ontology differences by cell line/exposure and focused on the top clusters for further investigation 23 .Live G. vaginalis and G. vaginalis culture supernatant upregulated genes were mostly associated with in ammation functional pathways (Fig. 2B-D).By contrast, exposure to L. crispatus culture supernatant was associated with dysregulation of transcriptional functional pathways, including histone modi cations, RNA polymerase II transcription and DNA binding (Fig. 2A, 2E-G).Live L. crispatus modulated the expression of too few genes to perform functional pathway analysis.Since 1) live G. vaginalis and its culture supernatant modulated many of the same immune functional pathways and 2) L. crispatus culture supernatant resulted in the highest number of differentially expressed genes while live L. crispatus had a minimal effect on epithelial cell gene expression, we chose to focus predominantly on the functional effects of live G. vaginalis and L. crispatus culture supernatant.

G. vaginalis activates an in ammasome-mediated immune response in cervicovaginal epithelial cells
As noted prior, cervicovaginal cells exposed to G. vaginalis resulted in the differential expression of genes signi cantly associated with in ammation-related functional pathways (Fig. 2A, 2B-2D).These genes, including those coding for multiple chemokines and cytokines (e.g., IL-8, IL-6, IL-1α and TNF), have been investigated previously 27 .Additionally, speci c anti-microbial peptides (AMPs), such as CCL20, SLPI, S100A8 and LCN2, were found to be upregulated in cervicovaginal epithelial cells after G. vaginalis but not L. crispatus exposure (Supplemental Fig. 1).Novel to this work, we identi ed several in ammasome-mediated genes that were signi cantly upregulated after G. vaginalis exposure including NLRP3, NLRP1, IL-1β and caspase-1 (adjusted p < 0.05) (Table 1).The gene expression changes for all detectable in ammasome-associated genes are listed in Supplemental Table 8.To determine if in ammasome activation contributes to the in ammatory response seen after G. vaginalis exposure in cervicovaginal epithelial cells, we performed a bacterial dose response experiment (1x10 5 -1x10 7 CFUs/ml) for both L. crispatus and G. vaginalis.Caspase-1, IL-1β, and lactate dehydrogenase (LDH, a marker of cell death) release were signi cantly increased (p < 0.0001) in a dose-dependent manner from ectocervical, endocervical and vaginal cells after exposure to G. vaginalis but not L. crispatus (Supplemental Fig. 2).

L. crispatus culture supernatants alter chromatin accessibility
After exposure to L. crispatus culture supernatants, the most differentially expressed genes in cervical and vaginal epithelial cells were related to alterations in histone and transcriptional regulation (Fig. 2A,   2E-G).To further assess the ability of L. crispatus to induce epigenomics shifts, we used assay for transposase-accessible chromatin high throughput sequencing (ATAC-seq) to ascertain whether these gene expression changes were associated with differences in chromatin accessibility [30][31][32] .
The number and percentage of aligned/unaligned reads were equal across all three cell types (Supplemental Fig. 3).Notably, the samples from the ectocervical cells had lower transcription start site (TSS) enrichment scores despite demonstrating similar quality control characteristics in alignment (Supplemental Fig. 4).There was a strong correlation between normalized read counts between conditions for each cell type (Supplemental Fig. 5), indicating good quality samples.
We obtained a consensus peak set for each cell type (Ecto: 55,917 peaks, Endo: 46,535 peaks, VK2: 48,164 peaks).As expected, the majority of open chromatin peaks for all cell types were found in proximal promoter regions (< 1 kb) or intergenic regions (Fig. 6A).To determine if L. crispatus culture supernatant leads to cell-speci c differences in chromatic organization, we compared the normalized read counts between different genomic regions.There was no difference in quantile normalized counts between the TSS and gene bodies for each cell type (Supplemental Fig. 6).However, only in the ectocervical cells, the consensus peak regions demonstrated nearly two discrete clusters (Fig. 6B).
Endocervical and vaginal cells did not exhibit these two clusters of accessible sites (Fig. 6B).As evidenced by this nding, as well as RNA sequencing demonstrating the most profound shift in epigenetic-based pathways in ectocervical cells, we continued with the ATAC-seq analysis focused on ectocervical cells.
We then tested whether there were regions of differential accessibility between the treatment conditions in the cell type-speci c consensus peak sites (Supplemental Tables 9a-c).We detected 8,147 regions with differential accessibility pro les in L. crispatus supernatant-treated ectocervical cells, almost all showing reduced accessibility (8,125 with decreased accessibility, 22 with increased accessibility).
There was only a fraction of sites differentially regulated in endocervical and vaginal cells (21 and 109 total sites, respectively).Notably, the distribution of differential accessibility sites was skewed towards an increase in the percentage of distal intergenic and intronic regions and a decrease in proximal promoter sites compared to the general distribution of all detected ectocervical peaks (Fig. 6C).We identi ed the genes neighboring the differentially regulated sites and sought to determine whether these genes overlapped with differentially expressed genes we had detected by RNA-sequencing after L. crispatus supernatant.Indeed, most genes with differential expression in ectocervical cells treated with L. crispatus supernatant had differential accessibility by ATAC-sequencing (683/762 genes) [33][34][35] .
Given the strong relationship between differentially expressed genes and differentially accessible chromatin, we sought to identify whether differentially accessible chromatin was more likely to be associated with previously identi ed putative tissue-speci c regulatory regions or all identi ed human enhancer regions (enhancer-like sequences of different primary tissues described by the Encyclopedia of DNA Elements (ENCODE) Project or open accessibility regions identi ed in different primary cancer specimens) [36][37][38] .Surprisingly, we discovered little overlap between differential accessibility in ectocervical cells and putative enhancer regions (Fig. 6D, Chi-squared test for trend p = 0.8334 in differentially accessible ectocervical peaks vs random matched control regions).However, ENCODE does not have specimens from primary cervical tissue; as such our comparative analyses was limited to specimens collected from the vagina and uterus.By contrast, there was a substantial overlap between all open chromatin regions identi ed in ectocervical cells and published likely enhancers (Fig. 6E, Chisquared test for trend p = 0.0028 in all accessible ectocervical peaks vs random matched control regions).
Despite not being associated with known enhancer regions, motif analysis of the downregulated differentially accessible sites demonstrated enrichment for 497 transcription factor motifs with an FDR of 0.05 (Supplemental Table 10, Fig. 6F).The position of the top 5 identi ed motifs was at the center of the peak, consistent with the expected location of any true transcription factor binding 39 .Gene-disease enrichment analysis of these transcription factors demonstrated a marked enrichment of multiple pathways related to neoplasms, endometriosis and infertility, all pathologies potentially alleviated or protected by Lactobacillus-dominated microbiota (Fig. 6G) 40,41 .To determine if this enrichment was speci c, we generated random lists of transcription factors from the same motif database and demonstrated only an overlapping enrichment in unrelated craniofacial anomalies, suggesting our ndings are indeed speci c (Supplemental Fig. 7).

Discussion
Host-microbe interactions have been shown to govern health and disease in multiple biological systems.
With this study, we now reveal unique molecular mechanisms involved in host-microbe interactions in the cervicovaginal space, thus addressing a large knowledge gap in reproductive health.These ndings con rm the induction of diverse immune pathways in response to G. vaginalis, a facultative anaerobic bacteria associated with many gynecological disorders, including STIs 12,13 , cervical cancer 42,43 , infertility 10,11 , and preterm birth 7,15 Importantly, we found that G. vaginalis activates the in ammasome via NLRP3 in cervicovaginal epithelial and immune cells.Suggesting a molecular bene t from colonization with Lactobacillus species, we report that exposure to L. crispatus culture supernatants results in epigenetic modi cations in ectocervical cells.Collectively, these studies demonstrate the complexity of host-microbe interactions in the lower reproductive tract and demonstrate speci c molecular mechanisms by which optimal and non-optimal bacteria contribute to reproductive health and disease.
While many studies have used high throughput sequencing technologies to identify and characterize microbial communities present in the cervicovaginal space [44][45][46] , few have leveraged whole genome transcriptomic sequencing to investigate the genes and functional pathways altered by host-microbe interactions between different vaginal bacteria and their supernatants on all three cervicovaginal epithelial cells types.The RNA-seq results from this study reveal host epithelial cell-speci c genes and functional pathways that are modulated by G. vaginalis and L. crispatus.Whole-genome transcriptomic evaluations of ectocervical, endocervical and vaginal cells exposed to G. vaginalis or L. crispatus (or their culture supernatants) revealed unexpected ndings in that L. crispatus supernatant exposure modulated a high number of genes while live L. crispatus modulated very few genes.As studies have shown that L. crispatus-dominated microbiota protect the cervicovaginal epithelium from viral and bacterial infection 45,47,48 , it has been hypothesized that L. crispatus secreted factors contribute to this protection.For example, one study identi ed that culture supernatant from Lactobacillus species, including L. crispatus, is protective against Chlamydia trachomatis infection due to the protective effects of D(-) lactic acid and a reduction in vaginal epithelial cell proliferation 45 .It is feasible that the epigenetic modi cation observed in our study was due to lactic acid in the L. crispatus microbial supernatants.
Additionally, we have recently demonstrated that L. crispatus supernatants contain functional extracellular vesicles that can interact with host epithelial cells providing another potential mechanism for host-microbial interactions in the cervicovaginal space 49 .Further work is needed to specify the component(s) of L. crispatus supernatants that drive the epigenetic shifts observed in this study.In contrast to live L. crispatus, RNA sequencing revealed that live G. vaginalis altered the expression of a high number of genes, with live G. vaginalis modulating more than its culture supernatant.The difference in gene expression pro les induced by live bacteria or their culture supernatants by these two common bacterial species demonstrates the complexity of cervicovaginal host-microbial interactions.
Understanding how individual microbes and their by-products, as well as how they function in larger microbial communities, is needed to fully ascribe mechanisms by which the vaginal microbiome promotes reproductive health or disease.
The lower reproductive tract is unique in the diversity of epithelium that line this biological niche.
Consistent with the different expected functions of these unique epithelial barriers, speci c transcriptomic signatures were observed for each type of cervicovaginal epithelial cell analyzed.After G. vaginalis exposure, distinct sets of genes were upregulated in each epithelial cell type, indicating that most altered genes appear to be cell type-speci c; however, many of these genes target similar innate immune pathways, providing evidence of the redundancy in immune regulation.The innate immune pathways most predominantly implicated in the cervical and vaginal epithelial cell response to G. vaginalis included the upregulation of NF-kB signaling, an increase in anti-microbial peptides (AMPs) and activation of the NLRP3 in ammasome.The nding that G. vaginalis initiates an NF-kB-mediated in ammatory response in cervicovaginal epithelial cells is not unexpected as previous studies by our laboratory, and many others, have shown signi cant increases in cervical and vaginal epithelial cell cytokine levels after G. vaginalis exposure 27,50,51 .Interestingly, a recent study investigated the vaginal epithelial proteome after exposure to G. vaginalis culture supernatants and identi ed an activation of the mTOR signaling pathway 52 .As mTOR is known to regulate many immune functions including the activation of NF-kB 53 , it is possible that mTOR may play a role in G. vaginalis host-microbial functions, however, we did not detect a change in mTOR gene expression in vaginal cells after exposure to either live bacteria or culture supernatants in our RNA-seq analysis.Even though many studies have identi ed an activated immune response in Lactobacillus-deplete microbial communities 27,54,55 , characterization of this response has been highly varied and complex putting even more emphasis on the need to clearly identify speci c immune factors/pathways activated by cell and microbe speci c interactions.
Interestingly, RNA-seq identi ed a signi cant upregulation of NOD-like receptor (NLR) family, pyrin domain-containing protein 3 (NLRP3) in ammasome-associated genes in cervicovaginal cells after G. vaginalis exposure.In ammasomes are a group of proteins that act as intracellular pathogen recognition receptors (PRRs) that respond to pathogen associated molecular patterns (PAMPs) from foreign microbes.As bacteria have evolved to avoid canonical host immune detection by cell surface TLRs and NF-kB-mediated signaling mechanisms, they have developed pathogen-speci c virulence factors such as pore-forming toxins (ex.vaginolysin for G. vaginalis) 56 , among others, which can be used to gain access into the host cell.In ammasomes, acting as an essential part of the innate immune defense, are located in the cytosol of host cells and are poised to act quickly upon the detection of bacterial pathogens or their by-products (proteins, nucleic acids, metabolites) 57 .In ammasome activation ultimately leads to the secretion of in ammatory factors and chemokines, leukocyte in ltrations, lymphocyte activation and often pyroptosis (cell death) of the host cell 58,59 .While in ammasomes are bene cial to pathogen clearance, constitutive in ammasome activation has been shown in a multitude of chronic in ammatory conditions including Alzheimer's disease, asthma, cancers, heart diseases, rheumatoid arthritis, and diabetes 60 .There is a paucity of data on microbe-mediated in ammasome activation in the lower reproductive tract.However, the role of in ammasomes in both pathogen clearance and activation of the in ammatory response in other body sites suggests that the in ammasome could play a signi cant role in host-microbe immune interactions at the cervicovaginal epithelial barrier.Recent work has demonstrated that a known pathogen (Neisseria gonorrhoeae) and a vaginal anaerobe (G.vaginalis) activate the NLRP3 in ammasome in immune cells [61][62][63] .Adding to this knowledge, we con rm the ability of G. vaginalis to activate the in ammasome in immune cells and note the speci city of G. vaginalis (but not L. crispatus) to induce this in ammatory pathway.
Novel to this work, we now show that that G. vaginalis activates a canonical in ammasome pathway, not just in immune cells, but also in cervicovaginal epithelial cells.G. vaginalis activation of the in ammasome pathway seems to be mediated, at least in part, through NLRP3, as the NLRP3 speci c inhibitor, glybenclamide, was able to reduce the G. vaginalis-mediated increases in caspase-1 and IL-1β.While prior reports demonstrate that NLRP3 activation is essential for the induction of the in ammasome by G. vaginalis for immune cells, we nd that the induction of the in ammasome in cervicovaginal epithelial cells is more complex 62,64 .In contrast to what was observed in immune cells, the NLRP3 inhibitor only partially abrogated G. vaginalis activation of the in ammasome in epithelial cells.As RNA-seq showed signi cant increases in expression of several genes associated with the in ammasome pathway (e.g.NLRP1), it is possible that in epithelial cells, the in ammasome can be activated through non-NLRP3 mechanisms (Supplemental Table 8).As evidence of this, one study showed that Yersinia pseudotuberculosis activated the NLRP3 in ammasome in macrophages, however, in intestinal epithelial cells, it activated the non-canonical caspase-4 in ammasome suggesting that the same bacteria can induce different in ammasome-mediated pathways in different cell types 65 .The nding that G. vaginalis can activate the NLRP3 in ammasome in both cervicovaginal epithelial cells and monocytes suggests that in ammasome activation may be an important and common contributor to the in ammatory response in a G. vaginalis-dominated cervicovaginal space.
While cell death (pyroptosis) is an integral part of in ammasome activation, G. vaginalis also secrets a cholesterol-dependent cytolysin, vaginolysin, that creates pores in the epithelial cell membrane causing cell death 56 .Therefore, it is di cult to discern if the cell death observed after G. vaginalis exposure is due to in ammasome activation directly, an indirect result of vaginolysin or a combination of both.While our results show that cell death in cervicovaginal cells is mediated by caspase-1 activation, in agreement with studies that show caspase-1 is required for Gasdermin cleavage and activation (necessary for NLRP3-induced pyroptosis) 66 , we did not see a reduction in cell death with the NLRP3 inhibitor in either cervicovaginal cells or monocytes.Further elucidation of the mechanisms leading to cell death would be important for therapeutic strategies to limit G. vaginalis-mediated epithelial dysfunction.
Demonstrating the capabilities of different vaginal bacteria species to induce unique molecular pro les, cervicovaginal cells exposed to L. crispatus culture supernatant predominantly showed alterations in RNA polymerase, histone, and transcription-related genes.These ndings suggest an important role for L. crispatus culture supernatant in inducing epigenomic modi cations of host epithelial cells.While an abundance of Lactobacillus species (including L. crispatus, L. jensenii and L. gasseri) in the cervicovaginal space have been shown to be protective against microbial pathogens, vaginal infection, and cervical cancer, 42 the biological mechanisms contributing to Lactobacillus protection remain largely unknown thus limiting the ability to leverage this molecular effect for therapeutic bene t.Unique to this study, we nd that L. crispatus supernatant caused substantial modulation of genes globally regulating the transcriptome and epigenome.These changes corresponded to a substantial unexpected reorganization of the epigenome as shown by ATAC-seq.L. crispatus largely reduced the number of open chromatin regions in ectocervical cells.An intriguing hypothesis is that L. crispatus may be important for increasing the resilience of cells to infection (ex.Chlamydia, HIV, HPV) by modulating the epigenetic susceptibility of the cells to a disease or pathogen 67,68 .Providing evidence for the ability of L. crispatus to alter the epigenome to protect from disease, a previous study has shown that supernatant from Lactobacillus species can decrease C. trachomatis infection by inhibiting cell proliferation through epigenetic mechanisms.Speci cally, L. crispatus culture supernatants decreased histone deacetylase 4 (HDAC4) and increased histone acetylase EP300 through mechanisms thought to involve D(-) lactic acid in vaginal cells 45 .Based on this report, we hypothesize that L. crispatus culture supernatants induce these epigenomic effects due to the presence of lactic acid isomers.However, it is also plausible that other metabolites and/or proteins within the L. crispatus culture supernatant are responsible for inducing these epigenomic changes.Additionally, we have recently shown that bacterial extracellular vesicles (bEVs), released from the bacteria into the culture supernatant, can carry proteomic (and other biological cargo) and can interact with host cervical epithelial cells 49 .These bEVs could also carry modulating factors resulting in epigenomic alterations to the host cell.Importantly, understanding the ability of L. crispatus supernatants to modify cervicovaginal epithelial function may provide new therapeutic strategies to optimize reproductive health.
Intriguingly, in this study the regions of differential chromatin accessibility did not overlap with regions previously described in published putative enhancers.This discrepancy could be due to multiple reasons.First, we do not have ENCODE or ATAC pro les of healthy primary cervical tissue for comparison.In addition, the pro led cervical cancer specimens were derived from only four specimens and likely have different biological pro les than the ectocervical lines used here.Therefore, it is possible that these differentially regulated regions are cell type-speci c enhancers in ectocervical cells.Second, it is notable that there is a substantial increase in the proportion of intronic sites among the ectocervical differentially accessible sites compared to all the ectocervical peaks (48.2% vs 26.7% respectively), suggesting a potential role for Lactobacillus in regulating isoform transcriptions through chromatin modulation.However, RNA-seq data in this study is limited in its ability to directly answer this question; a repeat of this study with long-read RNA-seq would be more suitable 69 .Lastly, there is a growing literature about the roles of lactate as a precursor to acetyl-CoA needed for histone acetylation as well as a direct modi er of histones through histone lactylation [70][71][72][73] .Histone lactylation is generally described to lead to increased expression, putatively associated with more open chromatin.Yet, the role of lactate in modifying histones and other epigenetic regulators is a new area of research and further work would be needed to elucidate its role in decreasing chromatin accessibility in speci c contexts.Nonetheless, while we cannot fully identify how these differentially regulated sites may be contributing to gene regulation, it is notable that disease gene enrichment analysis of transcription factors associated with the motifs at these sites demonstrated multiple pathologies related to women's health, including fertility and endometriosis, posing an intriguing avenue for better understanding the molecular underpinnings of these common but poorly understood disorders.
A limitation of this study is that it relies on single strains of common vaginal bacteria.While studies investigating broader microbial communities are needed to more accurately mimic the cervicovaginal microbiota present in vivo, by focusing on single microbes of interest we were able to identify speci c functional pathways that may drive adverse outcomes.An in-depth understanding of single bacteria lays the foundation for these future studies.In addition to G. vaginalis, studies investigating other high-risk anaerobic bacteria for adverse reproductive outcomes 7 , including Sneathia species, Mobiluncus species and Prevotella species would be necessary to elucidate the combinatorial effects of anaerobes more broadly on cervicovaginal function.Now that speci c biological functions of both G. vaginalis and L. crispatus have been identi ed, these pathways can be used to develop mechanistic hypotheses to evaluate more complex microbiota or clinical populations.For the analysis of ATAC data, publicly available databases, such as the ENCODE database, do not contain comparable ATAC data in cervical tissue or vaginal tissue for comparison of normal/control human tissue 36 .Developing these resources should be a priority as they will be critical for further understanding of the host-microbial interactions in the reproductive tract and their role in diverse reproductive outcomes.
Overall, the results of this study identify novel transcriptomic and epigenomic pathways altered by microbes within the cervicovaginal space that are most commonly associated with reproductive health and disease.Utilizing whole genome RNA-sequencing, we identi ed microbe-speci c functional gene pathways including activation of the innate immune response by G. vaginalis and increased RNA transcription and histone modi cations by L. crispatus that are regulated by host-microbe interactions within the cervicovaginal space.These results suggest that complex interactions between host cells and live bacteria or the factors they secrete have distinct and speci c functions in modifying host epithelial cell responses.Novel to this study, we identi ed two key areas for potential therapeutic targets: 1) activation of the in ammasome as part of the innate immune response to G. vaginalis and 2) epigenetic regulation by L. crispatus culture supernatants.Targeting speci c G. vaginalis-mediated innate immune pathways may serve to modulate the in ammatory response associated with G. vaginalis, and thus, could impact STIs, bacterial vaginosis and preterm birth.Likewise, leveraging the potential of L. crispatus to alter the epigenetic landscape may provide new opportunities to optimize the cervicovaginal epithelial barrier and prevent pathogenic (e.g.chlamydia, HPV) microbes from harming and infecting the epithelial barrier.Continued elucidation of host-microbial interactions in the female reproductive tract will undoubtedly serve to optimize reproductive health.

Bacterial Cultures and Preparation of Bacteria-Free Supernatants
Human clinical isolates of L. crispatus (ATCC 33197) or G. vaginalis (ATCC 14018), were obtained from the American Type Culture Collection (Manassas, VA).G. vaginalis was grown on Tryptic Soy Agar with 5% Sheep Blood plates (Hardy Diagnostics) and L. crispatus was grown on De Man, Rogosa and Sharpe agar (Fisher Scienti c); both strains were grown in New York City III (NYCIII) broth.Bacteria were grown at 37 o C in an anaerobic glove box (Coy Labs, Grass Lake, MI).
For each experiment the following bacterial growth protocol was followed: L. crispatus and G. vaginalis glycerol stocks were streaked on agar plates, as well as into broth tubes and grown overnight.The broth starter cultures were diluted to an optical density of 0.2 and then used to inoculate 20ml working cultures, which were grown for 20 hours (G.vaginalis) to 48 hours (L.crispatus) prior to use in experiments.Bacterial densities of the working cultures were estimated the day of the experiment based on optical density readings at 600 nm using an Epoch2 plate reader (Biotek, Winooski, VT), and the appropriate volume was centrifuged at 13,000 x g for 3 min.The bacterial pellets were resuspended in the appropriate cell culture media without antibiotics and added to epithelial cells at 10 5 -10 7 CFUs/well.Precise bacterial densities of the working cultures were determined by plating serial dilutions of the working cultures and counting CFUs.For all experiments, reported bacterial densities are +/-0.5 log of the noted bacterial density (CFU/well).
To obtain bacteria-free culture supernatants, the working cultures were centrifuged at 13,000 x g for 3 min and the supernatant was ltered through a 0.22 µm lter (Fisher Scienti c) to remove any remaining live bacteria.Bacteria-free culture supernatants were diluted to 1% v/v in the appropriate cell culture media without antibiotics.
In vitro Epithelial and Immune Cell -Bacteria Interactions Ectocervical, endocervical and vaginal cells were plated at 1.5 x 10 5 cells/well in twenty-four well plates containing KSFM without antibiotics.THP1-Null2 or THP1-KO-NLRP3 cells were plated at 1.5 x 10 5 cells/well in twenty-four well plates containing RPMI 1640 media without antibiotics.The next day, the cells were exposed to either live L. crispatus or G. vaginalis or 1% (v/v) bacteria-free supernatants (generated from a 1x10 7 CFU/mL culture) for 24 hr.For live bacteria experiments, a 1x10 5 CFU/well dose of bacteria was used for epithelial cell and THP1-Null or KO experiments.This dose was chosen based on LDH dose responses in epithelial 27 and THP1-Null/KO (Supplemental Fig. 8) cells.Both L. crispatus and G. vaginalis bacteria survive in a 5% CO 2 culture with cervicovaginal epithelial cells in KSFM media for up to 24 hours (Supplemental Fig. 9) as shown by CFUs done at the end of the incubation period.For all experiments where cells are exposed to live bacteria, non-treated control (NTC) media (KSFM or RPMI) without bacteria is used as the control.Bacteria-free supernatant percentage was based on a dose response (1% vs 10%) (Supplemental Fig. 10).For cells exposed to 1% bacteria-free supernatants from L. crispatus, KSFM media was supplemented with 50mM HEPES and sodium bicarbonate (3000 mg/L total concentration) to bring the pH of the media up to a physiological level (7.2).In additional experiments, ectocervical (n = 6/treatment), endocervical (n = 6/treatment), vaginal cells (n = 6/treatment) or THP1-Null2 (n = 6/treatment) cells were pre-treated with Glybenclamide (100 µg/ml, Invivogen), a speci c NLRP3 inhibitor, for 30 mins prior to live bacteria exposure (1x10 7 CFU/well) or Ac-YVAD-cmk (25 µM, Sigma-Aldrich), an irreversible caspase-1 inhibitor, for 1 hour prior to live bacteria exposure (1x10 7 CFU/well).For in ammasome experiments, LPS (10 ng/ml) exposure for 24 hours was used as a positive control of in ammasome activation.For all supernatant experiments, cells were also exposed to 1% (v/v) NYCIII bacterial growth media alone (diluted in KSFM) to determine any baseline effects of the bacterial growth media on the outcomes of interest.1% NYCIII (NYC) acted as the control for all bacteria-free supernatant exposures.At the end of each experiment, cell culture media was collected for cell death, ELISA assays and/or the cells were collected in Trizol (Invitrogen, Thermo-Fisher Scienti c) for RNA extraction.Supplemental Fig. 11 details the methodology of live bacteria and bacteria-free supernatant exposures of cervicovaginal epithelial and THP-1 cells.

RNA Sequencing and Analysis
RNA was extracted from ectocervical, endocervical and vaginal cells after exposure to bacteria or culture supernatants from L. crispatus and G. vaginalis (n = 3/treatment group) collected in Trizol using the Qiagen-RNeasy Plus Mini kit by the Penn Next-Generation Sequencing Core.The resulting RNA had RIN values > 9. Illumina sequencing libraries were prepared using the Illumina TruSeq mRNA stranded library prep kit according to the manufacturer recommendations.The resulting libraries had an average molarity of 69 nM +/1 27 nM.Libraries were sequenced to a median depth of 41 million 100 bp single reads on an Illumina NovaSeq 6000.Transcript quanti cation from RNA-seq data was performed using Salmon and release 38 (GRCh38.p13) of the human genome 74,75 .Several Bioconductor packages in R were used for subsequent steps 76,77 .The output was annotated and summarized using tximeta and further annotation was completed with biomaRt 78,79 .Principle Component Plots (PCA) were created using pcaExplorer 80,81 .The normalizations and statistical analyses were done with DESeq2 82 .The full RNA-seq dataset was submitted to Gene Expression Omnibus (accession # GSE234837).

Cell Death Assay
Ectocervical, endocervical, vaginal, THP1-Null2 and THP1-KO-NLRP3 cells were grown and exposed to bacteria or bacterial culture supernatants as described above.Lactate dehydrogenase (LDH) released upon cell lysis (n = 3-9 independent experiments per cell type) was quanti ed using the CytoTox 96 Nonradioactive cytotoxicity assay (Promega, Madison, WI), a coupled enzymatic assay that results in the conversion of a tetrazolium salt into a red formazan product.The colorimetric output was measured using a plate reader at 490 nm and absorbance values were recorded.

ELISA
Ectocervical, endocervical, vaginal, THP1-Null2 and THP1-KO-NLRP3 cells were cultured in 24-well plates and exposed to either live bacteria or bacterial supernatants as stated above.Anti-microbial peptides, CCL20, SLPI, LCN2, S100A8/A9, or in ammasome-associated cytokines, IL-1β and caspase-1, were measured in cell culture media after 24 hours of exposure (n = 6/group).The expression of these analytes was measured by a ligand-speci c commercially available ELISA kit that utilizes a quantitative sandwich enzyme immunoassay technique using reagents from R&D Systems (Minneapolis, MN).
ATAC-seq Nuclei Extraction, Tagmentation, Puri cation and Library Ampli cation ATAC-seq was performed on ectocervical, endocervical and vaginal cells after exposure to L. crispatus bacteria-free supernatants (n = 3/treatment group).ATAC-seq libraries were generated using the ATACseq Kit from Diagenode (Diagenode, A Hologic Company) according to manufacturer instructions.Brie y, nuclei were extracted from 50,000 cells.Tagmentation was completed by resuspending the isolated nuclei in transposase reaction mix and the samples were puri ed using the kit's provided columns.Following puri cation, library fragments were ampli ed by PCR according to the manufacturer recommendations.Unique Dual Indexes Primer Pairs were incorporated for multiplexed sequencing.To reduce ampli cation bias, after the rst 5 cycles of the PCR reaction, qPCR was used to determine how many additional cycles were needed to produce enough library to meet the required amount for sequencing.For this, an aliquot of the PCR reaction was added to Sybr Green and ampli ed for 20 cycles.Libraries were ampli ed for a total of 11-13 cycles (with one library requiring 17 cycles for ampli cation).Final libraries were puri ed using bead puri cation (Beckman Coulter), then assessed for size distribution and concentration using a BioAnalyzer High Sensitivity DNA Kit (Agilent Technologies).
The resulting libraries were pooled.The pool was diluted to 2 nM, denatured, and the 13 libraries were loaded onto an S1-100 (2x50) ow cell on an Illumina NovaSeq 6000 (Illumina, Inc.) according to the manufacturer's instructions.The average read number per sample was 50M+/-20%.De-multiplexed and adapter-trimmed sequencing reads were generated using bcl2fastq.The full ATAC-seq dataset was submitted to Gene Expression Omnibus (accession # GSE233444).

ATAC-seq Mapping and Peak Calling
ATAC-seq data analysis was adapted from a previously published approach using PEPATAC (v.0.10.3) 30.Peaks for each cervicovaginal line were called separately to allow for cell type speci c differences in chromatin accessibility pattern.In brief, raw FASTQ les were processed and mapped to release 38 (GRCh38.p13) of the human genome using the PEPATAC pipeline 83 .Reads were trimmed with skewer and then aligned with bowtie2 using default settings 84,85 .Duplicate reads were removed using samblaster 86 .
An iterative overlap peak calling strategy on xed-sized peaks of 501 bp was used to de ne a set number of peaks for each cell type for downstream differential accessibility comparison 30 .First, for each biological replicate, MACS2 was used to call peaks with the parameters as follows: --peak-type xed -extend 250 87 .Biological replicates of each treatment and then both treatments together from each cell type were merged using an iterative overlap approach previously described 30,38 .Blacklisted regions were excluded from called peaks (accessed 4 November 2022 at https://github.com/Boyle-Lab/Blacklist) 88.

ATAC-seq Peak and Differential Accessibility Analysis
Peak location was annotated with CHIPseeker (v 1.30.3) 89.Counts for peaks were calculated using Rsubread (v.2.8.2) 90 .We determined the differential accessibility of peaks between treatments with DESeq2 (v.1.34.0) 91.We compared L. crispatus culture supernatant treated to NYCIII media controls for each cell type.A Wald test was used to determine signi cance.A peak was de ned as statistically signi cant in differential accessibility if |log2foldchange| > 1 and FDR < 0.05.We utilized the R package rGREAT (v.1.99.0) for the nearest gene analysis to access the Genome Regions Enrichment of Annotations Tool (GREAT) web service [33][34][35] .For GREAT, we used the parameters for "the two closest genes" to a differential accessible site as it is frequently not the closest genes that is differentially regulated.
Motif analysis was performed using Simple Enrichment Analysis version 5.5 as part of the MEME Suite (https://meme-suite.org/meme/tools/sea) 92,93 .Differentially accessible sites were inputted, and the CIS-BP 2.0 motifs database were used for the query 94 .Gene-disease enrichment analysis was performed using disgenet2R (v.0.99.2, https://www.disgenet.org/) 95.Random gene lists were generated for comparison by sampling 497 transcription factors from the CIS-BP 2.0 database to ascertain the baseline disease enrichment bias of the database.

ATAC-seq Chromatin Accessibility Visualization
EasSeq (v1) was utilized to visualize the data 96 .Biological replicates of BAM les were pooled for quanti cation of speci c regions.Quantile normalization was used for counts per region for visualization to minimize bias from sequencing depth.Calculation of overlap performed both by any amount of overlap and the exact overlap of base pairs between all comparisons.Random regions for comparisons to differentially accessible regions or all ectocervical open chromatin regions were generated by Regulatory Sequence Analysis Tools matched for each cell type by number of fragments, fragment size, and GC content (random genome fragments tool; http://rsat.sb-roscoff.fr/) 97.
ENCODE datasets for all human enhancer-like sequences (ELS, de ned as high DNAse-seq signal and high H3K37me3), or tissue-speci c regulators were obtained from https://screen.encodeproject.org/ 36,37.For uterus and vaginal specimens, "Low-DNase" were ltered out to enrich for sites that had any evidence of potential enhancer or regulator activity.However, strict enhancer-like signature criteria could not be applied because all sequencing modalities were not available for all the samples.Primary cancer cell data sets were obtained from supplemental of published ATAC pro ling 38 .Chi square analysis of number of overlapping sites was performed by Graph Pad.

Figures
Figures

Figure 3 Live
Figure 3

Figure 4 In
Figure 4

Figure 5 Live
Figure 5