Characteristics of the Cervicovaginal Microenvironment in Childbearing-Age Women with Different Degrees of Cervical Lesions and HR-HPV Positivity.

Persistent infection with high-risk human papillomavirus (HR-HPV) is the most important determinate in the development of cervical cancer, and cervical microecology can modulate cervical viral infection. However, few studies have been conducted on the microecological analysis of cervical diseases using strict physiological factors. This study investigated the characteristics and dynamics of cervical microecology in childbearing-age Chinese women with different degrees of HR-HPV-positive cervical lesions. A total of 168 subjects were selected according to the selection criteria, including healthy HPV-negative individuals (n = 29), HR-HPV-infected individuals (n = 29), low-grade squamous intraepithelial lesion individuals (LSIL, n = 32), high-grade squamous intraepithelial lesion individuals (HSIL, n = 40), and cervical cancer individuals (n = 38). We sampled cervical secretions from each subject and performed comparative analysis using the 16S rRNA sequencing method. Comparison analysis showed that Lactobacillus and Ignatzschineria were the dominant genera in the healthy group, while Gardnerella and Prevotella were more enriched in the disease groups. Based on the taxa composition, we roughly divided the development of cervical cancer into two phases: phase I was from healthy status to HR-HPV infection and LSIL; phase II was from LSIL to HSIL and cervical cancer. Different interactions among different genera were observed in different groups. Prevotella inhibited the abundance of Lactobacillus in the healthy group, while Prevotella inhabited the abundance of Gardnerella in the other groups. In the HR-HPV infection group, Ignatzschineria and Enterococcus showed a positive interaction but dissociated with the increase in cervical lesions, which might eventually lead to a continuous decrease in the abundances of Lactobacillus and Ignatzschineria.

Introduction

Cervical cancer is common cancer in women (Thun et al. 2010; Liang et al. 2019), and almost all cervical cancer cases are linked to human papillomavirus (HPV) infection, especially persistent infection with high-risk human papillomavirus (HR-HPV). Approximately 70–80% of women will develop at least one HPV infection in their lifetime (Ojesina et al. 2014; Bober et al. 2019), but most HPV infections can be cleared within 8–12 months. Fewer than 10% of patients develop persistent infection and even cervical cancer (Bober et al. 2019). Previous studies suggested that vaginal microecology can modulate HPV infection and is closely related to the progression of cervical intraepithelial neoplasia (Bober et al. 2019; Wiik et al. 2019). When histological changes occur in the cervical epithelium and vaginal mucosa, the abundance of Lactobacilli in vaginal microbes will decrease, and the rate of bacterial dysregulation will increase. These changes will also promote HPV infection and cervical lesions (Mitra et al. 2015; He et al. 2018). Therefore, vaginal microecology plays a vital role in developing cervical precancerous lesions and invasive cervical cancer and can be used as a biological indicator for cervical cancer tests (Curty et al. 2019). In addition, the risks of virus clearance and malignant transformation after HPV infection were also associated with vaginal flora (Zhang et al. 2018).

A large number of microorganisms gather in the female reproductive tract, including Gram-positive aerobic bacteria, Gram-negative aerobic bacteria, anaerobic bacteria, Candida spp., etc., and form a dynamic balance of mutual restraint with the host (Turnbaugh et al. 2007; Martin and Marrazzo 2016; Greenbaum et al. 2019). A healthy vaginal community can prevent vaginal infections, eliminate inflammation, and maintain the microecological balance (Green et al. 2015; Fettweis et al. 2019). When the balance of the vaginal community is destroyed, the incidence of bacterial vaginosis (BV) or vulvovaginal candidiasis (VVC) will increase (Green et al. 2015). Because of the special physiological structure of the vagina, the vagina is more suitable for the growth of anaerobic bacteria (McFall-Ngai et al. 2013), and thus BV is often incurable and prone to relapse, resulting in repeated vaginal infections. Previous studies showed significant differences in the proportion of BV-associated microorganisms among the HPV-negative group, HPV transient-positive group, and HPV persistent-positive group, suggesting that BV-associated microorganisms could affect women in both gynecological inflammation and tumors (Brotman et al. 2014; Bober et al. 2019). The unhealthy vaginal microbiome can promote the occurrence of inflammatory diseases and cervical lesions and participate in the pathogenesis of cervical tumors (Cancer Genome Atlas Research Network et al. 2017).

