Adaptation of key bacterial vaginosis-associated bacteria to a medium simulating genital tract secretions: a transcriptomic analysis

Lúcia G. V. Sousa ¹ Angela França ¹ Vânia Pinheiro ¹ Christina Muzny ² Nuno Cerca ^1*

• ¹ University of Minho, Braga, Portugal
• ² University of Alabama at Birmingham, Birmingham, Alabama, United States

Approximately 30% of reproductive-age women around the world are affected by bacterial vaginosis (BV), which is the most common vaginal infection and the main cause of vaginal discharge (Peebles et al., 2019). BV is a polymicrobial vaginal dysbiosis with substantial data supporting the role of sexual transmission of BV-associated bacteria (BVAB) (Fethers et al., 2008; Muzny et al., 2022). In a case of BV, several strict and facultative anaerobic bacteria are increased in number in the vaginal microbiota replacing protective Lactobacillus species, which are the major colonizers present in a healthy vaginal microbiota (Rosca et al., 2020b; Chen et al., 2021). These anaerobic species are recognized to be living in a polymicrobial biofilm adhered to vaginal epithelial cells creating the characteristic clue cells (Swidsinski et al., 2024). Despite the identification of multiple bacterial species in cases of BV (Ceccarani et al., 2019), some may have a greater impact on BV development than others (Gilbert et al., 2019; Randis and Ratner, 2019; Castro et al., 2021). Previously, Gardnerella vaginalis, Fannyhessea vaginae, and Prevotella bivia were detected, in in vivo specimens from women with BV, early in the development of BV (Muzny et al., 2018), and further in vitro experiments showed important synergisms between these species suggesting an essential role in the initial stages of BV biofilm formation (Castro et al., 2021). Having this in mind, we recently investigated how these three bacterial species interact when growing on a nutrient-rich medium, which allowed for polymicrobial biofilm growth. Notably, we observed that hundreds of genes were differentially expressed when comparing triple-species with single-species biofilms (Sousa et al., 2024). However, the adaptation of bacterial species to different growth media results in alterations in the bacteria transcriptome and consequent variations in gene expression (Blair et al., 2013; Smith et al., 2018). As such, the consideration of which culture medium is most suitable to grow a specific species is crucial depending on the experiments that will be performed and the factors that need to be analyzed. Although the use of a nutrient-rich culture medium allows a better growth of bacterial species in laboratory settings, especially when growing fastidious microorganisms, using culture media that mimic in vivo conditions is preferable when trying to evaluate how microorganisms behave in in vivo-like conditions (Pouget et al., 2022). For the study of BV it is useful to have a medium that simulates the vaginal environment, and the medium simulating genital tract secretions (mGTS) contains several amino acids, as well as albumin, urea, and mucin found in the vaginal fluid (Tietz and Klein, 2018).Thus, to further investigate how the vaginal environment affects the transcriptome of the BV polymicrobial biofilm, we aimed to explore how transcriptome adaptations observed in the nutrient-rich medium would compare with adaptations to mGTS, which partially simulates in vivo vaginal conditions.2Materials and methodsBacterial strains and growth conditionsG. vaginalis ATCC 14018, F. vaginae ATCC BAA-55, and P. bivia ATCC 29303 were used in this study. The bacteria were kept in Brain Heart Infusion medium (Liofilchem, Roseto degli Abruzzi, Italy) supplemented with 23% of glycerol (98%, Panreac, Darmstadt, Germany) at -80ºC. Before each experiment, they were grown on Columbia Base Agar (Liofilchem) plates supplemented with 5% of defibrinated horse blood (CBA) (Thermo Fisher Scientific, Lenexa, KS) for 48 h at 37°C and in anaerobic conditions (AnaerocultTM A, Merck Millipore, Taufkirchen, Germany). New York City III (NYCIII) medium, a nutrient rich medium used for the growth of fastidious microorganisms, was prepared as previously described (Rosca et al., 2020a). mGTS medium, which simulates vaginal tract secretions, was adapted from previous studies (Liu et al., 2011; Stingley et al., 2014), and was prepared for the purpose of this study as previously described (Sousa et al., 2025).Biofilm formation Triple-species biofilms of G. vaginalis, F. vaginae, and P. bivia were formed on 24-well culture plates (Orange Scientific, Braine-l’Alleud, Belgium), in either NYCIII medium or mGTS, using the competitive model previously described (Rosca et al., 2022a). The biofilm competitive model considers that all species are present on the initiation of BV biofilm, and thus the species are incubated at the same time in the plates to form biofilm. For the biofilms formed in NYCIII medium, a suspension of each species was prepared and incubated at 37°C in anaerobic conditions for 24 h. For the biofilms grown in mGTS, since this species have limited growth in planktonic cultures in this medium, a suspension of each species was prepared from bacterial biomass growing on CBA plates on the same day of biofilm formation. Thereafter, for both protocols, the bacterial concentration was adjusted to 9107 CFU/mL by reading the optical density at 620 nm (Castro et al., 2021). The bacterial suspensions were then dispensed on the wells of the plate for a total volume of 1 mL, with a final concentration of 1107 CFU/mL for each species, and incubated for 48 h at 37°C under anaerobic conditions. For the characterization of the composition of triple-species biofilms by qPCR, the biofilm medium was removed and the biofilms washed once with 0.9% NaCl (Sigma, Germany), followed by mechanical detachment from the plates in 1 mL of 0.9% NaCl. The content of the wells was subsequently combined. Afterwards, genomic DNA was extracted and qPCR quantification was performed using specific primers for each bacterial species, as previously described (Rosca et al., 2022b). Using calibration curves previously designed (Lameira et al., 2024), the concentration of each species in the triple-species biofilm was determined. For RNA-seq experiments, the biofilm medium was removed, the biofilms were washed once with 1 phosphate-buffered saline (PBS) (Thermo Fisher Scientific) and then suspended in RNA protect Bacteria Reagent (Qiagen, MD, USA), diluted 2:1 in 1 PBS. The samples were then centrifuged at 5000 g for 10 min at room temperature. These experiments were repeated at least three times. RNA extractionTwelve triple-species biofilms formed in either NYCIII or mGTS were pooled to obtain samples with enough RNA concentration for further analysis. RNA extraction was performed using the RNeasy Mini Kit (Qiagen), according to the manufacturer’s instructions, with some minor modifications (França et al., 2012). First, cells were suspended in 600 μL of lysis buffer RLT and the volume was transferred to a tube with 0.1 mm zirconium beads (Merck, Darmstadt, Germany). Cells were lysed using the BeadBug 6 Microtube Homogenizer (Benchmark Scientific, NJ, USA) at maximum speed for 35 s. The lysis cycle was repeated four times and the samples were kept on ice for 5 min between cycles. Then, the samples were centrifuged, and the supernatant was recovered into a new tube. Ethanol at 70% (Thermo Fisher Scientific) was added in the same proportion (vol:vol) to the supernatant and the solution was transferred into an RNeasy Mini spin column. After the washing steps, total RNA was eluted in RNase-free water (Grisp) and the samples were treated with Turbo DNase (Invitrogen, Waltham, Massachusetts, USA) to degrade genomic DNA following the manufacturer’s instructions for rigorous protocols.cDNA library preparation and sequencingRNA quality was assessed using the Agilent 2100 Bioanalyzer (Agilent, CA, USA) and only samples with RNA quality indicators above 7 were used. RNA-seq libraries were prepared using Lexogen’s CORALLTM Total RNA-seq kit (Lexogen, Vienna, Austria) with 100 ng of total RNA from each sample. Before RNA-seq, rRNA was removed using the RiboCop for Bacteria (Mixed Bacterial Samples META) rRNA Depletion kit (Lexogen). Sequencing libraries were evaluated for quality on a Fragment Analyzer System (Agilent) and quantified with QubitTM dsDNA HS Assay Kit (Invitrogen). Sequencing data were generated using Illumina NextSeq 2000 Sequencing from single-end reads (SR100). FastQ files were generated via Illumina bcl2fastq2 (v.2.17.1.14). The quality of individual sequences was evaluated using FastQC software after adapter trimming with cutadapt software (1.18). RNA-sequencing data analysisFastQ files were then analyzed using CLC Genomics Workbench software (Qiagen, version 21.99). Quality trimming, including both quality scores and nucleotide ambiguity, was performed using the CLC genomics workbench default settings (Supplementary Table 1). Alignment of each species’ sequences was performed using G. vaginalis NCTC10287 (NCBI reference sequence: NZ_LR134385.1), F. vaginae FDAARGOS_934 (NCBI reference sequence: CP065631.1), and P. bivia DSM 20514 (NCBI reference sequence: NZ_JH660658.1; NZ_JH660659.1; NZ_JH660660.1) as reference genomes, also using the software default settings (Supplementary Table 2). Differential expression analysis was performed using Reads Per Kilobase per Million (RPKM) mapped fragments as the normalization strategy with triple-species biofilms grown in NYCIII as controls. Baggerley's test was applied to identify statistically significant alterations. Fold changes > 2 or < -2 and with a false discovery rate (FDR) p-value < 0.05 were considered significant and used for further bioinformatics analyses. Raw and analyzed datasets were deposited in the Gene Expression Omnibus database under the reference GSE279623.Functional annotationFunctional enrichment of differentially expressed genes was assessed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING, version 11.5) based on Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. Classes with FDR-adjusted p-values < 0.05 were considered for enrichment. REVIGO (version 1.8.1) was used for removing redundant GO terms. UniProt was used to find the homology of hypothetical proteins.Statistical analysisPrincipal Component Analysis (PCA) graphs and heatmaps were created using the CLC Genomics Workbench (version 21.99). All other figures and analyses were performed using GraphPad Prism version 8.2 (La Jolla, CA, USA). Statistical analysis was performed using two-way ANOVA with Tukey’s multiple comparisons test. Statistical differences were considered when p < 0.05.3Data descriptionTo study the influence of simulated in vivo conditions on the interactions of key BVAB, we analyzed the transcriptome of a triple-species biofilm composed of G. vaginalis, F. vaginae, and P. bivia, when biofilms were grown in a rich medium versus a medium mimicking genital tract secretions. Triple-species biofilms formed in mGTS were also analyzed by qPCR to determine each species' composition, revealing that the great majority of the biofilm was composed of G. vaginalis, followed by P. bivia, and F. vaginae, as shown in Supplementary Figure 1. Regarding RNA-seq data, the quality of the sequencing data was assessed by evaluating some of the summary parameters from the trimming and mapping steps as described at Supplementary Tables 3 and 4, respectively. The number of sequencing reads ranged between 13137570 and 20716698. The results of mapping to the reference genome sequence revealed a percentage of mapping between 1.82% and 75.76% for the NYCIII medium and between 0.63% and 67.48% for the mGTS medium. Overall, the sequences from P. bivia were the ones with the lowest percentages of mapping.The PCA plots, depicted in Supplementary Figures 2 to 6, revealed differences among the triplicates analyzed for each condition. Triplicates of biofilms formed in mGTS were more similar than those formed in NYCIII medium for G. vaginalis and P. bivia. In contrast, more differences were found among the biofilm triplicates in the mGTS condition for F. vaginae. These variations are likely a result from biological heterogeneity found in biofilms (Sousa et al., 2014). However, a clear separation from the triplicates in NYCIII and mGTS was observed for all the species, indicating that the two conditions caused alterations in the transcriptome of the three species.The distribution of gene expression was evaluated in both conditions for each species, as shown in Supplementary Figure 7. The box plots, which display the distribution of the RPKM values, revealed a similar distribution of the data in NYCIII and mGTS for the three species. The scatter plots, which represent the correlation between RPKM values in the NYCIII and mGTS conditions, showed that more genes were upregulated in the mGTS condition for G. vaginalis and F. vaginae (more genes displayed above the line), whereas for P. bivia more genes were downregulated in the same condition (more genes displayed below the line). Finally, heat maps depicting the normalized gene expression pattern in the three species across conditions and triplicates showed a clear alteration in gene expression from NYCIII to mGTS condition (Supplementary Figures 8 to 12). For differential gene expression analysis, biofilms formed in the NYCIII medium were used as the control. A total of 1315, 1202, and 2184 differentially expressed genes were found for G. vaginalis, F. vaginae, and P. bivia, respectively. MA and volcano plots (Figure 1) showed that the expression of the majority of the genes has a fold change between -32 and 32 (-5 < Log2 > 5). Moreover, for F. vaginae and P. bivia, a wide range of FDR p-values was found in comparison with G. vaginalis. For further analysis, only the genes whose fold change values were > 2 or < -2 and with FDR p-value < 0.05 were considered. This analysis resulted in 927, 492, and 166 upregulated genes for G. vaginalis, F. vaginae, and P. bivia, respectively. Moreover, 37, 104, and 490 genes were found to be downregulated in biofilms formed in mGTS for G. vaginalis, F. vaginae, and P. bivia, respectively. The top ten most upregulated and downregulated genes for each species and their respective functions are detailed in Table 1. To complete the analysis, a GO enrichment analysis was performed for the upregulated and downregulated genes of the three species and the results are reported in Supplementary Figures 13-15. The main categories of GO including biological processes, molecular functions, and cellular components were considered. For G. vaginalis, most of the upregulated genes were found associated with biological processes, mainly cellular and metabolic processes. Among the downregulated genes, enrichment was only found for molecular functions, namely for the active transmembrane transporter activity. For F. vaginae, more enrichment terms were found in the cellular components category, followed by molecular functions and biological process. The cellular anatomical entity and cellular process were the terms with the highest percentage of genes identified for both upregulated and downregulated genes, however, more terms of enrichment were found among the upregulated genes. For P. bivia, more GO terms were found for the category of biological process, followed by molecular functions and cellular components. GO enrichment was mostly observed within the downregulated genes which were mainly associated with metabolic processes. The upregulated genes were mainly associated with transport activity. The highest percentage of genes was observed for the term cellular anatomical entity for both upregulated and downregulated genes. Regarding KEGG pathways analysis (Supplementary Figure 16), different pathways were found enriched within the upregulated and downregulated genes for each of the three species. KEGG pathways among upregulated genes were only found for G. vaginalis with the identification of metabolic pathways, while for downregulated genes ABC transporters were identified. For F. vaginae and P. bivia, KEGG pathways were only found enriched within the downregulated genes and were mainly associated with metabolic pathways. Overall, this work highlights transcriptomic alterations in three key BVAB when growing in in vivo-simulating conditions compared to a nutrient-rich growth medium. It is important to note that, while these differences may be influenced by the growth medium, they may also reflect the shift in species abundance within the biofilm. Additionally, since the growth of some species, particularly F. vaginae, is very limited in the mGTS medium, the observed transcriptomic alterations for this species may result from stress responses caused by the lack of capacity to survive in this medium, rather than a response to the functional adaptation to the medium conditions. Further investigation is needed to elucidate the mechanisms involved in the transcriptomic alterations in these conditions and the contributions to the development of incident BV, including the specific interactions between key bacterial species.4Conflict of InterestCAM reports receiving grants to her institution from BioNTech, Lupin, Abbott, Visby, and Gilead Sciences, Inc. She also reports honorarium and/or consulting fees from Cepheid, BioNTech, biomerieux, Visby, Elsevier, Merck Manuals, UpToDate, Abbott, and Roche. The other authors have no disclosures. 5Author ContributionsCAM and NC designed the study. LS and VP performed the biofilm growth and qPCR and RNA extraction experiments. AF performed the analysis of RNA-seq data on CLC genomics. LS and NC drafted the manuscript. All authors participated in revising the manuscript. In summary, all authors contributed to this study and approved the submitted version of the manuscript.6FundingThis work was supported by the National Institute of Allergy and Infectious Diseases [R01AI146065-01A1 to CAM]. It was also partially funded by the Portuguese Foundation for Science and Technology (FCT), under the scope of the strategic funding of the unit [UIDB/04469/2020].7AcknowledgmentsLS acknowledges FCT for the financial support of individual Grant [2020.04912.BD] and AF acknowledge the financial support of FCT through the program DL 57/2016 – Norma transitória (DL57/2016/CP1377/CT0032).