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 and indicates the greatest differences (Fig. 1A-C) between exposures. 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). Figure 1G-J depicts the number of overlapping differentially expressed genes between cervicovaginal cell types for each bacterial exposure. For G. vaginalis exposures (live or supernatant), endocervical and vaginal epithelial cells have the highest number of overlapping differentially expressed genes, while exposure to L. crispatus culture supernatants illicit 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). In summary, genes that were upregulated by exposure to G. vaginalis and downregulated by L. crispatus point towards a pro-inflammatory response. On the other hand, those upregulated by exposure to L. crispatus and downregulated by G. vaginalis seem to be involved in transcription, epigenetic modifications, cell-cell adherence, and anti-inflammatory functions.
Functional pathway analysis using Gene Ontology terms identified specific pathways of interest between bacterial exposures. Live G. vaginalis and G. vaginalis culture supernatant upregulated genes were mostly associated with inflammation functional pathways (Fig. 2A, 2B-D). Exposure to L. crispatus culture supernatant was associated with upregulation 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 the same 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 on exposures of live G. vaginalis and L. crispatus culture supernatant for further analysis.
Live G. vaginalis and its culture supernatant differentially up-regulate anti-microbial peptide gene expression and increase protein levels
RNA-seq functional pathway analysis identified "defense response to a bacterium" along with various inflammation-based pathways as Gene Ontology terms of significance after live G. vaginalis exposure (Fig. 2B-D). The specific genes associated with these pathways included the anti-microbial peptides, Chemokine Ligand 20 (CCL20), Secretory Leukocyte Peptidase Inhibitor (SLPI), Lipocalin 2 (LCN2) and S100 Calcium Binding Protein 8 (S100A8/A9, Calgranulin). Our data show that all four of these genes were significantly upregulated (adjusted p < 0.05) by live G. vaginalis exposure in ectocervical and endocervical cells with only CCL20 and SLPI being upregulated in vaginal cells (Fig. 3A, 3E, 3I). After exposure to G. vaginalis culture supernatants, CCL20 was the only transcript upregulated in all three cell lines (adjusted p < 0.05, Fig. 3B, 3F, 3J). In ectocervical cells, in addition to CCL20, S100A8 gene expression was upregulated after G. vaginalis culture supernatant exposure (adjusted p < 0.05, Fig. 3B). In endocervical cells, all four anti-microbial peptide genes were significantly upregulated by G. vaginalis culture supernatant (adjusted p < 0.05, Fig. 3F). Only expression of CCL20 was upregulated in vaginal cells after G. vaginalis culture supernatant exposure (adjusted p < 0.05, Fig. 3J). ELISAs for these antimicrobial peptides confirmed overexpression of CLL20 and S100A8 (both p < 0.05), but not SLPI and LCN2, which remained unchanged, after exposure to G. vaginalis (Fig. 3C, 3G, 3K) and G. vaginalis culture supernatants (Fig. 3D, 3H, 3L).
Down-regulation of anti-microbial peptide gene expression, but not protein, by live L. crispatus and its culture supernatant
Live L. crispatus did not affect the gene expression of the four anti-microbial peptides measured (Fig. 3A, 3E, 3I), except that CCL20 mRNA was downregulated (p < 0.0001) in endocervical cells; however, the amount of CCL20 protein was not altered at this same time point (Fig. 3G). L. crispatus culture supernatants downregulated the gene expression of LCN2 (p < 0.05) and S100A8 (p < 0.001) in endocervical and vaginal cells (Fig. 3F and 3J), but not ectocervical cells (Fig. 3B). Additionally, CCL20 (p < 0.0001) gene expression was downregulated in endocervical cells (Fig. 3F). However, this decrease in anti-microbial gene mRNA expression in endocervical and vaginal cells was not associated with a change in their associated proteins (Fig. 3H and 3L).
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, included those coding for multiple chemokines and cytokines (e.g., IL-8, IL-6, IL-1α and TNF), have been investigated previously 23. 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 7. 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. 1).
G. vaginalis activates the NLRP3 inflammasome in THP monocytes
To determine the specificity of NLRP3 inflammasome activation in mediating G. vaginalis-induced increases in IL-1β, caspase-1 or cell death 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. 6A) and IL-1β (p < 0.0001, Fig. 6B) in THP1-Null2 cells but not in THP-NLRP3-KO cells (p < 0.001). Cell death (LDH, Fig. 6C) remained unchanged after bacterial exposure in both THP1 cell types. Additionally, we treated THP1-Null2 cells with glybenclamide (glyburide) 24, 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. 6D and 6E). As with THP1-NLRP3-KO, no reduction in G. vaginalis-mediated cell death (LDH, Fig. 6E) 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 activates the NLRP3 inflammasome in THP1 monocytes, we determined if G. vaginalis could also activate an NLRP3-specific inflammasome in cervicovaginal epithelial cells. Glybenclamide treatment resulted in a significant reduction in G. vaginalis-mediated increases in caspase-1 (p < 0.01, Fig. 5A, 5D, 5G) and IL-1β (p < 0.05, Fig. 5B, 5E, 5H) in ectocervical, endocervical and vaginal cells, while G. vaginalis-induced cell death, as measured by LDH release, was unchanged (Fig. 5C, 5F, 5I).
G. vaginalis activates the inflammasome-mediated immune response through caspase-1
Since inhibiting NLRP3 activation did not fully mitigate the G. vaginalis-mediated increases in IL-1β, caspase-1 or cell death in cervicovaginal epithelial cells, we sought to determine if inhibiting caspase-1 would fully block inflammasome activation. Using the irreversible caspase-1 inhibitor, Ac-YVAD-cmk 25, we determined that caspase-1 is an essential regulator of inflammasome activation by G. vaginalis. Inhibition of caspase-1 resulted in a significant reduction in the G. vaginalis-mediated increase in caspase-1 (p < 0.0001) (Fig. 4A, 4D, 4G), IL-1β (p < 0.0001) (Fig. 4B, 4E, 4H) and cell death (p < 0.0001) (LDH, Fig. 4C, 4F, 4I) in ectocervical, endocervical and vaginal cells after G. vaginalis but not L. crispatus exposure.
L. crispatus culture supernatants alter chromatin accessibility
As most of the enriched pathways associated with differentially expressed genes in all three cervicovaginal cell types exposed to L. crispatus culture supernatants were related to alterations in histone and transcriptional regulation (Fig. 2A, 2E-G), 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 26–28. While all three cell types demonstrated substantial modulation of gene expression after exposure to L. crispatus culture supernatants, the prominence of upregulated epigenomic-associated functional pathways in ectocervical cells pointed to potential changes in chromatin accessibility in this cell type.
The number and percentage of aligned/unaligned reads were equal across all three cell types (Supplemental Fig. 2). 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. 3). There was a strong correlation between normalized read counts between conditions for each cell type (Supplemental Fig. 4).
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. 7A). 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. 5). However, only in the ectocervical cells, the consensus peak regions demonstrated nearly two discrete clusters (Fig. 7B). Endocervical and vaginal cells did not exhibit these two clusters of accessible sites (Fig. 7B).
We then tested whether there were regions of differential accessibility between the treatment conditions in the cell type-specific consensus peak sites. 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. 7C). 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)29–31.
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 Encylopedia of DNA Elements (ENCODE) Project or open accessibility regions identified in different primary cancer specimens)32–34. Surprisingly, we discovered little overlap between differential accessibility in ectocervical cells and putative enhancer regions (Fig. 7D, 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 so we could only compare the ectocervical signatures 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. 7E, 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 8, Fig. 7F). 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 binding35. Gene-disease enrichment analysis of these transcription factors demonstrated a marked enrichment of multiple pathways related to neoplasms, endometriosis and infertility, all pathologies potentially associated with Lactobacillus-dominated microbiota (Fig. 7G) 36,37. 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. 6).