An increasing number of studies have focused on the vaginal microbiome in gynecological diseases, especially the correlation among vaginal microbiome, cervical lesions, and HR-HPV infection (Ojesina et al. 2014; Bober et al. 2019; Romero-Morelos et al. 2019; Zheng et al. 2019), and revealed that the diversities of the vaginal microbiome in cervical lesion patients with HPV infection were higher than those in healthy people (Silva et al. 2014; Bober et al. 2019). However, few studies have been conducted on the different degrees of cervical intraepithelial neoplasia before cancer appearance, and the influence of physiological factors on the microbiome is usually ignored.

This study collected cervical microecology samples from patients with HR-HPV infection and different degrees of cervical lesions under strict physiological conditions. 16S rRNA sequencing analysis was performed for all samples to detect the composition of the vaginal microbiome. To investigate the dominant microbiome in each group, a comparison analysis was performed between different groups. Our study may provide new insight into treating HR-HPV infection and blocking cervical lesions by speculating microecological regulation.

Experimental

Materials and Methods

Sample collection. In the present study, all samples were obtained from Chinese PLA General Hospital from January 2019 to December 2019. This study was approved by the Ethics Committee of Chinese PLA General Hospital (NO. S2018-221-01), and all subjects signed informed consent. A total of 168 subjects were involved, including 29 healthy subjects, 29 subjects with HR-HPV infection, 32 subjects with a low-grade squamous intraepithelial lesion (LSIL), 40 subjects with a high-grade squamous intraepithelial lesion (HSIL), and 38 subjects with cervical cancer. All subjects met the following criteria: 1) age from 30 to 50 with regular menstruation; 2) no sexual activity or vaginal medication within 3 days; 3) no inflammatory treatment within 1 week; and 4) body mass index (BMI) from 18.5 to 24. Cases were excluded with the following criteria: 1) diagnosis of cervical cancer before colposcopy; 2) vaginal disinfection before sampling; 3) pregnancy or breastfeeding; 4) pruritus vulvae and abnormal increase in vaginal secretions; and 5) treatment with oral sex hormone. The Kruskal-Wallis test was used for comparisons between different groups, which showed that the difference in the basic information of the subjects was not statistically significant (p > 0.05), and all selected subjects were comparable (Table I).

Table I

Analysis of essential information in each group.

Group	Cases number	Age	BMI
Group 1	29	40.08 ± 4.83	23.20 ± 3.47
Group 2	29	42.17 ± 5.18	22.00 ± 2.16
Group 3	32	40.63 ± 4.55	22.32 ± 2.93
Group 4	40	40.64 ± 5.57	22.37 ± 1.62
Group 5	38	42.43 ± 5.31	22.88 ± 1.96
K-W test (H)		5.336	7.640
K-W test (P)		0.255	0.106

Group 1 – healthy women; Group 2 – high-risk HPV infection; Group 3 – low-grade squamous intraepithelial lesion; Group 4 – high-grade squamous intraepithelial lesion; Group 5 – cervical cancer

