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 significant differences in gene expression profiles. 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 profiles appeared to be epithelial cell-specific 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 identified 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 specific 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 reflect 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 specific trends. Then, a word cloud of the Gene Ontology terms was generated to describe top terms associated with each cluster 25,26. Specific themes of inflammatory 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 inflammation functional pathways (Fig. 2B-D). By contrast, exposure to L. crispatus culture supernatant was associated with dysregulation of transcriptional functional pathways, including histone modifications, 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 inflammasome-mediated immune response in cervicovaginal epithelial cells
As noted prior, cervicovaginal cells exposed to G. vaginalis resulted in the differential expression of genes significantly associated with inflammation-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, specific 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 identified several inflammasome-mediated genes that were significantly 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 inflammasome-associated genes are listed in Supplemental Table 8. To determine if inflammasome activation contributes to the inflammatory response seen after G. vaginalis exposure in cervicovaginal epithelial cells, we performed a bacterial dose response experiment (1x105-1x107 CFUs/ml) for both L. crispatus and G. vaginalis. Caspase-1, IL-1β, and lactate dehydrogenase (LDH, a marker of cell death) release were significantly 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).
G. vaginalis activates the NLRP3 inflammasome in THP monocytes
To determine the specificity of NLRP3 in G. vaginalis-induced inflammasome activation in cervicovaginal epithelial cells, we treated monocytes, THP1-Null2 cells or THP1-NLRP3-KO (knock out) cells, with LPS (positive control), live L. crispatus (control for live G. vaginalis) or live G. vaginalis. Exposure to live G. vaginalis, but not L. crispatus, resulted in a significant increase in caspase-1 (p < 0.001, Fig. 3A) and IL-1β (p < 0.0001, Fig. 3B) in THP1-Null2 cells but not in THP-NLRP3-KO cells (p < 0.001). Cell death (LDH, Fig. 3C) remained unchanged after bacterial exposure in both THP1 cell types. Additionally, we treated THP1-Null2 cells with glybenclamide (glyburide) 28, a specific inhibitor of NLRP3 (with no effects on NLRC4 or NLRP1). Glybenclamide treatment significantly reduced caspase-1 (p < 0.05) and IL-1β (p < 0.0001, Fig. 3D and 3E). As with THP1-NLRP3-KO, no reduction in G. vaginalis-mediated cell death (LDH, Fig. 3F) was observed in THP1-Null2 cells treated with the NLRP3 inhibitor.
G. vaginalis -mediated inflammasome activation in cervicovaginal epithelial cells is NLRP3 specific
Since G. vaginalis activation of the inflammasome appeared to be NLRP3-dependent in THP1 monocytes, we assessed if NLRP3 was essential to G. vaginalis activation of the inflammasome in cervicovaginal epithelial cells. Glybenclamide treatment resulted in a significant reduction in G. vaginalis-mediated increases in caspase-1 (p < 0.01, Fig. 4A, 4D, 4G) and IL-1β (p < 0.05, Fig. 4B, 4E, 4H) in ectocervical, endocervical and vaginal cells, while G. vaginalis-induced cell death, as measured by LDH release, was unchanged (Fig. 4C, 4F, 4I).
G. vaginalis activates the inflammasome-mediated immune response through caspase-1
Since inhibiting NLRP3 activation did not fully mitigate G. vaginalis activation of the inflammasome in cervicovaginal epithelial cells, we sought to determine if inhibiting caspase-1 could prevent inflammasome activation. The irreversible caspase-1 inhibitor, Ac-YVAD-cmk 29, resulted in a significant reduction in the G. vaginalis-mediated increase in caspase-1 (p < 0.0001) (Fig. 5A, 5D, 5G), IL-1β (p < 0.0001) (Fig. 5B, 5E, 5H) and cell death (p < 0.0001) (LDH, Fig. 5C, 5F, 5I) in ectocervical, endocervical and vaginal cells after G. vaginalis but not L. crispatus exposure.
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–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-specific 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 finding, 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-specific consensus peak sites (Supplemental Tables 9a-c). We detected 8,147 regions with differential accessibility profiles 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 identified 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–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 identified putative tissue-specific regulatory regions or all identified human enhancer regions (enhancer-like sequences of different primary tissues described by the Encyclopedia of DNA Elements (ENCODE) Project or open accessibility regions identified in different primary cancer specimens) 36–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 identified in ectocervical cells and published likely enhancers (Fig. 6E, Chi-squared 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 identified 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 specific, we generated random lists of transcription factors from the same motif database and demonstrated only an overlapping enrichment in unrelated craniofacial anomalies, suggesting our findings are indeed specific (Supplemental Fig. 7).