Cervical Squamous Intraepithelial Lesions Are Associated with Differences in the Vaginal Microbiota of Mexican Women

ABSTRACT Cervical cancer is an important health concern worldwide and is one of the leading causes of death in Mexican women. Previous studies have shown changes in the female genital tract microbe community related to human papillomavirus (HPV) infection and cervical cancer; yet, this link remains unexplored in many human populations. This study evaluated the vaginal bacterial community among Mexican women with precancerous squamous intraepithelial lesions (SIL). We sequenced the V3 region of the 16S rRNA gene in cervical samples from 228 Mexican women, including 121 participants with SIL, most of which were HPV positive, and 107 healthy women without HPV infection or SIL. The presence of SIL was associated with changes in composition (beta diversity) and with a higher species richness (Chao1). A comparison of HPV-positive women with and without SIL showed that microbiota changes occurred even in the absence of SIL. Multivariate association with linear models (MaAsLin) analysis yielded independent associations between HPV infection and an increase in the relative abundance of Brachybacterium conglomeratum and Brevibacterium aureum as well as a decrease in two Lactobacillus iners operational taxonomic units (OTUs). We also identified a positive independent association between HPV-16, the most common HPV subtype linked to SIL, and Brachybacterium conglomeratum. Our work indicates that HPV infection leading to SIL is primarily associated with shifts in vaginal microbiota composition, some of which may be specific to this human population. IMPORTANCE Human papillomavirus (HPV) plays a critical role in cervical carcinogenesis but is not sufficient for cervical cancer development, indicating the involvement of other factors. The vaginal microbiota is an important factor in controlling infections caused by HPV, and, depending on its composition, it can modulate the microenvironment in vaginal mucosa against viral infections. Ethnic and sociodemographic factors influence differences in vaginal microbiome composition, which underlies the dysbiotic patterns linked to HPV infection and cervical cancer across different populations of women. Here, we provide evidence for associations between vaginal microbiota patterns and HPV infection linked to ethnic and sociodemographic factors. To our knowledge, this is the first report of the species Brevibacterium aureum and Brachybacterium conglomeratum linked to HPV infection or squamous intraepithelial lesions (SIL).