DNA extraction and 16S rRNA sequencing. The cervical secretions of subjects were collected by sterile cotton swabs during gynecological examinations and placed in bacterial stabilizers (Xinyan Biological 50030) at –80°C. DNA was extracted using a Swab Genomic DNA kit (Kangwei Century) according to the instructions. DNA quality was evaluated using 1% agarose gel electrophoresis. To construct 16S rDNA sequencing libraries, 16S V3 and V4 regions were amplified using 341F (5’-CCTAYGGGRBGCASCAG-3’) and 806R (5’-GGACTCNNGGGTATCTAAT-3’) primers. The PCR products were mixed in equal amounts according to the concentration and purified using a 2% agarose gel by electrophoresis and gel cutting. Gel recovery was performed with a GeneJET Gel Recovery Kit (Thermo Scientific), and the products were used for library construction with an Ion Plus Fragment Library Kit 48 rxns (Thermo Fisher). Libraries were assessed by Qubit quantification and then sequenced on the Ion S5TMXL (Thermo Fisher) platform.

Data analysis and statistics. The low-quality reads were removed from raw sequencing data using Cutadapt (v.1.9.1) (Langille et al. 2013) and assigned to different samples according to the barcode. Then, the barcode and primer sequences were cut off from the sequencing reads (Martin 2011). To remove the insert sequences, clean reads were compared with the annotation database (Rognes et al. 2016). Tags were clustered to operational taxonomic units (OTUs) using Uparse (v.7.0.1001) (Haas et al. 2011) with 97% identity. Taxonomic assignment of OTU sequences was classified using the Mothur method and SSU rRNA database of SILVA132 (Edgar 2013) with a 0.8 ~ 1 confidence value as a cutoff.

To obtain normalized data for each sample, OTU sequences were aligned with MUSCLE (v.3.8.31) (Quast et al. 2013). To perform the alpha diversity analysis, QIIME (v.1.9.1) (Caporaso et al. 2010) was used to calculate the observed species index, Chao index, Shannon index, Simpson index, ACE index, and PD whole tree index, and R software (v.2.15.3) was used to draw dilution curves, rank abundance curves, and species accumulation curves. Bray-Curtis distance was also calculated using QIIME software (Version 1.9.1) to perform beta diversity analysis. The detection of different species between different groups was analyzed using a t-test. Based on the abundance of species, the Pearson correlation coefficient value (PCC) between each genus was calculated (cutoff value > 0.6 and abundance >0.005%), and the PCC value was used to perform network analysis using GraphViz-2.38.0.

Results

16S rRNA sequencing. A total of 168 samples were collected and divided into five groups: healthy group, HR-HPV infection group, LSIL group, HSIL group, and cervical cancer group. All samples were subjected to 16S rRNA sequencing, and an average of 77,213 original tags was obtained for each sample. After quality control, approximately 76,316 tags were obtained for each sample. There were 7,594 OTUs with 97% consistency clustered in all samples, with an average of 748 OTUs for each sample. In addition, 7,523 OTUs could be annotated by the SILVA132 database, including 49 phyla, 874 genera, and 723 species. The alpha diversity and beta diversity were carried out based on the species annotation, and the t-test was performed between different groups to find the species with significant differences (p-value < 0.05) in each group.

Community composition in different groups. The phyla Firmicutes and Actinobacteria were the most dominant in the five groups, followed by the phyla Proteobacteria and Bacteroidetes (Fig. 1A). Among the dominant phyla, the abundance of Firmicutes was the highest in the healthy group, and it fell in the HR-HPV infection group, followed by a slight increase in the LSIL group, and gradually decreased to its lowest level with the aggravation of lesions. The changing trend of Actinobacteria was opposite to that of Firmicutes. The abundance of Actinobacteria was the lowest in the healthy group, and then it increased in the HR-HPV infection group, slightly decreased in the LSIL group, and finally increased to the highest level in the cervical cancer group.

Fig. 1.

Microbiome communities of five groups.

The dominant genera mainly included Lactobacillus (Firmicutes) and Gardnerella (Actinobacteria), in which Lactobacillus was the most dominant genus. Similarly, the changing trend of Lactobacillus also contrasted with that of Gardnerella (Fig. 1C). With the development of lesions, the abundance of Lactobacillus slightly decreased at the beginning and then slightly increased and finally declined to the lowest level, while the abundance of Gardnerella first increased, then decreased, and finally increased to the highest level. The dominant species mainly included Lactobacillus iners, Lactobacillus jensenii, Prevotella bivia, Sneathia amnii, and Bacteroides fragilis. Among these dominant species, the relative abundance of L. iners was higher than that of other species; it was the highest in the healthy group and reached its lowest level in the cervical cancer group.

