Study specimens were obtained from participants in the Females Rising through Education, Support, and Health (FRESH) cohort, a prospective observational cohort study that enrolls 18- to 23-year-old, HIV-uninfected, non-pregnant women in Umlazi, South Africa. The cohort and its inclusion and exclusion criteria have been previously described in detail (Anahtar et al., 2015; Gosmann et al., 2017; Dong et al., 2018). The study protocol was approved by the Biomedical Research Ethics Committee of the University of KwaZulu-Natal (UKZN; Ethics Reference Number BF131/11) and the Massachusetts General Hospital Institutional Review Board (2012P001812/MGH). Participants provided informed consent.

All genital tract sampling was performed at the FRESH clinical research site in an exam room dedicated to pelvic examination. Sampling was performed by the same two trained research nurses throughout the study. Sampling was performed outside the menstrual period. Sample collection was performed under direct visualization by speculum exam. Lubricants were not introduced into the vagina prior to or during the speculum exam. The plastic single-use speculum was wet with tap water prior to insertion. Upon speculum insertion, swab samples (Puritan 6” Sterile Standard Foam Swab 126 w/Polystyrene Handle) were collected from the mid-vaginal wall, then from the ectocervical mucosa. A cervicovaginal lavage with 5 mL of sterile 0.9% saline (Adcock Ingram) was then performed, followed by collection of a cellular cytobrush (CooperSurgical, Medscand Cytobrush Plus, Catalog #C0005), which was rotated 360° in the cervical os. Cytobrushes were put into a 15 mL conical tube containing 5 mL of RPMI-1640 media supplemented with 10% v/v heat-inactivated fetal calf serum, 1% v/v Penicillin-Streptomycin-Amphotericin B 10K/10K/25ug mixture (Lonza Bioscience), L-glutamine (Lonza Bioscience, 2 mM final concentration), and HEPES (Lonza Bioscience, 10 mM final concentration). The tube was immediately placed on ice, and transferred to the processing lab at the HIV Pathogenesis Programme (HPP), Doris Duke Medical Research Institute, UKZN, for further processing (see details below). One cervical swab was used to make a microscope slide preparation for Gram stain analysis. The remaining cervical and vaginal swabs were placed in individual sterile vials on ice. A cervical swab was transferred, along with the prepared slide, to Neuberg Global Laboratories, Durban, South Africa, an accredited commercial laboratory diagnostics company, where polymerase chain reaction (PCR) testing was performed from the cervical swab for the sexually transmitted infections (STIs) Neisseria gonorrheae, Chlamydia trachomatis, Trichomonas vaginalis, and Mycoplasma genitalium. BV status was determined from the prepared slide by trained laboratory technologists at Neuberg Global Laboratories using the Gram-stain-based Nugent scoring method (Nugent et al., 1991). The remaining cervical and vaginal swabs were transferred to the HPP processing lab and stored at -80°C for subsequent sequence-based microbiome analysis.

Samples for host transcriptomic analysis were selected from participants who were STI-negative (chlamydia, gonorrhea, trichomonas, mycoplasma) at the time of sampling and who represented a distribution across the four CTs ( Table 1 and Appendix 2 ). Participants and samples were also selected based on maximal number of sorted cervical cells count, yielding comparable cell counts for each cell type across participants.

At the time of sampling, participants were asked about their current family planning methods and their last menstrual period (LMP). Women using no hormonal contraceptive were determined to be in the follicular phase of the menstrual cycle based on a self-reported LMP less than 14 days prior to sampling.

Total nucleic acids from cervical or vaginal swab samples were extracted with a phenol-chloroform method, which included a bead beating process to disrupt bacteria as previously described (Anahtar et al., 2015; Anahtar et al., 2016). The V4 region of the bacterial 16S rRNA gene was PCR-amplified following standard protocols (Caporaso et al., 2012; Anahtar et al., 2015; Hoang et al., 2020). Amplicons were pooled, purified, and prepared according to standard Illumina protocols, and single-end sequenced on an Illumina MiSeq using a v2 300-cycle sequencing kit with addition of custom Earth Microbiome Project sequencing primers (Caporaso et al., 2012).

Sequence demultiplexing was performed as previously described (Hoang et al., 2020). Demultiplexed sequences were then processed using dada2 version 1.6.0 (Callahan et al., 2016) with taxonomy assignments performed using the RDP training database supplemented by manual curation (Anahtar et al., 2015). Tables with amplicon sequence variant (ASV) taxonomy, per-sample ASV read counts, species read counts, and sequencing metadata are detailed in Supplementary File 1 . The denoised dada2 results with final taxonomic assignment were analyzed in R using phyloseq version 1.30.0 (McMurdie and Holmes, 2013) and custom R scripts (available in Supplementary File 2 .

Sample read counts after dada2 processing all exceeded 13,000 reads per sample. For additional analysis, we excluded ASVs that could not be defined at least to the taxonomic level of class, that were not represented by >10 reads in >2 different samples, or that had an abundance of <50 reads within the entire cohort. To explore variation between bacterial communities, counts from distinct ASVs belonging to the same species were merged ( Supplementary File 1 ). Bray-Curtis distances (β diversity) between all samples were calculated and between-sample variation was examined by performing principal coordinate analysis (PCoA) in phyloseq (McMurdie and Holmes, 2013), then plotting the ordinations in R using ggplot2 (Wickham, 2016) from the tidyverse package (Wickham et al., 2019).

Raw sequence read files for genital tract 16S rRNA-gene profiling are available in the NCBI Sequence Read Archive (SRA) under BioProject PRJNA738803 (all 29 BioSamples in the BioProject), plus BioSample SAMN19246318 from BioProject PRJNA730929. Custom R code with associated data files sufficient to reproduce the 16S-based microbiota analysis is available as a compressed supplementary file ( Supplementary File 2 ).

Cells were gently dislodged from cytobrushes without use of enzymatic digestion, and samples were processed and sorted fresh on the day of collection. We used fluorescence-activated cell sorting to isolate four specific live (Invitrogen blue viability dye) cell populations from the cytobrush samples on a BD FACSAria Fusion using BD FACSDiva software: antigen presenting cells (APCs; CD19- (clone H1B19, BD), CD45+ (H130, BD), CD66b- (G10F5, Biolegend), HLADR+ (G46-6, BD), CD3- (UCHT1, BD), and CD11c+ (B-ly6, BD) and/or CD14+ (M5E2, BD)), CD4+ T cells (CD19-, CD45+, CD66b-, CD3+, CD4+ (SK3, BD)), neutrophils (CD19-, CD45+, CD66b+), and epithelial cells (CD19-, CD45-, large by FSC-A, EpCAM+ (94C, Biolegend)) (See Appendix 1 for sorting strategy and Appendix 3 for cell counts). Sorted cell populations were stored in Trizol Reagent at -80C. We then performed bulk RNA-seq on each sorted cell population from each of the 30 participants, totaling 117 samples across the four cell types. RNA was extracted using a Chloroform method and quality was assessed by TapeStation. Bulk RNA library prep was performed on these RNA samples using the SmartSeq2 protocol (Picelli et al., 2014). Size and concentration of the libraries were assessed with TapeStation and Qubit.

Pooled cDNA was sequenced three times on an Illumina NexSeq using paired-end sequencing with a 75-cycle kit to achieve adequate sequencing depth. We used fastqc (0.11.5) and multiqc (1.5) for initial quality control to flag samples with poor sequencing quality (Ewels et al., 2016). We then used STAR (2.5.2b) for alignment and mapping, eliminating samples with less than 20% uniquely mapped reads (Dobin et al., 2013). HTSeq (0.9.1) was used to count the number of reads per gene (Anders et al., 2015). All samples, across all three sequencing runs, were visualized using PCA to ensure there was no data skewing attributable to any run. There was no evidence of batch effect across sequencing runs (data not shown), so the data from the three runs was merged for analysis. Next, the mitochondrial transcriptional content of each sample was assessed using EnsDb.Hsapiens.v75 in R. Because the percentage of mitochondrial transcript content was significantly different between cell types, different mitochondrial content cutoffs were used in different cell types to filter out samples with a high proportion of dead or dying cells. For epithelial cells, which had significantly higher percentage of mitochondrial transcript, we used a cutoff of <50%. For all other cell types, we used a cutoff of <30% of mitochondrial transcript as a cutoff for samples to be included in downstream analysis. The gene count table was filtered to retain only genes with ≥0.5 counts per million (CPM) in at least two samples. Genes that were not approved in HGNC were also excluded. After performing all QC filtering, the transcriptome of 107 (out of 117) samples was analyzed across 18248 genes.

Differential expression analysis was conducted in R using DESeq2 version 1.30.0 (Love et al., 2014). A ranked gene list was compiled from this DESeq2 analysis and analyzed using GSEA v4.1.0 (Mootha et al., 2003; Subramanian et al., 2005). The ranked gene list was assessed for enriched gene sets using the hallmark gene list version 7.4 (Liberzon et al., 2015) and gene ontology (GO) biological processes version 7.4 (Ashburner et al., 2000; Gene Ontology, 2021). An FDR-adjusted p-value of 0.1 was used as the statistical significance cutoff. R scripts for host RNAseq analysis are available in Supplementary File 3 .

All samples and data were collected from participants enrolled in the FRESH (Females Rising through Education, Support, and Health) study in Umlazi, South Africa, who were between the ages of 18 and 23 years at the time of study enrollment and were HIV-uninfected, non-pregnant, and did not have sexually transmitted infections (STIs) at the time of sampling. We selected samples from a subset of 30 participants for paired transcriptome-microbiome analysis to ensure distribution of participants across the four vaginal microbiota community types as determined by bacterial 16S rRNA gene sequencing ( Figure 1 and Table 1 ). We have previously established a classification system that divides the FGT microbial communities into four cervicotypes (CTs) (Gosmann et al., 2017; Munoz et al., 2021). CT1 is dominated by Lactobacillus crispatus and CT2 is dominated by Lactobacillus iners. CT3 and CT4 both have higher diversity. In CT3, Gardnerella is the most abundant taxon, while CT4 is dominated by other taxa, typically with a significant prevalence of Prevotella species ( Figure 1 ). Prior work has shown that women with CT3 and CT4 communities exhibit higher levels of FGT inflammation, as measured by both inflammatory cytokines in cervicovaginal lavage and by an increased frequency of inflammatory cells, including HIV target cells (activated CCR5+ CD4+ T cells), in cervical cytobrush samples (Anahtar et al., 2015). CT3 and CT4 are also associated with an increased risk of acquiring HIV (Gosmann et al., 2017).

The 30 participants examined in this analysis included eight in CT1, five in CT2, five in CT3, and twelve in CT4 ( Figure 1A ). We performed principal coordinates analysis (PCoA) based on Bray-Curtis distances to investigate the range of variation between participant vaginal microbial communities. The first two PCoA axes represented 42.4% of the total variance in the population ( Figure 1B ). Axis 1 fully segregated Lactobacillus-dominant (CT1 and CT2) communities from non-Lactobacillus-dominant (CT3 and CT4) communities, while Axis 2 fully segregated CT1 from CT2. The assignment of the four CTs ( Figure 1B ) was consistent with previous published work that studied larger sample sizes (Anahtar et al., 2015; Gosmann et al., 2017). The clear separation of FGT microbial communities based on Lactobacillus dominance or depletion (Axis 1, Figure 1B ) highlighted the relevance of investigating host transcriptional differences between CT1/2 and CT3/4.

We divided the cohort based on FGT Lactobacillus dominance, comparing Lactobacillus dominant (CT1/2) to Lactobacillus deficient (CT3/4) communities. No factors related to sexual practices differed between groups, including the number of days since last having sex, the number of sexual encounters over the previous 30 days, and the number of sexual partners over the past 30 days ( Table 1 ). There was no significant difference in reported use of condoms, vaginal drying agents, or contraceptives; however, there was a non-significant trend toward increased use of the progestin-based long-acting reversible contraceptives (LARCs) DMPA, Nuristerate, and Implanon, in the CT1/2 group. The only significant difference was in days since last reported menstrual period (LMP) (p=0.009): women with CT1/2 had a median of 45 days elapsed since LMP, longer than a regular menstrual cycle and likely reflecting the amenorrhea associated with progestin-only contraceptive methods. Women with CT3/4 had a median of 16 days since LMP, and a range within that of regular menstrual cycles. No sampling occurred during menstruation.

We sorted four cell lineages from each of participant: epithelial cells, APCs, CD4+ T cells, and neutrophils ( Appendix 1 , Appendix 2 ). These populations were selected based on their role in genital immune surveillance, barrier function, and HIV acquisition. We then performed bulk-RNA sequencing on each sorted cell population from each participant. In unsupervised clustering, the CTs did not cluster together (data not shown). Although we knew we were underpowered in this analysis, we proceeded with a hypothesis-generating comparison between Lactobacillus-dominant (CT1/2) and Lactobacillus-deficient (CT3/4) CTs. There were no significant differences in cell count per sample between CT groups for any cell type ( Appendix 3 ) that might impact the differential gene expression analysis.

APCs demonstrated the greatest number of significantly differentially expressed genes (n=111) and the most inflammatory genes and pathways associated with non-Lactobacillus-dominant CTs ( Figure 2A and Appendix 4 ; Appendix 5 ). The APC genes most strongly associated with CT3/4 include those involved in the interferon response (e.g., IFIT2, IFIT3, IFIT5, OASL), class I major histocompatibility complex (e.g., HLA-S, HLA-G, HLA-H, HLA-L), class II major histocompatibility complex (e.g., HLA-DRB5), genes related to innate immunity (e.g., MEFV, PELI1, PTGS2, IL1RN), and genes promoting T-cell mediated adaptive immunity (e.g., TNFSF9). Genes related to regulation of neoplastic cell growth were also upregulated in the CT3/4 group (e.g., MYCL, DDIT4, PLAC8, and PMAIP1); the genes upregulated in CT3/4 predominantly suppress cell growth and/or induce apoptosis. Hierarchical clustering of samples based on genes with significant differences in expression nearly completely segregated CT1/2 from CT3/4 samples ( Appendix 4A ).

In accordance with these observations, pathways identified as significantly enriched through gene set enrichment analysis (GSEA) using hallmark gene sets in CT3/4 included multiple inflammatory pathways (e.g., TNFα signaling via NFkB, IFNγ response, inflammatory response, interferon alpha response, complement, allograft rejection, IL2 STAT5 signaling, and IL6 JAK STAT3 signaling) ( Figure 2B ). Pathways related to neoplastic suppression were also highlighted in this group (e.g., apoptosis, p53 pathway), as well as pro-proliferative pathways (KRAS signaling up, MTORC1 signaling, PI3K AKT MTOR signaling). Conversely, genes significantly enriched in CT1/2 included innate immune response genes (BPIFB1, PELI3), fractalkine (CX3CL1), genes related to mucosal epithelial integrity (TFF3, CLDN10, MUC5B, LAMB1, RIPK4), genes related to neoplastic invasion (CTSK, KLK11, UCA1), and few related to tumor suppression (S100A14). GSEA pathway analysis also highlighted proliferative signatures, such as progress through cell cycle checkpoints and Myc signaling, in CT1/2 as compared to CT3/4 ( Figure 2B ). Additional GSEA analysis performed with the gene ontology (GO) biological processes gene sets further supported APCs as drivers of dysbiosis-associated FGT inflammation, with 607 pathways upregulated in CT3/4, many of which were related to inflammatory signaling ( Appendix 5A ).

Epithelial cells represent another possible mediator of FGT inflammation. These cells exhibited very few differentially expressed genes based on microbiota composition ( Figure 2C ). No GSEA hallmark-derived gene sets appeared to have any relevance for the inflammatory signal seen grossly in the setting of CT3/4 ( Figure 2D ). A few inflammatory pathways were significantly upregulated in CT3/4 based on GO terms, but they were more related to innate inflammation, such as response to fungus, eosinophil migration, monocyte chemotaxis, and granulocyte chemotaxis, and all had relatively low normalized enrichment scores (NES 1.96-2.36) ( Appendix 5B ).

CD4+ T cells are thought to be a final common pathway of inflammation that increases risk of HIV acquisition (Haase, 2005; Haase, 2010). In this analysis, CD4+ T cells had the second greatest number of significantly differentially expressed genes after APCs (15) ( Figure 2E ). However, hierarchical clustering of samples based on differentially expressed genes was unable to segregate samples from different CTs ( Appendix 4C ), unlike the clustering observed for APCs ( Appendix 4A ). Furthermore, although a non-Lactobacillus-dominant FGT microbiome is known to be associated with inflammation, inflammatory pathways were surprisingly upregulated in CD4+ T cells from women in CT1/2 as compared to CT3/4 ( Figure 2F ). GO term analysis also showed some inflammatory pathways in CT1/2, such as TLR2 signaling and myeloid leukocyte mediated immunity ( Appendix 5C ), possibly suggesting a host-mediated role in curating the FGT microbiome away from BV-associated bacteria (Mares et al., 2008). The most strongly upregulated pathways in CT3/4 were related to enhanced translation capabilities (e.g. cotranslational targeting to membrane (NES 2.79), establishment of protein localization to ER (NES 2.46), translation initiation (NES 2.28)), which might relate to a CD4+ T cell activation state previously observed with dysbiosis (Anahtar et al., 2015).

Neutrophils showed no significantly differentially expressed genes between CT groups ( Appendix 4D ). The hallmark gene set “TNFα signaling via NFkB” was significantly enriched in CT3/4, but IFNα response was significantly enriched in CT1/2 ( Appendix 4E ), and both had relatively low normalized enrichment scores. GO term analysis similarly showed a small number of inflammatory pathways upregulated in CT3/4 ( Appendix 5D ).

Overall, these results indicate a strong inflammatory signal through APCs in the setting of Lactobacillus-depleted vaginal microbial communities. While other cell types may contribute to the inflammatory milieu, the prominent influence of APCs in helping to orchestrate a dysbiosis-associated inflammatory milieu is clear.

In order to explore whether hormonal effects also contribute to the transcriptional differences we observed between women with different vaginal microbiota, we compared samples from 7 women using DMPA (i.e., a high-progestin, low-estrogen state) with those from 6 women in the follicular phase of the menstrual cycle (i.e., a low-progestin, high-estrogen state; Appendix 6 ). Each hormone group was split relatively evenly by CT, resulting in 2-4 participants from each hormonal state in each CT group ( Appendix 6 ). No demographic variables examined were significantly different between groups ( Appendix 6 ). Due to the small sample size of an integrated analysis, we did not explicitly address the modifying effect of hormones on microbiota-associated inflammation.

Numerous transcriptional differences were observed across all cell types by hormonal state ( Figure 3 and Appendix 8 ). In APCs, epithelial cells, and CD4+ T cells, the pattern of differential gene expression nearly completely segregated women in the follicular phase from DMPA users ( Figures 3C, F, I ), suggesting a true and clear difference in transcriptional programs between groups.

In APCs, the hallmark gene set “inflammatory response” was significantly enriched ( Figure 3B ). In epithelial cells, one of the most highly differentially expressed gene across DMPA users was CCL3L3, which is a chemotactic cytokine for CCR5+ CD4+ T cells ( Figures 3D, F ). GO term analysis highlighted epithelial barrier differences in DMPA users (cornification, keratinization), while decreased leukocyte mediated immunity was observed in DMPA users, though with a low NES ( Appendix 9B ). Only five genes were significantly differentially expressed in CD4+ T cells ( Figures 3G, I ), but these genes were able, again, to perfectly segregate samples by hormone state.

These results demonstrate that DMPA is associated with significant changes in cellular pathways within FGT mucosal cell populations. Previous observations regarding DMPA-associated CD4+ T cell infiltrates may be partially explained by inflammation in APCs, chemotactic cytokines secreted by epithelial cells, and/or reduced epithelial barrier integrity.