From November 2023 to April 2024, 123 Caucasian women of reproductive age were enrolled in the study. The participants were volunteer students attending degree courses at the University of Bologna, Italy. At enrolment, the exclusion criteria were: (i) antibiotic use in the month prior to sampling; (ii) use of vaginal douches or topical agents in the last two weeks; (iii) age < 18 years; (iv) pregnancy; (v) menstruating at the time of sampling; (vi) HIV infection; (vii) presence of chronic conditions (e.g., diabetes, autoimmune diseases, malignancies); (viii) drug addiction or heavy smoking (>15 cigarettes/day). Moreover, women with urogenital infections due to sexually transmitted pathogens (i.e., Chlamydia trachomatis, Neisseria gonorrhoeae, Trichomonas vaginalis, Mycoplasma genitalium) were further excluded when diagnostic test results were available.

Several demographic (e.g., BMI, educational level, marital status), health-related (e.g., presence of Candida spp., smoking habits) and antibiotic-related behavioral (e.g., individual practice and awareness about antibiotics) data were collected (see specific paragraphs on exposures).

After obtaining written informed consent from all participants, they underwent two self-collected vaginal samplings at the same time point. The first swab (E-swab, Copan, Brescia, Italy) was used for diagnostic tests to exclude the presence of sexually-transmitted infections (STIs) by NAAT (Nucleic Acid Amplification Test, Alinity m STI Assay, Abbott Molecular Inc, Des Plaines, IL, USA). The second was collected with a sterile cotton bud, re-suspended in 1 mL of sterile saline, and stored at -80 °C until use.

The study was conducted in accordance with the Declaration of Helsinki and approved by the Institutional Review Board (IRB) of the University of Bologna, protocol number 0122421.

Nucleic acids were extracted from vaginal swabs by means of the DNeasy Blood & Tissue Kit (QIAGEN GmbH, Hilden, Germany) according to the manufacturer’s instructions. From the DNA eluate, each sample was screened for the presence of selected resistance genes and specific transposons that confer resistance to different classes of antibiotics by means of end-point PCR assays. The following genes were investigated: erm(A), erm(B) and erm(F) conferring resistance to macrolides; tet(M), tet(O), tet(W), and tet(Q), conferring resistance to tetracyclines; blaOXA-2, blaTEM, blaZ, blaSHV and blaCTX-M associated with β-lactam resistance; and qnrA and qepA linked to quinolone resistance. Additionally, the presence of two transposons associated with AR, namely tet(M)-Tn916 and tet(Q)-rteA, was also assessed. Detailed information on primers used and PCR conditions are displayed in the Supplementary Materials ( Supplementary Table S1 ). These results were dichotomized as 0=absence of antibiotic resistance and 1=presence of antibiotic resistance.

An overall composite score for the ARGs was developed by summing the presence/absence of all the genes detected in at least one sample. The total possible score ranged from 0 to 10, with higher scores indicating greater antibiotic resistance. The score was also categorized into tertiles based on its distribution: low (≤3), medium (4-5), and high (≥6).