ABSTRACT Cervical cancer is an important health concern worldwide and is one of the leading causes of death in Mexican women. Previous studies have shown changes in the female genital tract microbe community related to human papillomavirus (HPV) infection and cervical cancer; yet, this link remains unexplored in many human populations. This study evaluated the vaginal bacterial community among Mexican women with precancerous squamous intraepithelial lesions (SIL). We sequenced the V3 region of the 16S rRNA gene in cervical samples from 228 Mexican women, including 121 participants with SIL, most of which were HPV positive, and 107 healthy women without HPV infection or SIL. The presence of SIL was associated with changes in composition (beta diversity) and with a higher species richness (Chao1). A comparison of HPV-positive women with and without SIL showed that microbiota changes occurred even in the absence of SIL. Multivariate association with linear models (MaAsLin) analysis yielded independent associations between HPV infection and an increase in the relative abundance of Brachybacterium conglomeratum and Brevibacterium aureum as well as a decrease in two Lactobacillus iners operational taxonomic units (OTUs). We also identified a positive independent association between HPV-16, the most common HPV subtype linked to SIL, and Brachybacterium conglomeratum. Our work indicates that HPV infection leading to SIL is primarily associated with shifts in vaginal microbiota composition, some of which may be specific to this human population. IMPORTANCE Human papillomavirus (HPV) plays a critical role in cervical carcinogenesis but is not sufficient for cervical cancer development, indicating the involvement of other factors. The vaginal microbiota is an important factor in controlling infections caused by HPV, and, depending on its composition, it can modulate the microenvironment in vaginal mucosa against viral infections. Ethnic and sociodemographic factors influence differences in vaginal microbiome composition, which underlies the dysbiotic patterns linked to HPV infection and cervical cancer across different populations of women. Here, we provide evidence for associations between vaginal microbiota patterns and HPV infection linked to ethnic and sociodemographic factors. To our knowledge, this is the first report of the species Brevibacterium aureum and Brachybacterium with HPV infection and high squamous intraepithelial lesion (HSIL) (women diagnosed with CIN2 or CIN3). All women were cancer free, but cases exhibited altered results in cytological, histological, and colposcopy analyses, and control, healthy women exhibited normal results.
By molecular analysis, the HPV subtype frequency was highest for HPV-16, HPV-58, and HPV-18 types, with 53.72%, 15.7%, and 9.1%, respectively, observed within cases. Within women without SIL, 32.7% were positive for HPV by PCR, and the distribution of HPV subtypes was similar to that observed in cases. Most of the women in both groups did not smoke (75.21% of cases versus 71.96% of controls), had regular menstrual periods (69.42% of cases versus 68.22% of controls), and most did not have intermenstrual bleeding (85.95% of cases versus 88.79% of controls). Statistically significant differences between groups were detected in relation to an active sexual life at the time of the study (73.55% of cases versus 92.53% of controls), use of contraceptives (66.12% of cases versus 52.34% of controls), and genital hygiene, recorded by the frequency of vaginal douching (81.82% of cases versus 52.34% of controls). More details of the characteristics of each group are described in Table 1  controls). There were also significant differences in the number of births and miscarriages, but these values are likely skewed by age differences between the two groups ( Table 2).
Associations between the vaginal microbiota SIL status. We determined the bacterial community by amplification and sequencing of the 16S rRNA gene (V3 region). The presence of SILs was associated with changes in bacterial alpha and beta diversity ( Fig. 1), with notable compositional differences at the family and genus level (Fig. 2). Beta diversity analysis, measured by principal-coordinate analysis (PCoA; Bray Curtis distance) (Fig. 1A) indicated that cervical SILs explain 1.4% of the variation in vaginal bacterial community structure (N = 228; Adonis P = 0.002). The presence of SILs was also associated with significantly higher species richness than an absence of SILs (Chao1; P = 1.15E204) (Fig. 1B) as well as a smaller but still significant increase in the Shannon index, which considers both richness and evenness (P = 0.041) (Fig. 1C), suggesting that the broadest diversity change is explained by bacterial community richness.
Within this study, we identified 35 participants who had not developed SILs, yet were HPV positive. This allowed us to compare the association between HPV infection and an absence of SILs. Interestingly, beta diversity, measured by principal-coordinate analysis, detected a very similar proportion of variation in composition to that detected for SIL status ( Fig. 2A). Moreover, there were no differences in beta or alpha diversities between HPV-positive samples with SILs and samples without SILs ( Fig. 2B to D), suggesting that the detected vaginal microbiota shifts are primarily associated with HPV infection, which, in many cases, subsequently leads to SILs. To further investigate the effect of SILs on the vaginal microbiota, we also compared samples from women with high-grade SILs with samples from women with low-grade SILs and did not detect any differences in alpha or beta diversity (Mann-Whitney P = 0.41; Adonis P = 0.11), further strengthening the link between HPV infection, rather than SILs, on the vaginal microbiota. Because of this, we further investigated the differences in vaginal microbiota explained by HPV infection. While results are similar to when data were stratified by SIL ( Fig. 1), the differences in species richness (Chao1) were more pronounced ( Fig. 3A to C). Relative abundance plots of the individual samples showed that while there is great interindividual variability, there are noticeable differences between women with and without HPV infection that align with the observed increase in species richness detected by Chao1 in women with SILs and/or HPV infection (Fig. 3D).
To study associations between specific taxa, SILs, and HPV status, given the importance of controlling for potential confounding variables that could explain or correlate with the detected associations between SILs or HPV status and the microbiota, we used multivariate association with linear models (MaAsLin). MaAsLin is a multivariate linear modeling tool with boosting that tests for associations between specific microbial taxa and continuous and/or Boolean metadata. This method reduces the total amount of correlations to be tested, therefore allowing for improvements in the robustness of the additive general linear models. With MaAsLin, we found significant independent associations between HPV-positive status and an increase in the relative abundance of Brachybacterium conglomeratum and Brevibacterium aureum as well as a decrease in two Lactobacillus iners operational taxonomic units (OTUs) ( Table 3). This indicates that no other variable, including SILs, explained the taxonomic differences influenced by HPV infection status. Interestingly, other independent associations were also detected between HPV subtypes or contraception use and several bacterial taxa (Table 3). For example, we detected a positive association between Brachybacterium conglomeratum and HPV 16, the most common HPV subtype associated with SILs, further suggesting the link between this bacterial species and HPV infection ( Table 3).