Comparison of microbiome diversity. The alpha diversity analysis showed that the healthy group had the highest value, and the HR-HPV infection group had the lowest value. The values of the intraepithelial lesion groups and the cervical cancer group were between the healthy and HR-HPV infection groups. The alpha diversity indices in the healthy group and the HR-HPV infection group were significantly different (p ≤ 0.05), except for the PD whole tree index. At the same time, there was no significant difference between the squamous intraepithelial lesion groups and the cervical cancer group (Fig. 2). The beta diversity value was the lowest in the healthy group, increased in the HR-HPV infection group, declined in the LSIL group, and finally reached the highest value with aggravation of the lesion (Fig. 3). Combined with the microbiome composition and diversity analysis, these five groups could be roughly divided into two phases. The first phase was from healthy to HR-HPV infection and LSIL, and the second phase was from LSIL to HSIL and cervical cancer.

Fig. 2.

Alpha diversity index of cervical microflora between groups. Group 1 – healthy group; Group 2 – HR-HPV infection group; Group 3 – LSIL group; Group 4 – HSIL group; Group 5 – cervical cancer group.

Fig. 3.

Beta diversity between five groups; A) beta diversity between groups; B) box plot based on weighted UniFrac beta diversity.

Phase I. From health to LSIL. A total of five significantly different species with three abundance change models were detected at the phylum level in phase I. Among them, the LDA values of two phyla were larger than four, including the phyla Firmicutes (LDA = 5.9) and Actinobacteria (LDA = 5.3). The abundance of Firmicutes was the highest in the healthy group and the lowest in the HR-HPV infection group, and Firmicutes was also the dominant phylum in the healthy and LSIL groups. The abundance of Actinobacteria was different from that of Firmicutes, and it was the highest in the HR-HPV infection group. Actinobacteria was a dominant phylum in the HR-HPV infection group (Table SI).

At the genus level, 75 significantly different species with five abundance change models were detected between different groups. The most abundant species continuously decreased, and seven genera showed LDA values larger than four. From the healthy group to the LSIL group, the abundance of Streptococcus (LDA = 4.2) was the highest in the healthy group and gradually decreased, and the abundances of Lactobacillus (LDA = 5.8) and Stenotrophomonas (LDA = 4.0) first increased and then decreased. The abundances of Ignatzschineria (LDA=4.9) and Atopobium (LDA=4.3) were consistent. The abundances of Gardnerella (LDA = 5.3) and Prevotella (LDA = 4.6) were the highest in the HR-HPV infection group, and they increased at the beginning and dropped with the infection (Fig. 4, Table SII). The microbial compositions of the healthy group and LSIL group were slightly different. The dominant genera of the healthy group were Lactobacillus (Firmicutes), Streptococcus (Firmicutes), and Stenotrophomonas (Proteobacteria). The genera Gardnerella (Actinobacteria) and Prevotella (Bacteroidetes) were the dominant genera in the HR-HPV infection group, while the dominant genera in the LSIL group were Lactobacillus (Firmicutes), Ignatzschineria (Proteobacteria), Atopobium (Actinobacteria), and Stenotrophomonas (Proteobacteria).

Fig. 4.

Abundance heatmap of different speciesat the genus level in phase 1.

Phase II: From LSIL to cervical cancer. Five significantly different species were detected with two abundance change models at the phylum level at phase II. Among them, the LDA values of the phyla Actinobacteria (LDA=5.5), Bacteroidetes (LDA=5.2), Firmicutes (LDA=5.8), and Proteobacteria (LDA=5.1) were larger than four. In phase II, the abundances of Actinobacteria and Bacteroidetes increased consistently and were the highest in the cervical cancer group. The abundances of Firmicutes and Proteobacteria steadily decreased and were the lowest in the cervical cancer group, showing that they had a negative correlation with the aggravation of the lesion (Table SIII).