The hypervariable V3-V4 regions of the bacterial 16S rRNA gene were amplified from genomic DNA extracted from vaginal swabs. The PCR conditions and primer sequences were obtained from the Illumina 16S Sample Preparation Guide (https://support.illumina.com/documents/documentation/chemistry_documentation/16s/16s-metagenomic-library-prep-guide-15044223-b.pdf) (Illumina, San Diego, CA, USA), with the primers originally described in Klindworth et al., 2013. Final indexed libraries were prepared by equimolar (4 nmol/L) pooling, denaturation, and dilution to 6 pmol/L before loading on a MiSeq flow cell (Illumina) for a 2 × 300 bp paired-end run. 16S rRNA sequences were processed according to the same methods as described in Severgnini et al., 2021, using the latest release of the SILVA database (v. 138, Quast et al., 2013) for classification.

Depth of sequencing was set to the lowest sequenced sample (n=6,627 reads), in order to compensate for the sequencing unevenness of the samples and to provide a consistent minimum amount for the downstream analysis, carried out through the QIIME pipeline (version 1.9.0, Caporaso et al., 2010). Alpha-diversity was estimated according to Chao1, Shannon Index, Observed Species, and the Faith’s phylogenetic tree diversity metrics (“PD whole tree”), whereas beta-diversity was evaluated using unweighted and weighted UniFrac distances (Lozupone et al., 2011). Since all alpha-diversity indices performed similarly, we focused on the PD whole tree for the univariate and multivariate analyses (see below).

Lactobacillus species-level characterization was performed as in Severgnini et al., 2022 by BLAST-aligning all reads belonging to the Lactobacillaceae family to a custom reference database made up of all available reference sequences in the NIH-NCBI database (https://ftp.ncbi.nlm.nih.gov/genomes/GENOME_REPORTS/prokaryotes.txt) of 17 species commonly found in the vaginal environment, including a total of 3,392 genomes. In case of multiple matches with the same confidence, the taxonomy was reset to “Unclassified Lactobacillus” at the genus level.

Community-state types (CST) of the vaginal microbial communities were determined from the taxonomic profiles using VALENCIA, a nearest centroid-based tool that classifies samples into 5 major CST according to the similarity to a set of about 13,000 reference microbial profiles (France et al, 2020).

In order to highlight common patterns of abundance among taxa, a co-abundance groups (CAGs) analysis was performed as previously described (Claesson et al, 2012). Briefly, only the bacterial taxa present at >0.5% of abundance in at least 10% of the samples (n≥11) were selected, in order to exclude minor and transient contributors of the gut microbiota. This resulted in a subset of 21 taxa on the whole dataset. The co-abundance between each pair of taxa was evaluated by calculating Spearman’s correlation coefficient and displayed in a heatmap, hierarchically clustered using Pearson’s correlation metric and average linkage. Only associations having a Benjamini-Hochberg adjusted p<0.05 were used to build the hierarchical clustering. Permutational multivariate analysis of variance (P-MANOVA) was used to assess that the computed CAGs were significantly different from each other, using p<0.05 as a threshold.

For the present analysis, exposures were selected based on theoretical knowledge because they are crucial for contextualizing the results, as they could potentially influence both antibiotic use and health outcomes (Morsli et al., 2024). These factors included the following demographic, health- and antibiotic-related behavioral data of the participants: age (continuous), educational level (categorized as high school vs university degree or higher), marital status (married/cohabiting vs single/separated), body mass index (BMI, calculated as weight in kg divided by height in cm squared), contraceptive use, smoking habits (non-smokers vs smokers), presence/absence of Candida spp., adherence to a healthy diet such as the Mediterranean diet (measured calculating the MEDILITE score, Sofi et al., 2017). In addition, we investigated individual practice and awareness about antibiotics and AR, and participants’ family of origin approach and behavior toward antibiotics, through ad hoc questions, to identify possible misuse of these drugs. An adapted version of the “Antibiotic Resistance: Multi-Country Public Awareness Survey,” developed by the World Health Organization (WHO, 2015), was employed. This revised survey consisted of 25 items organized into several sections. The first section contained two questions regarding the amount of antibiotics taken in the last year and month, along with four questions focused on best practices for antibiotic use where participants’ individual scores were derived from the sum of their answers. The second section included nine statements aimed at assessing participants’ understanding of antibiotic resistance, with responses recorded as True or False. Individual scores were calculated by summing the responses. The final section comprised ten statements evaluating participants’ beliefs and attitudes towards antibiotic use and resistance, utilizing a 5-point Likert scale (1 = strongly disagree to 5 = strongly agree). The total scores from these items were used as individual scores to gauge “awareness,” which refers to the level of self-awareness and understanding regarding strategies to combat AR. Eight ad hoc questions were included to investigate participants’ family background and attitudes towards antibiotics, with the aim to identify family models which used antibiotics appropriately (e.g., avoiding self-prescribing and adhering to prescribed dosages). In order to better explore the variable related to family antibiotic improper use, we categorized the total score in tertiles based on the distribution as low (≤10), medium (11–17) and high (≥18).

Sample characteristics were summarized using mean and standard deviation (SD) for continuous variables and counts and percentages for categorical variables.

Statistical evaluation of alpha-diversity indices for the microbiota profiles was performed by non-parametric Monte Carlo-based tests. Beta-diversity differences were assessed by a permutation test with pseudo F-ratios using the “adonis” function from R package “vegan” (version 2.0-10, https://cran.r-project.org/package=vegan). Relative abundances were analyzed first setting all “0” to a value of 10^–6 in order to compensate for missing data due to the normalized ASV table (where the “0” can be due to insufficient depth of coverage for assessing the ASV presence) and, then, performing a centered log-ratio (CLR) transformation in order to deal with the compositional nature of the data; finally, a Kruskal-Wallis test, followed by a Dunn post hoc evaluation was performed. Correlation between presence of specific ARGs and bacterial taxa abundances was performed by point-biserial correlation, whereas correlation between the total ARG score and bacterial taxa abundances was based on the Spearman’s correlation coefficient (with a p-value ≤ 0.05 considered as significant).

Two regression models were performed to identify factors associated with antibiotic-resistance genes. Binary logistic regression models were constructed to calculate odds ratios (ORs) and 95% confidence intervals (CIs) to estimate the association between exposures of interest and the presence of each of the antibiotic resistance genes.

Linear regression models were used to calculate the exponentiated β-coefficients with 95% CIs to examine the associations between exposures and the continuous antibiotic-resistance gene score.

We first performed univariate analyses including all the following variables in the model separately: age, educational level, marital status, smoking habits, BMI, presence of Candida, family history of inappropriate antibiotic use, antibiotics taken in the past year, past compliance with antibiotic use, awareness of antibiotic resistance, adherence to MEDILITE, and biodiversity index (PD whole tree). Multivariate models were performed, including factors associated with a p-value ≤ 0.10 in univariate analyses. All analyses were performed with STATA (version 15.0, StataCorp LP, College Station, Texas, USA) and Matlab (v 2008b, Natick, MA, USA) and the functions from the Fathom Toolbox (Jones, 2017). Two-sided p-values of 0.05 and 0.10 were considered statistically significant or borderline significant, respectively.

Raw sequencing data for this project are available in NCBI Short-Read Archive (SRA) under accession number PRJNA1188525 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA1188525).

A total of 123 participants were initially recruited. After excluding individuals with missing vaginal microbiota data (n = 7) and/or incomplete antibiotic behavior questionnaires (n = 14), 105 women were included in the final analysis. The mean age of participants was 21.6 ± 2.6 years (range: 19–30), and the mean BMI was 21.5 ± 2.8 kg/m² (range: 16.2–30.1).

Vaginal samples were categorized into the five main CSTs according to their microbial composition: CST I (41.9%, n = 44), CST II (5.7%, n = 6), CST III (30.5%, n = 32), CST IV (17.1%, n = 18), and CST V (4.8%, n = 5).

Many participants reported some level of inappropriate antibiotic use within their family of origin, with 32.4% indicating a high tendency toward misuse. Regarding individual antibiotic consumption, 42.9% of participants reported not having taken antibiotics in the previous year.

Overall, participants demonstrated high awareness of antibiotic resistance (mean score: 42 ± 4), along with strong adherence to recommended antibiotic use practices (mean score: 9 ± 2). Baseline characteristics of the sample, including behavioral factors and antibiotic-related variables, are presented in Table 1 .

Regarding ARGs distribution, erm(F) was the most detected one, with 72.2% of positivity (80 cases), followed by: tet(M) (74.3%, 78 cases), erm(B) (68.6%,72 cases), erm(A) (66.7%, 70 cases), tet(W) (65.7%, 69 cases), blaZ (36.2%, 38 cases), tet(Q) (35.2%, 37 cases), and blaTEM (16.2%, 17 cases). Furthermore, the less frequently detected genes were blaOXA-2 and tet(O), both with 2.9% of positivity (3 cases), qnrA with 1.9% (2 cases) and qepA with only 1% (1 case). No samples tested positive for blaSHV or blaCTX-M.

Positivity for specific ARGs was frequently associated with the presence of the corresponding mobile genetic elements. Specifically, tet(Q)-rteA was identified in 100% of the tet(Q)-positive samples (37/37), and tet(M)-Tn916 was detected in 96.2% of the tet(M)-positive samples (75/78).

We initially performed univariate analyses to explore potential associations between selected factors and the presence of ARGs (data not shown). Variables with a p-value ≤ 0.10 were subsequently included in multivariate models to assess their independent association with the outcomes of interest.

Several associations remained statistically significant in the multivariate analyses ( Table 2 ). In detail, we found that the microbial alpha-diversity, estimated by the Faith’s phylogenetic distance metric (PD whole tree), was positively associated with tet(M), tet(O), tet(Q), erm(A), and blaZ, with odds ratios ranging from 1.57 to 2.22. Faith’s PD was also strongly associated with the total ARG score (Exp(β-coefficient) 1.82, 95%CI 1.53-2.17), whereas borderline significance was observed with blaTEM (p-value ≤ 0.10).

Contraceptive use reduced the odds of tet(M) by 70% compared to non-use. Smoking increased the odds of erm(B) by 3.47 times and the Exp(β-coefficient) of the total score increased by 2.09 times. A higher BMI was positively associated with erm(A) and a higher total ARG score, with a borderline positive association with tet(O). Diet quality, as assessed by the MEDILITE score, was inversely associated with tet(O) detection, suggesting a potential protective effect of adherence to a Mediterranean diet. The presence of Candida spp. was significantly associated with a higher total ARG score (OR = 3.85).

Regarding antibiotic behavior, a higher degree of inappropriate antibiotic use in the family was associated with a slight increase in the odds of erm(A) and a borderline association with blaOXA-2. Antibiotic use in the past year was positively associated with tet(M). Past compliance of antibiotic usage mainly reduced the odds of tet(M) and blaZ. Moreover, a borderline inverse association with erm(B) and the total score was also noticed (with odds ratios ranging from 0.47 to 0.80). Awareness about antibiotic resistance decreased the odds of tet(Q) by 14%.

The multivariate analysis indicated that the alpha-diversity of the vaginal microbiota (according to the PD whole tree metric) was positively associated with the presence of several ARGs. At the same time, the microbial profiles between samples positive or negative to several ARGs were found to be statistically different. In particular, we confirmed the different microbial composition of samples positive/negative for erm(A), erm(B), tet(M), tet(O), tet(Q), and blaZ (unweighted UniFrac distance, p ≤ 0.008, adonis test of pseudo F-ratios). Moreover, also other ARGs (i.e.: erm(F) and tet(W)) had different microbial profiles (p=0.001 and p=0.002, respectively, unweighted UniFrac), whereas blaTEM, borderline associated to biodiversity, showed only a trend towards separation (p=0.191 unweighted UniFrac). Finally, microbial profiles for blaOXA-2, qnrA and qepA were not significantly different, due to the very low fraction of samples positives for these ARGs (n=3, n=2 and n=1, respectively) ( Supplementary Figure S1 in the Supplementary Materials ).

Interestingly, significant changes in both alpha- and beta-diversity were found when analyzing the vaginal microbiota and the number of ARGs per sample. In fact, a clear trend towards an increase in alpha-diversity with increasing ARGs was evident, both considering the total ARG score and its categorization into tertiles. Moreover, the beta-diversity analysis showed a significant separation among the samples according to the ARG score tertiles (p<0.008 and p<0.024 for all pairwise comparisons, unweighted and weighted UniFrac distances, respectively) and also a pattern along principal coordinate 1 moving from the right to the left of the x-axis with increasing total ARG scores ( Figure 1 ).

Considering the bacterial taxa involved in the vaginal microbiota, the abundance of some genera such as Prevotella, Porphyromonas, Dialister, Finegoldia, and Peptoniphilus were significantly (p<0.05, Kruskal-Wallis and Dunn’s post hoc pairwise tests) increased in both the “medium” and the “high” categorized ARG score tertiles as compared to “low”, whereas Gardnerella, Streptococcus, Megasphaera, Haemophilus, Anaerococcus, and Fusobacterium increased in “high” score samples only ( Supplementary Figure S2 in the Supplementary Materials ). Species from the Lactobacillus genus did not show a significant decrease with increasing total ARG scores but, as a whole, they accounted for about 90% of the relative abundance in “low” samples, while shifting to about 70% in “medium” and to about 55% in “high” ones ( Figure 2 ).

Co-abundance analysis highlighted the presence of 4 groups (CAGs) with a similar pattern of abundance of their members over all the samples ( Supplementary Figure S3 in the Supplementary Materials ), which were labelled according to the most representative taxon/taxa ( Supplementary Table S2 in the Supplementary Materials ). Notably, while L. crispatus, L. jensenii and L. gasseri all belonged to the same CAG, L. iners stood apart together with Ureaplasma; the two other CAGs were those of taxa more frequently associated to vaginal dysbiosis, comprising Gardnerella, Prevotella, Atopobium, Dialister and Megasphaera in one and Streptococcus, Porphyromonas, Anaerococcus, and Finegoldia in the other.

Correlating together the relative abundances of the microbial CAGs and the presence/absence of specific ARGs over all the samples ( Figure 3 ) highlighted some interesting features, such as the positive association between several genes related to resistance against macrolides (i.e.: erm(A), erm(F)) and tetracycline (i.e.: tet(M), tet(Q)) and the Gardnerella-Prevotella CAG, whereas the Streptococcus CAG showed some weak positive correlations, in particular, to tet(W) and blaZ. On the other hand, the L. crispatus/jenesenii/gasseri CAG was negatively associated with all ARGs, particularly to tet(M) and tet(Q), whereas the L. iners CAG was negatively correlated to blaZ resistance. Notably, a positive correlation to quite all the ARGs was observed for the “Other” CAG, which groups together all the low-abundance taxa (i.e.: relative abundance >0.5% in <10% of the samples). At the same time, the Gardnerella-Prevotella, the Streptococcus and the “Other” CAGs were positively correlated to the total ARG score, whereas the L. crispatus/jenesenii/gasseri CAG was negatively correlated.

When evaluating single taxa, we found that, among the Gardnerella-Prevotella CAG, Prevotella and Dialister showed the higher number of ARGs (erm(A), erm(B), erm(F), tet(M), tet(Q)), whereas Atopobium, Gardnerella, and Megasphaera all showed a correlation to tet(Q); all these genera were positively correlated to the total ARG score. At the same time, in the Streptococcus CAG, several positive correlations were found for Finegoldia (i.e.: erm(A), erm(F), tet(M), tet(Q)), Porphyromonas and Anaerococcus (i.e.: erm(A), erm(B), erm(F), tet(M), tet(Q), tet(W), tet(O), blaZ); again, all of them were directly correlated to the total ARG score. Conversely, in the L. crispatus/jenesenii/gasseri CAG, L. crispatus was negatively associated with all ARGs, particularly to tet(M) and tet(Q), while the other lactobacilli in this CAG showed a similar behavior, except with a positive association between L. gasseri and erm(A), and between unclassified members of the Lactobacillus genus and qnrA. Finally, among the bacteria comprised in the “Other” CAG, in particular, Sneathia showed a positive correlation to qepA and tet(Q), whereas Escherichia-Shigella and Haemophilus were correlated to blaTEM and blaZ genes.