DISCUSSION
Several factors are known to play a role in cervical carcinogenesis, with HPV infection being one of the most important in the development of the disease (1). There are more than 100 types of HPV, of which at least 14 high-risk HPV types have been defined as carcinogenic (28). In this study, we found that over half of HPV infections are caused by the HPV-16 type, followed by HPV-58 and HPV-18, and all of them are considered high-risk HPVs worldwide (29). This predominance of the HPV-16 type was expected because it is generally accepted that HPV-16 is the major high-risk genotype in Mexico and in the world (30,31). We also found HPV-58 as the second most prevalent genotype in 15.7% of the cases, which is aligned with what has been reported in Asia (14.36 to 15.90%) (31).
Our study also revealed several other factors associated with SIL status, some of which reaffirm previously reported links (32). Factors positively associated with SILs included younger age, HPV infection, younger age of sexual debut, and the use of contraceptives, with the biggest difference explained by intrauterine device (IUD) use. In contrast, being sexually active at the time of the study and vaginal douching were linked to a reduced risk of SILs in this group of women.
Regarding contraceptive use, our results showed that an increased risk of SILs was linked to IUD use, which differs from that reported by Cortessis et al. (33) where they indicated that invasive cervical cancer can be approximately 30% less frequent in women who have used an IUD. Likewise, Agenjo et al. (34) described an inverse relationship between IUD use and cervical cancer risk, with women using IUDs reporting half the risk of developing this type of cancer. Our contrasting results, however, are in line with previous microbiome correlations with cervical cancer. We found significant correlations with IUD use and the increased abundance of Acinetobacter lwoffii, which has been previously reported in HPV-positive women (35). In addition, we detected an independent positive correlation with the use of IUDs and Fusobacterium sp. Fusobacterium has been studied as a possible diagnostic biomarker of cervical cancer because it is positively correlated with tumor differentiation (36). Furthermore, Fusobacteriaceae have been reported as the most abundant microorganisms in cervical carcinoma (37). Future work should study how IUD use may modulate the associations between the vaginal microbiota and SILs or cervical cancer.
While it is unclear why younger age was linked to SILs in our study, it is likely that it relates to the common age of onset of SILs, which occurs between 25 and 35 years of age (38,39). In contrast, healthy women would be less likely to visit the Instituto Mexicano del Seguro Social (IMSS) for a routine gynecological visit. Our microbiome results did not find any differences associated with age, suggesting that age did not confound our results. Several study variables linked sexual activity with SILs, including younger age of sexual debut. These and other related sexual behavioral factors have been previously linked with SILs, HPV infection, and cervical cancer risk (40). Interestingly, our study revealed that vaginal douching was linked to a reduced risk of SILs (odds ratio [OR] of 0.24; confidence interval [CI] of 0.13 to 0.44). Studies on cervical cancer and vaginal douching have reported positive, negative, and no associations (41). Although it is unlikely that SILs would lead to symptoms that would motivate genital douching, this practice is more common among women with other risk factors linked to sexually transmitted infections, which are a common cause of symptoms.
An important finding from our study is that the vaginal microbiota differences were primarily attributed to HPV infection (or subtype) and not SILs, indicating that infection itself may lead to changes in the vaginal microbial community. Among the predominant components of a healthy vaginal microbiome is the presence of Lactobacillus species, including L. crispatus, L. iners, L. jensenii, and L. gasseri (17,42), which results in reduced community diversity. Indeed, bacterial richness increases as Lactobacillus spp. levels are reduced in association with precursor lesions of cervical cancer (17) and with HPV infection itself (2,43,44). In support of this, our results showed higher species richness in cases as well as shifts in beta diversity. Compositional differences involved several taxa, including lactobacilli. Two L. iners OTUs were decreased in women with HPV infection. L. iners has been previously associated with a dysbiotic community and displays a series of characteristics that make this species different from other known vaginal lactobacilli (45)(46)(47). For instance, L. iners is a lower producer of D-lactic acid and induces interleukin-8 (IL-8) secretion, causing proinflammatory activity in the cervix, which may influence the progression of cervical intraepithelial neoplasia (15). In other studies, the dominance of L. iners and interactions with other vaginal anaerobic microorganisms alters the balance of the vaginal microbiota in association with cervical intraepithelial neoplasia (13). The most discriminant microbial differences between cases and controls involved Brevibacterium aureum and Brachybacterium conglomeratum (increased in cases). While these differences were very significant, these species were not uniformly present among either group (Fig. 3), suggesting that interindividual compositional differences may prevent the identification of microbiota species with biomarker potential for HPV infection or SILs. However, our study identified Brachybacterium conglomeratum as independently associated with HPV and with HPV-16, the most common subtype detected in our study. This finding prompts future investigation of the link between this bacterial species and SILs or cervical cancer risk associated with this specific HPV subtype and raises the possibility that microbiome links with HPV infection are subtype specific. To our knowledge, this is the first time this species has been linked to HPV infection. B. conglomeratum has not been readily reported in vaginal microbiome studies, which have surveyed mainly North American and European populations (48)(49)(50). This finding underlines the importance of considering ethnicity-and geography-driven differences in human microbiome studies, as dysbiotic patterns may be population specific.