A total of 63 significantly different species with four abundance change models were detected between other groups at the genus level. Among them, the abundances of most species gradually increased, and six genera showed LDA values larger than four. From the LSIL group to the cervical cancer group, the abundances of Gardnerella (LDA = 5.4), Prevotella (LDA = 4.8), and Faecalibacterium (LDA=4.1) increased consistently and were highest in the cervical cancer group. The abundances of Lactobacillus (LDA = 5.8) and Ignatzschineria (LDA = 5.0) steadily decreased and were lowest in the cervical cancer group. The abundance of Enterococcus (LDA = 4.1) initially increased and then decreased, and it reached its highest level in the HSIL group (Fig. 5, Table SIV). Some genera showed a positive correlation with the disease, which mainly included Gardnerella (Actinobacteria), Prevotella (Bacteroidetes), and Faecalibacterium (Firmicutes). In contrast, the negative correlation genera mainly included Lactobacillus (Firmicutes) and Ignatzschineria (Proteobacteria). The genus Enterococcus (Firmicutes) had a high abundance in the HSIL group, which might inhibit the lesions.

Fig. 5.

Abundance heatmap of different species at the genus level in phase 2.

Network analysis. Based on the dominant genera of the five groups, we found that the genera Lactobacillus and Ignatzschineria were the most connected to health, and the genera Gardnerella and Prevotella were the most related to the disease. In the healthy group, the abundance of Prevotella was negative with that of Lactobacillus, meaning that the abundance of Lactobacillus was limited, while the abundance of Prevotella was positive with that of Ignatzschineria. In addition, the genera Stenotrophomonas, Streptococcus, and Enterococcus were positive for each other to maintain microbiome balance. After infection, the abundances of Gardnerella and Atopobium were increased, and Prevotella suppressed the abundance of Lactobacillus, causing the abundance of Lactobacillus to drop sharply. Still, the abundance of Ignatzschineria increased and formed a positive interaction with Enterococcus, which plays an important role in maintaining the health of the microecology. The genus Gardnerella inhibited Lactobacillus, while the abundance of Gardnerella slightly dropped. Compared with the HR-HPV infection group, the abundance of Lactobacillus increased in the LSIL group because of the poor inhibition of Atopobium and Prevotella by Lactobacillus. Following the aggravation of lesions, the abundance of Gardnerella increased consistently and showed a stronger inhibitory effect on Lactobacillus, causing the abundance of Lactobacillus to decrease gradually. At the same time, the abundance of Enterococcus increased, which inhibited Ignatzschineria and made its abundance fall consistently. When the abundances of Lactobacillus and Ignatzschineria were minimized, the lesions gradually transformed to cervical cancer (Fig. 6).

Fig. 6.

Co-occurrence network diagram.

Discussion

Cervical cancer is the fourth most common cancer for women worldwide, but it is considered preventable because HR-HPV infection is necessary for cervical cancer (Borgdorff et al. 2014). The vaginal microenvironment plays an important role in HPV infection (Mitra et al. 2015). The vaginal environment is a comprehensive environment home to various microorganisms, in which the abundance of Lactobacillus is the highest (Salas-Jara et al. 2016). Histological changes in the cervical epithelium and vaginal mucosa could impact the microbial composition of the vagina (Zheng et al. 2019) and promote HPV infection and cervical lesions (Song et al. 2015; Di Paola et al. 2017). Vaginal microecology is closely related to cervical intraepithelial neoplasia (Mitra et al. 2015; 2016). Thus it could be a biological indicator for cervical cancer patients in clinical practice (Łaniewski et al. 2019). Here, we performed 16S rRNA sequencing of cervical secretions of patients with different degrees of cervical lesions and healthy controls and tried to identify the key factors of the cervical microecological system to provide a theoretical basis for preventing cervical cancer.

Most investigations of vaginal microecology have only focused on disease and ignored the impact of physiological factors (Klein et al. 2019). However, vaginal microecology can be impacted by human hormone fluctuations, microbial biomass, vaginal operations, etc. (Silva et al. 2014). In this study, strict standards of sample selection were used to exclude pathological factors as much as possible. We performed homogenization analysis on the age and BMI of subjects to show comparability among these groups. In addition, the sequencing depth was sufficient based on the sparse sample curve. Therefore, the fundamentals of our analysis were adequate and credible.

There was no significant difference between the disease groups and the healthy group based on the alpha diversity results, which was inconsistent with previous studies showing that community diversity was higher in the disease group than in the healthy group (Audirac-Chalifour et al. 2016; Salas-Jara et al. 2016; Kyrgiou et al. 2017). Our results might be significantly influenced by the sample types and cervical secretions, in which community diversity was usually high, which would reduce the difference between different groups. At the same time, the strict criteria for subject selection decreased the influence of physiological factors on the vaginal microbiome and increased the difficulty in sample collection (168 cases in a year) (Łaniewski et al. 2018). Furthermore, there may be a temporal relationship between the microbiome and the progression of the infection. In the study of (Audirac-Chalifour et al. 2016), the sample’s age ranged from 22 to 61, which ignored women’s menopausal status, which can impact the composition of the vaginal microbiome (Di Paola et al. 2017) and affect the accuracy of the results. Additionally, in contrast to previous studies, the degrees of cervical lesions were divided in more detail. We first reported that the community composition of the LSIL group was closer to that of the healthy group, which was consistent with the clinical prognosis of the disease and further confirmed that the cervical microbiome is related to disease prognosis.

Compared with the healthy group, the HR-HPV infection group had more diversities at the genus level and showed an obvious disordered vaginal microenvironment based on the beta diversity results consistent with previous reports (Klein et al. 2019). Among these different genera, the genera Lactobacillus, Gardnerella, Prevotella, and Streptococcus showed high abundance. Gardnerella vaginalis, as the only species of the Gardnerella genus, was reported to overgrow in patients with BV, and it was also related to health complications associated with BV (Morrill et al. 2020). In our study, the abundance of the Gardnerella genus was increased in the disease groups (Fig. 1C), which was consistent with previous studies. When the abundances of Lactobacillus and Streptococcus decreased, the abundances of other genera increased, which was also reported in previous studies (Oh et al. 2015).

Moreover, the healthy group had no difference from the LSIL group in weighted UniFrac beta analysis, and 60% of LSIL group patients could heal themselves in clinical recorders. However, we could not speculate whether these changes in the microbiome were the result of immune regulation or the effect of microecological regulation, and it should be discussed in the future by expanding the sample size. Our study provided reference data on the cervical microbiome and could play a key role in further studies for blocking HR-HPV infection and preventing cervical cancer.

Conclusion

In conclusion, the most related genera to the healthy group were Lactobacillus and Ignatzschineria, and the genera Gardnerella and Prevotella were most related to the disease groups. According to the microbiome’s composition, the development of cervical cancer could be roughly divided into two phases. The first phase was from health to HR-HPV infection and LSIL, and the second phase was from LSIL to HSIL and cervical cancer. At different phases, different genera had other interactions. Prevotella inhibited the abundance of Lactobacillus in the healthy group, while Prevotella inhibited the abundance of Gardnerella from the HR-HPV infection group to the disease groups. After infection, Ignatzschineria and Enterococcus showed a positive interaction but dissociated with the aggravation of the lesions, which eventually led to a consistent drop in the abundances of Lactobacillus and Ignatzschineria. The vaginal environment with low abundances of Lactobacillus and Ignatzschineria might promote the transformation of lesions to cancer.