MATERIALS AND METHODS
Study design. Healthy women and women infected with HPV, regardless of the degree of cervical squamous lesion, over 18 years of age (adult age in Mexico) and attending the Instituto Mexicano del Seguro Social (IMSS) in Mexico City were invited to participate as volunteers in this study. Written informed consent was obtained from all volunteers after providing detailed information about the study and its characteristics. The clinical research protocol and letter of informed consent were evaluated and approved by the Comité Local de Investigación y de Bioética de la División de Educación e Investigación Médica de la Unidad de Alta Especialidad Médica Pediatría del IMSS. All participants completed a study questionnaire that was used to obtain sociodemographic and risk factor information. Data were registered in a secured database for subsequent statistical analysis.
A total of 228 Mexican women over 21 years of age who attended the IMSS from December 2003 to July 2006 were enrolled in the study. These women were divided in two groups. The first group (controls) consisted of 107 healthy women with a mean age of 42.8 years (67.9) with three previous Papanicolaou (Pap) tests negative for HPV infection for three consecutive years (a fourth negative Pap result occurred at the time participants were invited to join the study) and diagnosed without SIL with normal cytology and colposcopy results by the treating gynecologist. The second group (cases) consisted of 121 patients with a mean age of 37.3 years (610.9) with different degrees of SILs and were positive for HPV infection based on cytology, histology, and colposcopy examination. This group included women diagnosed with cervical intraepithelial neoplasia from 1 to 3 (CIN1, CIN2, and CIN3) according to the Bethesda classification (51). Participants who had received treatment for vaginal or urinary infections, who were currently pregnant or up to 2 months postpartum, who had a history of hysterectomy, or who had a severe chronic disease were excluded from the study.
Besides SIL and HPV status, the following study variables were included in the analysis: age, HPV subtype, smoking (daily number of cigarettes consumed at time of study), regularity of menstrual period, intermenstrual bleeding, current sexual activity status, use and type of contraceptives (IUD [all copper], tubal ligation, hormonal contraceptives, or condoms), vaginal douching, age of sexual debut, number of sexual partners, number of pregnancies, and number of births and miscarriages.
Samples of vaginal exudate. Samples of vaginal exudate were taken by swabbing the mucosa using sterile Teflon swabs that were placed in a sterile 15-ml conical plastic tube with sterile 0.9% sodium chloride (Baxter physiological saline solution). The sample was kept at 220°C until its use for microbiome sequencing analysis.
Cervical DNA extraction and HPV detection and typing. Cervical DNA was extracted directly from the brushing done by the gynecologist before the application of acetic acid for the colposcopy study. The brush was stored in 1 ml of saline solution at 4°C for transport. DNA was extracted immediately from a first centrifugation at 8,000 rpm for 5 min, and the pellet was resuspended in 200 ml of extraction buffer (SDS 1%) and proteinase K (200 mg/ml), making a modification in the digestion of the sample (53°C for 5 h) and inactivating the reaction by boiling for 5 min (52). Once the sample was digested and inactivated, it was centrifuged at 12,000 rpm for 5 min, and the supernatant was removed and frozen at 220°C until use. HPV was detected via PCR using two sets of oligonucleotides, MY09/MY11 (53) and GP5/GP6 (54). Cycling conditions were used as previously described for the detection of HPV DNA in cervical cells (53)(54)(55). HPV DNA obtained from HeLa cell cultures containing 10 to 50 copies of the HPV-18 open reading frame (ORF) L1 was used as a positive control (56). All positive samples for HPV were subsequently typed with the HPVFast 2.0 kit (Pharma Gen SA, Madrid, Spain) according to the manufacturer's instructions.
16S rRNA gene sequencing. From vaginal DNA samples, the 16S rRNA gene was amplified by PCR in triplicate using bar-coded primer pairs flanking the V3 region as previously described (57). Each 50-ml sample of PCR mixture contained 22 ml of water, 25 ml of TopTaq master mix (Qiagen), 0.5 ml of each forward and reverse bar-coded primer, and 2 ml of template DNA. The PCR program consisted of an initial DNA denaturation step at 95°C (5 min), 25 cycles of DNA denaturation at 95°C (1 min), an annealing step at 50°C (1 min), an elongation step at 72°C (1 min), and a final elongation step at 72°C (7 min). Controls without template DNA were included to ensure that no contamination occurred. Amplicons were run on a 2% agarose gel to ensure adequate amplification. Amplicons displaying bands at ;160 base pairs (bp) were purified using an Illustra GFX PCR DNA purification kit. Purified samples were diluted 1:50 and quantified using PicoGreen (Invitrogen) in a Tecan M200 plate reader (excitation at 480 nm and emission at 520 nm).
For 16S rRNA gene sequencing, each PCR pool was analyzed on an Agilent Bioanalyzer using a highsensitivity double-stranded DNA (dsDNA) assay to determine approximate library fragment size and verify library integrity. Pooled library concentrations were determined using the TruSeq DNA sample preparation kit, version 2 (Illumina). Library pools were diluted to 4 nM and denatured into single strands using fresh 0.2 N NaOH. The final library loading concentration was 8 pM, with an additional PhiX spikein of 20%. Sequencing was performed using a Hi-Seq 2000 bidirectional Illumina sequencing and cluster kit, version 4 (Macrogen, Inc.). PCR products were visualized on E-gels, quantified using an Invitrogen Qubit with PicoGreen, and pooled at equal concentrations, according to a previous report (58).
Bioinformatic analysis of 16S rRNA gene sequences. All sequences were processed using mothur v.1.43.0 according to the standard operating procedure as previously described (59). Quality sequences were obtained by removing sequences with ambiguous bases, with a low-quality read length, and/or chimeras identified using chimera uchime. Quality sequences were aligned and compared to the SILVA bacterial reference version 132, and OTUs were generated using a dissimilarity cutoff of 0.03. The sequences were classified using the classify seqs command. There was a total of 1,099,273 reads in this sequencing run, which decreased to 1,072,800 after removing singletons and those samples with less than 1,000 reads with a median of 4,720 reads per sample (IQR of 1,175). A cutoff of 1,000 reads per sample was applied, and 63 were excluded from the analysis, yielding a final N of 228. Further filtering was applied to remove mitochondrion and chloroplast sequences and OTUs present in fewer than three samples. The OTU table was included as a supplemental data set in the supplemental material.
Statistical analysis. Differences in frequencies for categorical variables between cases and controls were evaluated using the chi-square statistic with Yates correction. Risk was estimated and expressed as an odds ratio (OR) and a 95% confidence interval (CI). For numerical variables, Mann-Whitney (nonparametric) or Student t tests (parametric) were used based on the D¨Agostino-Pearson normality test. We assessed vaginal microbial diversity and the relative abundance of bacterial taxa using Phyloseq (60) along with additional R-based computational tools in R Studio (R Studio, Boston, MA). Principal-coordinate analysis (PCoA) was conducted using Phyloseq and was statistically confirmed by permutational multivariate analysis of variance (PERMANOVA; Adonis test). The Chao1 and Shannon diversity indices were calculated using Phyloseq and statistically confirmed by Mann-Whitney or Kruskal-Wallis and Dunn tests if more than two groups were compared (GraphPad Prism software, version 5c, San Diego, CA). Multivariate association with linear models (MaAsLin; see reference 61) were used to calculate differentially abundant OTUs between the cases and controls, including several other study variables available from the metadata. The following covariates were fitted into the MaAsLin model based on previously reported associations with HPV infection or with microbiome shifts: SIL grade, HPV infection, HPV type, current smoking status, intermenstrual bleeding, sexual activity status, use of contraceptives, type of contraceptive, genital hygiene, age, age of sexual debut, number of sexual partners, number of sexual partners by age, number of pregnancies, number of births, and number of miscarriages.
Data availability. Sequences have been deposited to BioProject (accession number PRJNA766648).

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, XLS file, 0.3 MB.