Front. Microbiol. Frontiers in Microbiology Front. Microbiol. 1664-302X Frontiers Media S.A. 10.3389/fmicb.2023.1264768 Microbiology Original Research The vulvar microbiome in lichen sclerosus and high-grade intraepithelial lesions Pagan Lisa1 2 Huisman Bertine W.1 2 van der Wurff Michelle3 Naafs Rosanne G. C.1 Schuren Frank H. J.3 Sanders Ingrid M. J. G.4 Smits Wiep Klaas4 5 Zwittink Romy D.4 5 Burggraaf Jacobus1 6 Rissmann Robert1 6 7 * Piek Jurgen M. J.8 Henderickx Jannie G. E.4 5 † van Poelgeest Mariëtte I. E.1 2 † 1Centre for Human Drug Research, Leiden, Netherlands 2Department of Gynaecology and Obstetrics, Leiden University Medical Center, Leiden, Netherlands 3Netherlands Organisation for Applied Scientific Research (TNO), Zeist, Netherlands 4Department of Medical Microbiology, Leiden University Center of Infectious Diseases (LU-CID), Leiden University Medical Center, Leiden, Netherlands 5Department of Medical Microbiology, Center for Microbiome Analyses and Therapeutics, Leiden University Medical Center, Leiden, Netherlands 6Leiden Amsterdam Center for Drug Research, Leiden University, Leiden, Netherlands 7Department of Dermatology, Leiden University Medical Center, Leiden, Netherlands 8Department of Obstetrics and Gynaecology, Catharina Cancer Institute, Eindhoven, Netherlands

Edited by: Atte Johannes Von Wright, University of Eastern Finland, Finland

Reviewed by: Amy Rebecca Sapkota, University of Maryland, United States; Karl Kashofer, Medical University of Graz, Austria

 *Correspondence: Robert Rissmann, rrissmann@chdr.nl

†These authors share last authorship

 29 11 2023 2023 14 1264768 24 07 2023 01 11 2023 Copyright © 2023 Pagan, Huisman, van der Wurff, Naafs, Schuren, Sanders, Smits, Zwittink, Burggraaf, Rissmann, Piek, Henderickx and van Poelgeest. 2023 Pagan, Huisman, van der Wurff, Naafs, Schuren, Sanders, Smits, Zwittink, Burggraaf, Rissmann, Piek, Henderickx and van Poelgeest

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

 Background

The role of the vulvar microbiome in the development of (pre)malignant vulvar disease is scarcely investigated. The aim of this exploratory study was to analyze vulvar microbiome composition in lichen sclerosus (LS) and vulvar high-grade squamous intraepithelial lesions (HSIL) compared to healthy controls.

 Methods

Women with vulvar lichen sclerosus (n = 10), HSIL (n = 5) and healthy controls (n = 10) were included. Swabs were collected from the vulva, vagina and anal region for microbiome characterization by metagenomic shotgun sequencing. Both lesional and non-lesional sites were examined. Biophysical assessments included trans-epidermal water loss for evaluation of the vulvar skin barrier function and vulvar and vaginal pH measurements.

 Results

Healthy vulvar skin resembled vaginal, anal and skin-like microbiome composition, including the genera Prevotella, Lactobacillus, Gardnerella, Staphylococcus, Cutibacterium, and Corynebacterium. Significant differences were observed in diversity between vulvar skin of healthy controls and LS patients. Compared to the healthy vulvar skin, vulvar microbiome composition of both LS and vulvar HSIL patients was characterized by significantly higher proportions of, respectively, Papillomaviridae (p = 0.045) and Alphapapillomavirus (p = 0.002). In contrast, the Prevotella genus (p = 0.031) and Bacteroidales orders (p = 0.038) were significantly less abundant in LS, as was the Actinobacteria class (p = 0.040) in vulvar HSIL. While bacteria and viruses were most abundant, fungal and archaeal taxa were scarcely observed. Trans-epidermal water loss was higher in vulvar HSIL compared to healthy vulvar skin (p = 0.043).

 Conclusion

This study is the first to examine the vulvar microbiome through metagenomic shotgun sequencing in LS and HSIL patients. Diseased vulvar skin presents a distinct signature compared to healthy vulvar skin with respect to bacterial and viral fractions of the microbiome. Key findings include the presence of papillomaviruses in LS as well as in vulvar HSIL, although LS is generally considered an HPV-independent risk factor for vulvar dysplasia. This exploratory study provides clues to the etiology of vulvar premalignancies and may act as a steppingstone for expanding the knowledge on potential drivers of disease progression.

 vulvar microbiome vaginal microbiome lichen sclerosus vulvar HSIL HPV vulvar cancer section-at-acceptance Antimicrobials, Resistance and Chemotherapy Introduction

Vulvar squamous cell carcinoma (VSCC) occurs in 1–2 per 100,000 women and has an increasing incidence with age (Hacker et al., 2012). VSCC is preceded by human papillomavirus (HPV)-related or HPV-independent precursor lesions (van de Nieuwenhof et al., 2009; de Sanjosé et al., 2013; Nitecki and Feltmate, 2018). Predominantly HPV type 16 and 18 can cause vulvar high-grade squamous intraepithelial lesions (HSIL), a premalignant condition responsible for approximately 20% of VSCC (de Sanjosé et al., 2013). These high-risk HPV (hrHPV) types are also notorious drivers of cervical dysplasia and carcinoma (Cubie, 2013). Vulvar lichen sclerosus (LS) is a chronic inflammatory condition that may promote development of differentiated vulvar intraepithelial neoplasia (dVIN), a premalignancy preceding the remaining 80% of all VSCC (Thuijs et al., 2021; Voss et al., 2021). dVIN has low disease incidence and poses a diagnostic challenge for both clinicians and pathologists. The central symptom in LS is pruritus, while scarring and anatomical changes of the labia minora and clitoral phimosis occur in severe cases. The etiology of LS remains debated, with indications of a genetic or autoimmune root cause (van de Nieuwenhof et al., 2009; Thuijs et al., 2021). A considerable amount of taboo is associated with vulvar disease, resulting in delays in clinical recognition and treatment, adding to substantial physical, sexual and psychological morbidity (Lockhart et al., 2013; Grimm et al., 2016; Pérez-López and Vieira-Baptista, 2017; Mullen et al., 2019).

Changes in microbiome composition have been associated with several disease conditions, including cancer (Feng et al., 2015; Byrd et al., 2018). The research field started off with recognition of single malignancy-driving micro-organisms (such as HPV for cervical cancer) and has expanded to include association of microbiome composition patterns to cancer development (Picardo et al., 2019). In vaginal diseases, higher grades of HPV-driven cervical dysplasia are correlated with a non-Lactobacillus dominated vaginal microbiome composition (Mitra et al., 2015; Piyathilake et al., 2016; Norenhag et al., 2020). Moreover, the presence of Sneathia spp. and Mycoplasma spp. has been correlated to co-infection with hrHPV types (Lee et al., 2013; Mitra et al., 2015; Audirac-Chalifour et al., 2016; Adebamowo et al., 2017; Klein et al., 2018; Łaniewski et al., 2019). These findings can expand our understanding of cervical diseases and serve as biomarkers and potential targets for drug development to treat cervical (pre)malignant disease.

In contrast to the vaginal microbiome, the vulvar microbiome and its role in the development of (pre)malignant vulvar disease is scarcely investigated and current knowledge is inconclusive (Pagan et al., 2021). The bacterial genera and species described on the healthy vulva include the genera Lactobacillus, Corynebacterium, Staphylococcus, and Prevotella, suggesting transfer from vaginal, cutaneous and intestinal origin. Current literature lacks longitudinal, case–controlled studies and elucidation of microbiome components other than bacteria, such as viruses, fungi and archaea. Therefore, the aim of this study was to describe and compare the vulvar microbiome composition by metagenomic shotgun sequencing in patients with lichen sclerosus and vulvar HSIL compared to healthy controls. Biophysical properties of the vulvar skin were additionally explored as the microbiome composition can influence the skin barrier function (Harris-Tryon and Grice, 2022).

 Methods Trial design and subjects

This study was part of an observational clinical trial to identify biomarkers for premalignant vulvar disease to increase the understanding of the etiology of VSCC (Huisman et al., 2023). The trial was conducted at the Centre of Human Drug Research in Leiden (the Netherlands) from February 2021 to October 2021. The Declaration of Helsinki was the guiding principle for trial execution and the study was approved by the medical ethics committee “Medisch-Ethische Toetsingscommissie Leiden Den Haag Delft” prior to initiation. Written informed consent from all participants was obtained prior to any study-related procedure.

Criteria for inclusion were women aged between 25 and 95 years with BMI <30 kg/m2. Patients with vulvar HSIL were required to have ≥1 demarcated lesion(s) ≥15 mm in diameter with confirmed histologic diagnosis. Patients with lichen sclerosus were considered eligible in case of a clinical and/or histological diagnosis of LS. Exclusion criteria were significant concomitant diseases, pregnancy, immunodeficiency, sexually transmitted disease, AIDS or hepatitis. Furthermore, individuals with other dermatological conditions in the genital area were excluded.

Lifestyle restrictions were incorporated to standardize vulvar conditions for microbiome sampling. A ≥ 28-day wash-out applied for systemic immunomodulatory medication and antibiotic use (topical or systemic). Wash-out for topical medication on the vulvar area was ≥14 days. Participants were instructed to refrain from sexual intercourse, vigorous exercise, applying vulvar products and shaving, waxing or other depilatory treatments at least 24 h before every visit. Additionally, they were instructed not to wash the vulvar area from midnight onwards on every visit.

In total, 10 healthy volunteers, 5 patients with HSIL and 10 patients with LS were enrolled in this observational study. Anal, vaginal and vulvar samples were obtained at a single time point for vulvar HSIL and LS patients, as well as healthy volunteers.

 Microbiome characterization Sample collection

The vulvar microbiome was characterized in LS and HSIL patients as well as in the healthy volunteers. Microbiome samples were obtained using a 1 mL DNA/RNA Shield™ Collection Tube with Swab (Zymo Research, Freiburg, Germany). The pre-wetted swab with saline was rubbed along the vulvar skin for 30 s while slowly rotating the swab. Sampling locations on the vulva (i.e., labia minora, labia majora or perineum) depended on the location of the vulvar lesions (Supplementary Table S1). Healthy-appearing vulvar skin contralateral to the HSIL lesion was selected as non-lesional HSIL reference site where possible. Of note for LS, the non-hirsute vulvar skin should be considered affected, even if the skin appears healthy. Therefore, non-lesional LS sites were selected on the distal side of the labia majora toward the groin. Reference samples of the vaginal microbiome were obtained by introducing a dry swab mid-vaginally and once rotating 360 degrees along the vaginal wall, without touching the vulvar area upon introduction or removal. Another reference sample of the anal microbiome was obtained by rubbing a dry swab along the anus five times. Samples were stored in DNA/RNA shield at −80°C until DNA extraction was performed.

 DNA extraction

DNA of vulvar, vaginal and anal swabs was extracted with the Quick-DNA Fecal/Soil Microbe Miniprep kit (Zymo Research; D6010). During DNA extraction, positive controls (D6300 ZymoBIOMICS Microbial Community standard, Zymo Research, United States) and negative controls (empty tubes) were included. In short, 600 μL BashingBead Buffer was added to the swabs and processed with Precellys 24 Homogenizer (Bertin Technologies) at 5,500 rpm for three rounds of 60 s each. Subsequently, samples were centrifuged (10,00 RCF; 1 min), 800 μL supernatant was transferred to a Zymo-Spin III-F Filter and centrifuged again (8,000 RCF; 1 min). Next, 1,200 μL Genomic Lysis Buffer was added to the filtrate. Of the resulting mixture, 800 μL was transferred to a Zymo-Spin IIC Column and centrifuged (10,000 RCF; 1 min). After discarding the flowthrough, 200 μL of DNA Pre-Wash Buffer was added to the Zymo-Spin IIC Column and centrifuged (10,000 RCF; 1 min). Five hundred microliter g-DNA Wash Buffer was added to the Zymo-Spin IIC Column and centrifuged (10,000 RCF; 1 min) after which 50 μL DNA Elution Buffer was added and centrifuged (10,000 RCF; 30 s). The eluted DNA was transferred over the same column and centrifuged (10,000 RCF; 30 s). Lastly, the eluted DNA was transferred to a prepared Zymo-Spin III-HRC Filter and centrifuged (16,000 RCF; 3 min). The resulting DNA was quantified with a Qubit 4 fluorometer (Invitrogen). In total, 89 out 90 clinical samples yielded measurable concentrations of DNA.

 Metagenomic shotgun sequencing

DNA and additional positive sequencing controls were analyzed using metagenomic shotgun sequencing by GenomeScan (Leiden, the Netherlands). Upon sample entry, quality of samples was assessed by the Fragment Analyzer (Advanced Analytical Technologies) according to GenomeScan’s protocol. Given their low biomass, 73/97 (75%) samples passed entry quality control of >30 pg/mL (of which 15/25 (60%) anal, 25/25 (100%) vaginal, 29/39 (74%) vulvar samples, 0/4 (0%) negative controls and 4/4 (100%) positive controls), yet all samples passed library preparation quality control and were included for metagenomic shotgun sequencing.

Sequencing libraries were prepared using Illumina’s DNA PCR-Free Prep kit and checked on quality with the Fragment Analyzer. Libraries were sequenced with the Illumina NovaSeq6000 platform to a target depth of 3.3 million reads per sample.

 Biophysical assessments Trans-epidermal water loss

Measurement of the trans-epidermal water loss (TEWL) determines the skin barrier function in a non-invasive manner (AquaFlux AF200 System, Biox, London, United Kingdom). The measurements were performed under standard environmental conditions and patients were acclimatized with removed clothing for ≥15 min before initiation of the measurements. All TEWL measurement conditions were constant during the study, with mean probe temperatures of 24.3°C and average humidity of 39.2%. A measurement was considered valid at the settling of the flux curves at a final steady level, as described previously (Niemeyer-van der Kolk et al., 2022).

 Vulvar and vaginal pH analysis

Vulvar skin pH was determined using an electronic pH probe (Skin-pH-Meter PH905, Courage and Khazaka, Cologne, Germany). The average of three consecutive readings was recorded. The vaginal pH was measured by collecting vaginal fluid using a sterile Puritan swab rotated once mid-vaginally and subsequently applied to color-coded pH paper (Macherey-Nagel, pH 4.0–7.0), as described previously (Donders et al., 2007).

 Bioinformatic processing Metagenomics pre-processing

Raw data were pre-processed with an in-house workflow.1 In short, the workflow removes the host genome reads and subsequently performs quality trimming of the reads. First, the host genome was removed using bowtie2 (version 2.4.2) by mapping reads to the human reference genome (Langmead and Salzberg, 2012). The parameters passed in bowtie2 included “--very-sensitive-local” and reference genome “GRCh38.p7.”2 Unmapped, paired reads were filtered from the output using samtools (version 1.11) (Li et al., 2009). Subsequent filtered reads were processed with fastp (version 0.20.1) performing quality trimming, adapter removal and low-complexity filtering (Chen et al., 2018). Fastp parameters included trim “--cut-right, −-cut_window_size 4 --cut_mean_quality 20”; minlen “-l 50”; adapter “--detected_adapter_for_pe”; complexity “-y.” The mean total reads before processing, after filtering for human reads and after quality trimming were calculated per study group and sample type (Supplementary Table S2).

 Microbial community profiling

Pre-processed reads were analyzed using MetaPhlAn (version 3.0.14) to profile the composition of the microbial communities and to predict read counts (Beghini et al., 2021). To profile the composition of the microbial communities, “--add_viruses” and “--unknown_estimation” were added besides default parameters. The outputs were merged with “merge_metaphlan_tables.py.” For the predicted read counts, parameters “-t rel_ab_w_read_stats,” “--add_viruses” and “--unknown_estimation” were included besides default parameters. Outputs were merged with an adapted version of the merge utility script. The pre-processed sequences mapped to 855 taxa (661 Bacteria, 152 Viruses, 33 Eukaryota and 9 Archaea).

 Data analysis and availability

The resulting abundance tables were analyzed and visualized using R version 4.1.2 (Vienna, Austria) (Team R, 2014). For data analyses, 89 clinical samples were available (LS, n = 39; HSIL, n = 20; Healthy, n = 30). Shannon diversity and Chao1 richness were computed on rarefied data with the phyloseq package (version 1.38.0) at species level (McMurdie and Holmes, 2013). The stat_compare_means from the ggpubr package (version 0.4.0) was used to compute overall significant differences with a Kruskal-Wallis test and to compare means between swab sites with the Wilcoxon Rank Sum Test (Kassambara, 2023). The abundance table was transformed to compositional data with the microbiome package (version 1.16.0). Subsequently, the mean relative abundance of the 10 most abundant genera was visualized with ggplot2 (version 3.3.6) for the bacterial and viral kingdoms on non-lesional and lesional skin of each study group (Lahti and Shetty, 2017; Wickham, 2023). For Linear discriminant analysis Effect Size (LEfSe), an object was created with the phyloseq2lefse function from the phyloseqCompanion package (version 1.1.) (Stagaman, 2023). Subsequent LEfSe analyses were performed until species level with default parameters (except LDA score > 4.0) on the Huttenhower lab Galaxy server to assess differences in relative abundance between lesional skin of LS patients and healthy vulvar skin, as well as lesional skin of vulvar HSIL patients and healthy vulvar skin (Segata et al., 2011). Aitchison distance was calculated for the Principal Coordinate Analyses (PCoA). On genus level data of the bacterial and viral kingdoms, Centered Log Ratio (CLR) transformation was performed using the transform function of the microbiome package. The distance function of the phyloseq package was used to generate a distance matrix with Euclidean distance. The betadisper function from the vegan package (version 2.6-4) was used to assess differences in variation between swab sites, while adonis2 was used to assess differences in centroids of the swab sites and study groups with constrained permutations for each patient.3 Alluvial plots were generated by calculating the mean relative abundance of the 10 most abundant bacterial and viral genera in each swab site for LS and vulvar HSIL. Since each swab site amounted to a total relative abundance of 100%, mean relative abundances for each swab site were normalized by dividing the relative abundances by the number of swab sites included. Alluvial plots were visualized using ggplot2 and ggalluvial (version 0.12.3) (Brunson, 2020).

 Quality control

Positive and negative controls were included during DNA extraction, and additional positive controls were included during sequencing. Positive and negative DNA extraction controls were compared for mean total reads. Mean total reads of the positive controls (4,742,574) were significantly higher (Wilcoxon Rank Sum Test, p = 0.03) compared to the mean total reads of the negative controls (12,565). Although having significantly lower number of reads, the negative controls contained skin-derived bacterium Cutibacterium acnes.

Additionally, the taxonomic species composition of the positive controls was compared to each other and to the theoretical mock community composition. All species of the mock community could be identified in the positive controls, except Bacillus subtilis. Instead, Bacillus intestinalis was identified, which has previously been reported as expected misclassification of Bacillus subtilis using the MetaPhlAn taxonomy database (Wright et al., 2023). The bacterial species are present in equal ratio’s, indicating overrepresentation (>15% increase) of Lactobacillus fermentum and underrepresentation (>15% decrease) of Pseudomonas aeruginosa, Enterococcus faecalis, Staphylococcus aureus, Listeria monocytogenes, Saccharomyces cerevisiae and Cryptococcus neoformans (Supplementary Figure S1).

 Data availability statement

All pre-processed metagenomics data have been deposited in the European Nucleotide Archive under accession number PRJEB61325.

 Results Cohort characteristics

In total, 25 women were included in the study. Baseline characteristics were comparable between groups (Table 1). Menopausal status and age were equally distributed. Most (24/25) participants were of Caucasian descent, with one healthy volunteer of mixed Caucasian and Latin American descent. All LS and vulvar HSIL patients had previously undergone one or multiple treatments for their vulvar condition, while healthy volunteers were naïve to any treatments of the vulvar skin.

Baseline characteristics.

Characteristics Healthy control(N = 10) Vulvar high grade squamous intraepithelial lesion(N = 5) Lichen sclerosus (N = 10)
Age in years—mean (range) 46.5 (25–73) 46.6 (32–66) 50.3 (25–72)
Pre-menopausal 5 3 5
Post-menopausal 5 2 5
Body mass index (BMI) in kg/m2—mean (range) 22.8 (19.6–27.6) 26.6 (21.6–30.0) 25.7 (18–30)
Ethnicity—N (%)
White 9 (90%) 5 (100%) 10 (100%)
Other 1 (10%)* 0 (0%) 0 (0%)
Smoking—N (%)
No 9 (90%) 1 (20%) 8 (80%)
Yes 1 (10%) 4 (80%) 2 (20%)
Disease duration in years—median (range) N/A 8 (7–25) 5.5 (1–12)
Vulvar squamous cell carcinoma in medical history 0 0 1 (10%)
Fitzpatrick skin type—N (%)
I 1 (10%) 1 (20%) 3 (30%)
II 4 (40%) 1 (20%) 4 (40%)
III 5 (50%) 3 (60%) 3 (30%)
HPV genotype biopsy—N (%)
HPV16 0 4 (80%) 0
HPV53 0 1 (20%)** 1 (10%)
Negative 10 (100%) 1 (20%) 9 (90%)
Previous vulvar treatments—N (%)
None 10 (100%) 0 0
Yes, 1 previous treatment 0 1 (20%) 7 (70%)
Yes, 2 previous treatments 0 1 (20%) 3 (30%)
Yes, 3 previous treatments 0 1 (20%) 0
Yes, 4 previous treatments 0 1 (20%) 0
Yes, 5 previous treatments 1 (20%)
Topical treatment 0 5*** 9****
Surgical 0 5 2*****
Coagulation 0 1 0
Laser 0 3 0
HPV vaccination (Gardasil) 0 2 0
Estriol/estradiol (vaginal) 0 0 2

*Mixed Latin American and Caucasian descent; **One HSIL patient was both HPV 16 and 53 positive; ***5x imiquimod, 1x 5-FU; ****2x triamcinolonacetonide, 8x dermovate; *****Surgical treatment encompassed vulvectomy and introitus plasty.

 Vulvar microbial skin diversity is decreased in LS compared to healthy skin

Diversity and richness of the bacterial and viral fraction were used to assess differences in microbial ecology between vulvar skin of healthy controls and of LS and HSIL patients. Healthy vulvar skin had a higher mean diversity compared to the lesional and non-lesional skin of LS patients (p = 0.029 and p = 0.001, respectively), while other comparisons to healthy skin were not significant (Figure 1). Interestingly, within the patient groups, lesional vulvar skin showed a non-significant rise in mean diversity and richness compared to non-lesional skin.

Diversity and richness of the vulvar skin microbiome. (A) Diversity measured by the Shannon index. (B) Richness measured by the Chao1 index. LS, lichen sclerosus; vHSIL, vulvar high-grade squamous intraepithelial lesions; NL, non-lesional (healthy appearing) vulvar skin; L, lesional vulvar skin. Only significant p-values are displayed.

 The bacterial fraction on diseased vulvar skin differs from healthy vulvar skin

To further assess differences in the microbial ecology of healthy vulvar skin and vulvar skin of LS and vulvar HSIL patients, taxonomic profiles of the bacterial and viral fraction were generated (Figure 2 and Supplementary Figures S2, S3). The bacterial and viral fraction are considered separately below.

Mean relative abundance on genus level of (A) the bacterial fraction and (B) the viral fraction of the microbiome in healthy controls, LS patients and HSIL patients. For each study group, the taxonomic profiles of non-lesional and lesional vulvar skin are displayed. No viral taxa were detected in some individuals (Supplementary Figure S1), therefore the bars do not accumulate to 100% in these summary graphs. LS, lichen sclerosus; vHSIL, vulvar high-grade squamous intraepithelial lesions; NL, non-lesional (healthy appearing) vulvar skin; L, lesional vulvar skin.

The main bacterial genera identified in non-lesional and lesional skin of LS and vulvar HSIL patients were similar as observed in healthy vulvar skin (Figure 2A). These bacterial genera included Prevotella, Lactobacillus, Staphylococcus, and Gardnerella.

LEfSe analyses were used to identify differentially abundant taxa between healthy vulvar skin and skin lesions in LS and vulvar HSIL. Regarding the bacterial fraction, lesional vulvar skin of LS patients was characterized by a depletion by taxa from the Prevotella genus and Bacteroidales order compared to healthy skin (Figures 3A,C,D). Meanwhile, lesional skin of HSIL patients showed an increase in the Fusobacteria phylum and depletion in the Actinobacteria class (Figures 4A,C,D).

Differentially abundant features in lesional vulvar skin of lichen sclerosus patients compared to healthy controls. (A) LEfSe analysis histogram of LDA scores (log10) showing overrepresented taxa (pink) and underrepresented taxa (green) in lesional vulvar skin of lichen sclerosus patients. Taxa with LDA scores above 4.0 were selected and unclassified taxa were excluded for viewing purposes. Bold taxa were differential abundant features classified to the lowest taxonomic level and used for (B–D). (B–D) Distribution of differential abundant features in individual samples per skin site with (B) showing overrepresented family Papillomaviridae, (C) the underrepresented genus Prevotella and (D) the underrepresented order Bacteroidales in lesional vulvar skin of lichen sclerosus patients. LS, lichen sclerosus; NL, non-lesional (healthy appearing) vulvar skin; L, lesional vulvar skin.

Differentially abundant features in lesional skin of vulvar HSIL patients compared to healthy controls. (A) LEfSe analysis histogram of LDA scores (log10) showing overrepresented taxa (red) and underrepresented taxa (green) in lesional vulvar skin of vulvar HSIL patients. Taxa with LDA scores above 4.0 were selected and unclassified taxa were excluded for viewing purposes. Bold taxa were differential abundant features classified to the lowest taxonomic level and used for (B–D). (B–D) Distribution of differential abundant features in individual samples per skin site with (B) showing the overrepresented genus Alphapapillomavirus, (C) showing the overrepresented phylum Fusobacteria and (D) showing the underrepresented class Actinobacteria in lesional vulvar skin of vulvar HSIL patients. vHSIL, vulvar high-grade squamous intraepithelial lesions; NL, non-lesional (healthy appearing) vulvar skin; L, lesional vulvar skin.

 <italic>Papillomaviridae</italic> are abundant both in LS and HSIL

For the viruses, various taxa were identified on the vulva of 6/10 healthy women, although predicted read counts were relatively low compared to the bacterial reads (565 ± 892 viral reads vs. 987,765 ± 1,140,330 bacterial reads) (Supplementary Figure S4). Viruses were more often identified in LS and vulvar HSIL patients (non-lesional LS: 9/9; lesional LS: 8/10; non-lesional HSIL: 4/5; lesional HSIL: 5/5; healthy controls: 6/10). The non-lesional and lesional skin in LS and HSIL patients mainly comprised viruses within the Papillomaviridae family, while these were present but not dominant in healthy vulvar skin (Figure 2B). Although viruses were prevalent, the abundance of the total viral fraction on diseased skin was not significantly increased compared to healthy vulvar skin (p = 0.15, Wilcoxon Rank Sum Test).

LEfSe results for the viral fraction showed that Papillomaviridae and Alphapapillomaviruses were significantly more abundant in vulvar lesional skin of LS and vulvar HSIL patients compared to healthy controls, respectively (Figures 3A,B, 4A,B). More specifically, various mucosal and cutaneous HPV types were identified in LS and vulvar HSIL (Supplementary Table S3). While non-lesional and lesional vulvar skin of HSIL patients mainly contained Alphapapillomaviruses 7 and 9 species (corresponding to clinical hrHPV types 18/45 and 16/31/33, respectively), the HPV profile of LS was more diverse (Cubie, 2013). Alphapapillomaviruses 3, 6 and 13 were detected in both non-lesional and lesional vulvar LS skin. Clinically, these species correspond to several high- and low-risk HPV types, but not HPV type 16 or 18.

 <italic>Prevotella</italic> spp. are shared between vulvar lesional skin, the vagina and anus in lichen sclerosus and vulvar high-grade squamous intraepithelial lesions

To further investigate the variation in the microbiome between (non-lesional and lesional) vulvar skin, vagina and anus, a Principal Coordinates Analysis was performed. The non-lesional and lesional vulvar skin showed overlap with the anus, suggesting the microbiome in vulvar skin is similar to the anal environment (Figure 5). However, PERMANOVA indicated the variation between sampling sites was significantly different (p = 0.001). Yet, the non-homogenous dispersion among the swab sites may have affected these PERMANOVA results.

Principal Coordinates Analysis (PCoA) of Aitchison distances between bacterial and viral communities associated to the swab sites including the vagina, non-lesional (NL) and lesional (L) vulva and anus. Betadisper p = 0.000; PERMANOVA p = 0.001 and r2 = 0.196.

To further assess horizontal transfer of the vaginal and anal environment to (lesional) vulvar skin and vice versa, the distribution and flow of genera was visualized (Figure 6). In both LS and HSIL patients, the Lactobacillus and Gardnerella genera were predominantly detected in the vagina and in lower relative abundances on the non-lesional and lesional skin and anus. Additionally, genera including Campylobacter, Corynebacterium, Finegoldia and Gardnerella were shared between the anus and non-lesional and lesional vulvar skin in LS and HSIL. Moreover, Alphapapillomavirus spp. were identified on all skin sites in vulvar HSIL patients. Interestingly, in LS patients Prevotella spp. was scarcely present on non-lesional vulvar skin while being detected on lesional vulvar, anal and vaginal skin sites.

Alluvial plots showing distribution of genera over the different sampling sites in (A) lichen sclerosus and (B) vulvar high-grade squamous intraepithelial lesions (HSIL) patients.

 Eukaryota and Archaea are rarely detected on vulvar skin

Compared to bacterial and viral taxa, Eukaryota and Archaea were rarely identified and in low abundance (Supplementary Table S4). The most observed fungus was Mallassezia globosa, which was detected in low relative abundance (up to 2%) on vulvar sites of 3/10 LS patients and 1/10 healthy volunteers while being absent in the vaginal or anal milieu. Methanobrevibacter smithii was the most frequently found Archaea, albeit in low relative abundance (up to 1%), in anal samples of 7 participants across groups.

 The barrier function of lesional vulvar skin in vulvar HSIL is impaired

The barrier function of the vulvar skin was assessed with TEWL. A significantly higher TEWL flux was observed between lesional vulvar HSIL [mean (SD) 79.2 g/m2/h (±44.9)] and healthy controls [mean (SD) 42.2 g/m2/h (±27.8); p = 0.043] (Figure 7). No significant differences in TEWL flux were observed between lesional LS [mean (SD) 57.5 g/m2/h (±36.8)] and HV, p = 0.309, nor between the non-lesional and the lesional sites of HSIL, p = 0.810. The lesional site of LS patients had a significant higher flux compared to the non-lesional LS site, p = 0.006.

Bar graph showing the trans epidermal water loss at baseline level for all groups on de non-biopsy site. The mean and standard deviations are displayed for each group. Trans epidermal water loss in g/m2/h (y-axis) plotted against measurements clustered per patient group (x-axis). LS, lichen sclerosus; vHSIL, vulvar high-grade squamous intraepithelial lesions; NL, non-lesional (healthy appearing) vulvar skin; L, lesional vulvar skin.

The mean vulvar pH of lesional LS and HSIL was 5.72 (SD 0.45) and 6.21 (SD 0.8), respectively, with no significant changes over time (Supplementary Figure S5). Non-lesional vulvar skin and vaginal pH of patients and healthy volunteers did not display group differences or changes over time (data not shown). Menopausal status was the stratifying factor for observed differences in pH in the LS and healthy control group.

 Discussion

This exploratory study shows that the vulvar skin of patients with LS and HSIL—both non-lesional and lesional—differs from healthy vulvar skin, which was most prominent for the viral fraction. Notably, papillomaviruses were characteristic for LS, despite being considered an HPV-independent condition. In the bacterial fraction, Prevotella spp. were less abundant in LS than healthy vulvar skin, but shared between vulvar skin, vaginal and anal sites, indicating potential horizontal transfer between sites. This study is the first to compare the vulvar microbiome composition of LS and HSIL to healthy controls using metagenomic sequencing techniques. The distinct bacteria and viruses found in both LS and vulvar HSIL warrant further investigation and may aid in future identification of driving factors or biomarkers for vulvar (pre)malignant diseases.

The bacterial fraction of vulvar LS skin shared many taxa with healthy vulvar skin, mainly comprising Lactobacillus, Prevotella, and Gardnerella genera, which were also previously identified in LS and vulvodynia (Pagan et al., 2021). Specifically in vulvar LS, depletion of Lactobacillus and gain of Prevotella species has previously been described in juvenile LS (Chattopadhyay et al., 2021). The opposite was observed in this study, where taxa from the Prevotella genus were reduced in LS compared to healthy vulvar skin. Meanwhile within the LS group, Prevotella spp. were shared between anal, vaginal and lesional LS skin sites, but almost absent on non-lesional vulvar skin. Thus, the exact role of Prevotella in LS remains inconclusive (Chattopadhyay et al., 2021; Liu et al., 2022). Relative abundance of Staphylococcus spp. seemed higher in non-lesional vulvar skin, which likely reflects the more skin-like sampling location toward the groin. Based on differential abundance analysis, a significant overrepresentation of papillomaviruses and an underrepresentation of taxa from the Prevotella genus and the Bacteriodales order were characteristic for LS compared to healthy controls. Interestingly, LS is generally considered a HPV-independent precursor of VSCC, although concurrent HPV-infections have been previously reported (Van Der Avoort et al., 2006; Aidé et al., 2010; Hald and Blaakaer, 2018). Co-occurrence of LS and papillomaviruses may be coincidental, given the high prevalence (~10%) of both high- and low-risk genital HPV infections in the general female population (Hald and Blaakaer, 2018; Kombe Kombe et al., 2020). Alternatively, our results may point toward an etiological role of HPV in LS, in which a subclinical HPV-infection could hypothetically induce LS in genetically and/or immunologically predisposed women. Such infection could occur through disruption of the skin barrier and local immune environment, known as the Koebner phenomenon (Tasker and Wojnarowska, 2013). This could apply for both high- and low-risk HPV types as observed in LS patients of our cohort. However, no significant increase of TEWL—i.e., decrease in skin barrier function—in LS skin was observed compared to healthy controls. Another explanation for the observed co-occurrence of HPV in LS could be that skin damage or immunosuppression from corticosteroid treatment in pre-existing LS facilitates colonization with papillomaviruses. These observations in LS lead to the hypothesis that women who eventually develop VSCC may display immunological variations resulting in less effective viral clearance.

Like LS, the vulvar skin of HSIL shared bacterial genera with healthy vulvar skin. The lesional skin of vulvar HSIL was characterized by significant increases in relative abundance of the Fusobacteria phylum and Alphapapillomavirus genus, while the Actinobacteria class was significantly reduced compared to healthy controls. The presence of high-risk Alphapapillomaviruses corroborates with the HPV-driven etiology of vulvar HSIL. Although no data on vulvar skin sites in HSIL is available hitherto, cervical dysplasia closely relates to vulvar HSIL and has been studied extensively. In cervical dysplasia, Sneathia, Mycoplasma and Prevotella species have been associated with co-infection and persistence of hrHPV types (Lee et al., 2013; Mitra et al., 2015; Audirac-Chalifour et al., 2016; Adebamowo et al., 2017; Klein et al., 2018; Łaniewski et al., 2019; Ritu et al., 2019; So et al., 2020). Prevotella has been designated as marker genus for cervical cancer, where it may influence HPV persistence through NOD-like receptor signaling (Wu et al., 2021). These findings indicate that functional processes, driven by the microbiome, potentially contribute to persistence or progression of cervical HPV-driven diseases. No significant overrepresentation of these taxa was identified in this study to vulvar HSIL. Furthermore, microbiome composition plays an integral role in the skin barrier, interacting with its other—physical, immunological and chemical—components (Harris-Tryon and Grice, 2022). For example, TEWL can change upon topical application of Lactobacillus and Corynebacterium formulations leading to increased and reduced TEWL, respectively (2013; Park et al., 2014; Kim et al., 2021; Li et al., 2022). Our results show that TEWL was significantly increased in affected vulvar HSIL skin, indicating a disrupted skin barrier function (Akdeniz et al., 2018; Alexander et al., 2018). Whether these observed differences represent a disease-driven disturbance of skin barrier function or are due to variability in sampling location will require confirmation in an expanded population.

Besides bacteria and viruses, metagenomic sequencing allowed for the identification of Archaea and Eukaryota. Solely two Archaea were identified in low prevalence and abundance. To date, no Archaeal taxa are associated with pathogenesis and are generally considered commensals (Aminov, 2013; Wirth and Young, 2020). Eukaryotic (i.e., fungal) Mallasezia globosa was also identified, albeit in low relative abundance, on non-lesional vulvar sites of LS patients, but not in healthy controls or vulvar HSIL. Previously, Mallasezia globosa was the most identified fungal species on healthy labia majora (Bruning et al., 2020). Candida taxa were detected once in this study, despite the reported C. albicans colonization rate of 20% in the general female population (Bradford and Ravel, 2017).

The main strengths of this study are the case–control trial design including patients and healthy controls with inter-participant lesional and non-lesional control. Vaginal and anal samples allowed for intra-individual comparison and correlation of results with literature. Also, this study is the first to investigate microbiome composition in vulvar HSIL. Limitations of this exploratory study mostly pertain to the low sample size and sequencing depth. As such, future studies are needed to confirm the findings described herein. In addition, our study did not include longitudinal analyses, although we assume that the vulvar and vaginal microbiome can be particularly subject to temporal changes. Factors dictating microbiome composition fluctuations may include cycle-related changes, demographic background and lifestyle choices including sexual activity, hair removal practices and intimate hygiene routines (Ravel et al., 2011; Shiraishi et al., 2011; Gajer et al., 2012; Hickey et al., 2015; Gupta et al., 2019; Song et al., 2020; Krog et al., 2022). Lastly, the absence of certain taxa or inability to identify a proportion of the sequences may be attributed to the low biomass samples as well as the non-amplification sequencing method.

The vulvar microbiome is a growing research field, with ongoing trials in LS (NCT05671263, NCT05147129), vulvar Paget’s disease (NCT03564483) and lichen planus (NCT05330572) (ClinicalTrials.gov, 2023a,b,c,d). Future studies should strive to include a variety of vulvar diseases, such as dVIN and VSCC (HPV-positive or HPV-negative) to capture the full disease spectrum from healthy vulvar skin to VSCC. Linking microbiome findings to changes in the tumor microenvironment may be further explored, as recently reported for VSCC (Rustetska et al., 2023). Novel treatments could be developed based on microbial targets, as previously attempted for bacterial vaginosis and genitourinary symptoms (Falagas et al., 2007; Armstrong et al., 2022; Yoshikata et al., 2022). Studies to new treatment modalities, such as the current clinical trial applying the topical JAK-inhibitor ruxolitinib in LS, could consider assessing the microbiome composition as exploratory biomarker (Rissmann et al., 2020; Papp et al., 2022; ClinicalTrials.gov, 2023e). The vulvar research field to date has mainly identified presence of taxa without appraisal their involvement in biophysical or pathologic processes, which should be the focus for future studies to unravel.

 Data availability statement

The data presented in the study are deposited in the European Nucleotide Archive repository, accession number PRJEB61325.

 Ethics statement

The studies involving humans were approved by the Medisch-Ethische Toetsingscommissie Leiden Den Haag Delft. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.

 Author contributions

LP: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Visualization, Writing – original draft. BH: Conceptualization, Data curation, Investigation, Methodology, Project administration, Writing – review & editing. MW: Data curation, Formal analysis, Software, Visualization, Writing – review & editing, Validation. RN: Formal analysis, Investigation, Methodology, Visualization, Writing – review & editing, Data curation. FS: Conceptualization, Investigation, Methodology, Supervision, Validation, Writing – review & editing. IS: Investigation, Methodology, Project administration, Writing – review & editing. WS: Investigation, Methodology, Supervision, Writing – review & editing. RZ: Conceptualization, Investigation, Methodology, Supervision, Writing – review & editing. JB: Conceptualization, Resources, Supervision, Writing – review & editing. RR: Conceptualization, Funding acquisition, Investigation, Methodology, Resources, Supervision, Writing – review & editing. JP: Conceptualization, Investigation, Resources, Supervision, Writing – review & editing. JH: Data curation, Formal analysis, Investigation, Methodology, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. MP: Conceptualization, Funding acquisition, Investigation, Methodology, Resources, Supervision, Writing – review & editing.

 Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This project was funded by the Research and Development fund of the Centre for Human Drug Research, in addition to financial contributions from the Bontiusstichting (project 8222-32146) of the Leiden University Medical Centre and individual donors of the Catharinaziekenhuis Onderzoeksfonds.

The authors would like to thank all patients and healthy volunteers for their participation. The patients in this manuscript have given written informed consent to publication of their case details. The authors thank Ellen Swanborn from the “Stichting Lichen Sclerosus” for aid in patient recruitment. The authors thank the donors of the Bontiusstichting of the Leiden University Medical Centre for their financial contribution to clinical research and the individual donors to the Catharina ziekenhuis Onderzoeksfonds for their personal efforts to support this project.

 Conflict of interest

JP has received grants for projects outside the scope of this article from Philips, Innosign, KWF, Catharina Onderzoeksfonds and Ruby and Rose.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2023.1264768/full#supplementary-material

1https://git.lumc.nl/snooij/metagenomics-preprocessing

2https://www.ncbi.nlm.nih.gov/assembly/GCF_000001405.33/

3https://cran.r-project.org/web/packages/vegan/index.html

References Adebamowo S. N. Ma B. Zella D. Famooto A. Ravel J. Adebamowo C. (2017). Mycoplasma hominis and Mycoplasma genitalium in the vaginal microbiota and persistent high-risk human papillomavirus infection. Front. Public Heal. 5, 1–10. doi: 10.3389/fpubh.2017.00140, PMID: 28695118 Aidé S. Lattario F. R. Almeida G. Do Val I. C. Da Costa C. M. (2010). Epstein-barr virus and human papillomavirus infection in vulvar lichen sclerosus. J. Low. Genit. Tract Dis. 14, 319–322. doi: 10.1097/LGT.0b013e3181d734f1, PMID: 20885159 Akdeniz M. Gabriel S. Lichterfeld-Kottner A. Blume-Peytavi U. Kottner J. (2018). Transepidermal water loss in healthy adults: a systematic review and meta-analysis update. Br. J. Dermatol. 179, 1049–1055. doi: 10.1111/BJD.17025, PMID: 30022486 Alexander H. Brown S. Danby S. Flohr C. (2018). Research techniques made simple: Transepidermal water loss measurement as a research tool. J. Invest. Dermatol. 138, 2295–2300.e1. doi: 10.1016/J.JID.2018.09.001, PMID: 30348333 Aminov R. I. (2013). Role of archaea in human disease. Front. Cell. Infect. Microbiol. 3:42. doi: 10.3389/fcimb.2013.00042, PMID: 23964350 Armstrong E. Kaul R. Cohen C. R. (2022). Optimizing the vaginal microbiome as a potential strategy to reduce heterosexual HIV transmission. J. Intern. Med. 293, 433–444. doi: 10.1111/joim.13600, PMID: 36544257 Audirac-Chalifour A. Torres-Poveda K. Bahena-Román M. Téllez-Sosa J. Martínez-Barnetche J. Cortina-Ceballos B. . (2016). Cervical microbiome and cytokine profile at various stages of cervical cancer: a pilot study. PLoS One 11:e0153274. doi: 10.1371/journal.pone.0153274, PMID: 27115350 Beghini F. McIver L. J. Blanco-Míguez A. Dubois L. Asnicar F. Maharjan S. . (2021). Integrating taxonomic, functional, and strain-level profiling of diverse microbial communities with bioBakery 3. Elife 10:10. doi: 10.7554/eLife.65088, PMID: 33944776 Bradford L. L. Ravel J. (2017). The vaginal mycobiome: a contemporary perspective on fungi in women’s health and diseases. Virulence 8, 342–351. doi: 10.1080/21505594.2016.1237332 Bruning E. Chen Y. McCue K. A. Rubino J. R. Wilkinson J. E. Brown A. D. G. (2020). A 28 day clinical assessment of a lactic acid-containing antimicrobial intimate gel wash formulation on skin tolerance and impact on the vulvar microbiome. Antibiot 9:55. doi: 10.3390/antibiotics9020055, PMID: 32024047 Brunson J. C. (2020). {ggalluvial}: layered grammar for alluvial plots. J. Open Source Softw. 5:2017. doi: 10.21105/joss.02017, PMID: 36919162 Byrd A. L. Belkaid Y. Segre J. A. (2018). The human skin microbiome. Nat. Rev. Microbiol. 16, 143–155. doi: 10.1038/nrmicro.2017.157 Chattopadhyay S. Arnold J. D. Malayil L. Hittle L. Mongodin E. F. Marathe K. S. . (2021). Potential role of the skin and gut microbiota in premenarchal vulvar lichen sclerosus: a pilot case-control study. PLoS One 16:e0245243. doi: 10.1371/journal.pone.0245243, PMID: 33444404 Chen S. Zhou Y. Chen Y. Gu J. (2018). Fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34, i884–i890. doi: 10.1093/bioinformatics/bty560, PMID: 30423086 ClinicalTrials.gov (2023a) Characterizing the microbiome in postmenopausal women with vulvar lichen sclerosus. Available at: https://clinicaltrials.gov/ct2/show/NCT05147129?cond=lichen+sclerosus&draw=2&rank=10 (Accessed January 11, 2023). ClinicalTrials.gov (2023b). Genital lichen sclerosus - epidemiology, comorbidities and the role of vulvar and penile microbiome. Available at: https://clinicaltrials.gov/ct2/show/NCT05671263?cond=lichen+sclerosus&draw=2&rank=2 (Accessed January 14, 2023. ClinicalTrials.gov (2023c). Gynecologic extramammary paget’s disease. Available at: https://clinicaltrials.gov/ct2/show/NCT03564483?term=microbiome&cond=vulva&draw=2&rank=3 (Accessed January 18, 2023). ClinicalTrials.gov (2023d). Clinical and molecular correlates of response to first line treatment in lichen planus. Available at: https://clinicaltrials.gov/ct2/show/NCT05330572?term=microbiome&cond=vulva&draw=2&rank=2 (Accessed January 18, 2023). ClinicalTrials.gov (2023e). A study to evaluate the efficacy and safety of ruxolitinib cream in participants with lichen sclerosus. Available at: https://clinicaltrials.gov/ct2/show/NCT05593445?cond=Lichen+Sclerosus&draw=2&rank=4 (Accessed January 3, 2023). Cubie H. A. (2013). Diseases associated with human papillomavirus infection. Virology 445, 21–34. doi: 10.1016/j.virol.2013.06.007 de Sanjosé S. Alemany L. Ordi J. Tous S. Alejo M. Bigby S. M. . (2013). Worldwide human papillomavirus genotype attribution in over 2000 cases of intraepithelial and invasive lesions of the vulva. Eur. J. Cancer 49, 3450–3461. doi: 10.1016/j.ejca.2013.06.033, PMID: 23886586 Donders G. G. G. G. Caeyers T. Tydhof P. Riphagen I. van den Bosch T. Bellen G. (2007). Comparison of two types of dipsticks to measure vaginal pH in clinical practice. Eur. J. Obstet. Gynecol. Reprod. Biol. 134, 220–224. doi: 10.1016/j.ejogrb.2006.07.016, PMID: 16952417 Falagas M. E. Betsi G. I. Athanasiou S. (2007). Probiotics for the treatment of women with bacterial vaginosis. Clin. Microbiol. Infect. 13, 657–664. doi: 10.1111/j.1469-0691.2007.01688.x Feng Q. Liang S. Jia H. Stadlmayr A. Tang L. Lan Z. . (2015). Gut microbiome development along the colorectal adenoma-carcinoma sequence. Nat. Commun. 6:6528. doi: 10.1038/ncomms7528 Gajer P. Brotman R. M. Bai G. Sakamoto J. Schütte U. M. Zhong X. . (2012). Temporal dynamics of the human vaginal microbiota. Sci. Transl. Med. 4:132ra52. doi: 10.1126/scitranslmed.3003605, PMID: 22553250 Grimm D. Eulenburg C. Brummer O. Schliedermann A. K. Trillsch F. Prieske K. . (2016). Sexual activity and function after surgical treatment in patients with (pre)invasive vulvar lesions. Support Care Cancer 24, 419–428. doi: 10.1007/s00520-015-2812-8, PMID: 26094599 Gupta S. Kakkar V. Bhushan I. (2019). Crosstalk between vaginal microbiome and female health: a review. Microb. Pathog. 136:103696. doi: 10.1016/j.micpath.2019.103696, PMID: 31449855 Hacker N. F. Eifel P. J. van der Velden J. (2012). Cancer of the vulva. Int. J. Gynecol. Obstet. 119, S90–S96. doi: 10.1016/S0020-7292(12)60021-6 Hald A. K. Blaakaer J. (2018). The possible role of human papillomavirus infection in the development of lichen sclerosus. Int. J. Dermatol. 57, 139–146. doi: 10.1111/ijd.13697, PMID: 28737238 Harris-Tryon T. A. Grice E. A. (2022). Microbiota and maintenance of skin barrier function. Science 376, 940–945. doi: 10.1126/SCIENCE.ABO0693 Hickey R. J. Zhou X. Settles M. L. Erb J. Malone K. Hansmann M. A. . (2015). Vaginal microbiota of adolescent girls prior to the onset of menarche resemble those of reproductive-age women. MBio 6:e00097-15. doi: 10.1128/mBio.00097-15, PMID: 25805726 Huisman B. W. Pagan L. Naafs R. G. C. ten Voorde W. Rissmann R. Piek J. M. J. . (2023). Dermatoscopy and optical coherence tomography in vulvar high-grade squamous intraepithelial lesions and lichen sclerosus: a prospective observational trial. J. Low. Genit. Tract Dis. 27, 255–261. doi: 10.1097/LGT.0000000000000731, PMID: 36924426 Kassambara A. Package ‘ggpubr.’ (2023). Available at: https://cran.r-project.org/web/packages/ggpubr/index.html Kim J.-H. Son S.-M. Park H. Kim B. K. Choi I. S. Kim H. . (2021). Taxonomic profiling of skin microbiome and correlation with clinical skin parameters in healthy Koreans. Sci. Rep. 11:16269. doi: 10.1038/s41598-021-95734-9, PMID: 34381083 Klein C. Gonzalez D. Samwel K. Kahesa C. Mwaiselage J. Aluthge N. . (2018). The relationship between the cervical microbiome, HIV status, and pre-cancerous lesions. bioRxiv 10, 1–12. doi: 10.1101/375592 Kombe Kombe A. J. Li B. Zahid A. Mengist H. M. Bounda G. A. Zhou Y. . (2020). Epidemiology and burden of human papillomavirus and related diseases, molecular pathogenesis, and vaccine evaluation. Front. Public Heal. 8:552028. doi: 10.3389/fpubh.2020.552028, PMID: 33553082 Krog M. C. Hugerth L. W. Fransson E. Bashir Z. Nyboe Andersen A. Edfeldt G. . (2022). The healthy female microbiome across body sites: effect of hormonal contraceptives and the menstrual cycle. Hum. Reprod. 37, 1525–1543. doi: 10.1093/humrep/deac094, PMID: 35553675 Lahti L. Shetty S. Microbiome R package. (2017). Available at: http://microbiome.github.io Langmead B. Salzberg S. L. (2012). Fast gapped-read alignment with bowtie 2. Nat. Methods 9, 357–359. doi: 10.1038/nmeth.1923, PMID: 22388286 Łaniewski P. Cui H. Roe D. J. Barnes D. Goulder A. Monk B. J. . (2019). Features of the cervicovaginal microenvironment drive cancer biomarker signatures in patients across cervical carcinogenesis. Sci. Rep. 9, 1–14. doi: 10.1038/s41598-019-43849-5 Lee J. E. Lee S. Lee H. Song Y. M. Lee K. Han M. J. . (2013). Association of the Vaginal Microbiota with human papillomavirus infection in a Korean twin cohort. PLoS One 8:e63514. doi: 10.1371/journal.pone.0063514, PMID: 23717441 Li H. Handsaker B. Wysoker A. Fennell T. Ruan J. Homer N. . (2009). The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079. doi: 10.1093/bioinformatics/btp352, PMID: 19505943 Li Z. Koban K. C. Schenck T. L. Giunta R. E. Li Q. Sun Y. (2022). Artificial intelligence in dermatology image analysis: current developments and future trends. J. Clin. Med. 11:6826. doi: 10.3390/jcm11226826, PMID: 36431301 Liu X. Zhuo Y. Zhou Y. Hu J. Wen H. Xiao C. (2022). Analysis of the vulvar skin microbiota in asymptomatic women and patients with vulvar lichen sclerosus based on 16S rRNA sequencing. Front. Cell Dev. Biol. 10:842031. doi: 10.3389/fcell.2022.842031, PMID: 35445011 Lockhart J. Gray N. M. Cruickshank M. E. (2013). The development and evaluation of a questionnaire to assess the impact of vulval intraepithelial neoplasia: a questionnaire study. BJOG 120, 1133–1143. doi: 10.1111/1471-0528.12229 McMurdie P. J. Holmes S. (2013). Phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One 8:e61217. doi: 10.1371/journal.pone.0061217 (2013). Management of symptomatic vulvovaginal atrophy: 2013 position statement of the North American Menopause Society. Menopause 20, 884–888. doi: 10.1097/GME.0b013e3182a122c2 Mitra A. MacIntyre D. A. Lee Y. S. Smith A. Marchesi J. R. Lehne B. . (2015). Cervical intraepithelial neoplasia disease progression is associated with increased vaginal microbiome diversity. Sci. Rep. 5, 1–11. doi: 10.1038/srep16865, PMID: 26574055 Mullen M. M. Merfeld E. C. Palisoul M. L. Massad L. S. Woolfolk C. Powell M. A. . (2019). Wound complication rates after vulvar excisions for premalignant lesions. Obstet. Gynecol. 133, 658–665. doi: 10.1097/AOG.0000000000003185, PMID: 30870300 Niemeyer-van der Kolk T. Buters T. P. Krouwels L. Boltjes J. de Kam M. L. van der Wall H. . (2022). Topical antimicrobial peptide omiganan recovers cutaneous dysbiosis but does not improve clinical symptoms in patients with mild to moderate atopic dermatitis in a phase 2 randomized controlled trial. J. Am. Acad. Dermatol. 86, 854–862. doi: 10.1016/j.jaad.2020.08.132, PMID: 33010325 Nitecki R. Feltmate C. M. (2018). Human papillomavirus and nonhuman papillomavirus pathways to vulvar squamous cell carcinoma: a review. Curr. Probl. Cancer 42, 476–485. doi: 10.1016/j.currproblcancer.2018.06.008, PMID: 30093165 Norenhag J. Du J. Olovsson M. Verstraelen H. Engstrand L. Brusselaers N. (2020). The vaginal microbiota, human papillomavirus and cervical dysplasia: a systematic review and network meta-analysis. BJOG 127, 171–180. doi: 10.1111/1471-0528.15854, PMID: 31237400 Pagan L. Ederveen R. A. M. Huisman B. W. Schoones J. W. Zwittink R. D. Schuren F. H. J. . (2021). The human vulvar microbiome: a systematic review. Microorganisms 9:2568. doi: 10.3390/microorganisms9122568, PMID: 34946169 Papp K. Szepietowski J. C. Kircik L. Toth D. Eichenfield L. F. Forman S. B. . (2022). Long-term safety and disease control with Ruxolitinib cream in atopic dermatitis: results from two phase 3 studies. J. Am. Acad. Dermatol. 88, 1008–1016. doi: 10.1016/j.jaad.2022.09.060, PMID: 36574595 Park S. B. Im M. Lee Y. Lee J. H. Lim J. Park Y. H. . (2014). Effect of emollients containing vegetable-derived lactobacillus in the treatment of atopic dermatitis symptoms: split-body clinical trial. Ann. Dermatol. 26, 150–155. doi: 10.5021/ad.2014.26.2.150, PMID: 24882967 Pérez-López F. R. Vieira-Baptista P. (2017). Lichen sclerosus in women: a review. Climacteric 20, 339–347. doi: 10.1080/13697137.2017.1343295, PMID: 28661705 Picardo S. L. Coburn B. Hansen A. R. (2019). The microbiome and cancer for clinicians. Crit. Rev. Oncol. Hematol. 141, 1–12. doi: 10.1016/j.critrevonc.2019.06.004 Piyathilake C. J. Ollberding N. J. Kumar R. Macaluso M. Alvarez R. D. Morrow C. D. (2016). Cervical microbiota associated with higher grade cervical intraepithelial neoplasia in women infected with high-risk human papillomaviruses. Cancer Prev. Res. 9, 357–366. doi: 10.1158/1940-6207.CAPR-15-0350, PMID: 26935422 Ravel J. Gajer P. Abdo Z. Schneider G. M. Koenig S. S. K. McCulle S. L. . (2011). Vaginal microbiome of reproductive-age women. Proc. Natl. Acad. Sci. U. S. A. 108, 4680–4687. doi: 10.1073/pnas.1002611107 Rissmann R. Moerland M. van Doorn M. B. A. (2020). Blueprint for mechanistic, data-rich early phase clinical pharmacology studies in dermatology. Br. J. Clin. Pharmacol. 86, 1011–1014. doi: 10.1111/bcp.14293, PMID: 32253783 Ritu W. Enqi W. Zheng S. Wang J. Ling Y. Wang Y. (2019). Evaluation of the associations between cervical microbiota and HPV infection, clearance, and persistence in cytologically normal women. Cancer Prev. Res. 12, 43–56. doi: 10.1158/1940-6207.CAPR-18-0233, PMID: 30463989 Rustetska N. Szczepaniak M. Goryca K. Zalewska E. B. Figat M. (2023). The intratumour microbiota and neutrophilic inflammation in squamous cell vulvar carcinoma microenvironment. J. Transl. Med. 21:285. doi: 10.1186/s12967-023-04113-7, PMID: 37118737 Segata N. Izard J. Waldron L. Gevers D. Miropolsky L. Garrett W. S. . (2011). Metagenomic biomarker discovery and explanation. Genome Biol. 12:R60. doi: 10.1186/gb-2011-12-6-r60, PMID: 21702898 Shiraishi T. Fukuda K. Morotomi N. Imamura Y. Mishima J. Imai S. . (2011). Influence of menstruation on the microbiota of healthy women’s labia minora as analyzed using a 16s rRNA gene-based clone library method. Jpn. J. Infect. Dis. 64, 76–80. doi: 10.7883/yoken.64.76, PMID: 21266747 So K. A. Yang E. J. Kim N. R. Hong S. R. Lee J. H. Hwang C. S. . (2020). Changes of vaginal microbiota during cervical carcinogenesis in women with human papillomavirus infection. PLoS One 15:e0238705. doi: 10.1371/journal.pone.0238705, PMID: 32941440 Song S. D. Acharya K. D. Zhu J. E. Deveney C. M. Walther-Antonio M. R. S. Tetel M. J. . (2020). Daily vaginal microbiota fluctuations associated with natural hormonal cycle, contraceptives, diet, and exercise. mSphere 5:e00593-20. doi: 10.1128/msphere.00593-20, PMID: 32641429 Stagaman K. phyloseqCompanion: provides additional functions to work with phyloseq objects. R package version 1.1. (2023). https://rdrr.io/github/kstagaman/phyloseqCompanion/ Tasker G. L. Wojnarowska F. (2013). Lichen sclerosus. J. Dermatol. Nurses Assoc. 5, 213–216. doi: 10.1097/JDN.0b013e31829cde87 Team R (2014). R: a language and environment for statistical computing MSOR Connections. Thuijs N. B. van Beurden M. Bruggink A. H. Steenbergen R. D. M. M. Berkhof J. Bleeker M. C. G. G. (2021). Vulvar intraepithelial neoplasia: incidence and long-term risk of vulvar squamous cell carcinoma. Int. J. Cancer 148, 90–98. doi: 10.1002/ijc.33198, PMID: 32638382 Van Der Avoort I. A. M. Shirango H. Hoevenaars B. M. Grefte J. M. M. de Hullu J. A. de Wilde P. C. M. . (2006). Vulvar squamous cell carcinoma is a multifactorial disease following two separate and independent pathways. Int. J. Gynecol. Pathol. 25, 22–29. doi: 10.1097/01.pgp.0000177646.38266.6a van de Nieuwenhof H. P. van Kempen L. C. L. T. L. T. De Hullu J. A. Bekkers R. L. M. Bulten J. Melchers W. J. G. . (2009). The etiologic role of HPV in vulvar squamous cell carcinoma fine tuned. Cancer Epidemiol. Biomarkers Prev. 18, 2061–2067. doi: 10.1158/1055-9965.EPI-09-0209, PMID: 19567503 Voss F. O. Thuijs N. B. Vermeulen R. F. M. Wilthagen E. A. van Beurden M. Bleeker M. C. G. (2021). The vulvar cancer risk in differentiated vulvar intraepithelial neoplasia: a systematic review. Cancers 13:6170. doi: 10.3390/CANCERS13246170/S1, PMID: 34944788 Wickham H. (2023) Ggplot2: elegant graphics for data analysis. Available at: https://books.google.com/books/about/ggplot2.html?hl=nl&id=XgFkDAAAQBAJ (Accessed January 16, 2023). Wirth J. Young M. (2020). The intriguing world of archaeal viruses. PLoS Pathog. 16:e1008574. doi: 10.1371/JOURNAL.PPAT.1008574, PMID: 32790746 Wright M. L. Podnar J. Longoria K. D. Nguyen T. C. Lim S. Garcia S. . (2023). Comparison of commercial DNA extraction kits for whole metagenome sequencing of human oral, vaginal, and rectal microbiome samples. bioRxiv. doi: 10.1101/2023.02.01.526597 Wu S. Ding X. Kong Y. Acharya S. Wu H. Huang C. . (2021). The feature of cervical microbiota associated with the progression of cervical cancer among reproductive females. Gynecol. Oncol. 163, 348–357. doi: 10.1016/J.YGYNO.2021.08.016, PMID: 34503848 Yoshikata R. Yamaguchi M. Mase Y. Tatsuyuki A. Myint K. Z. Y. Ohta H. (2022). Evaluation of the efficacy of Lactobacillus-containing feminine hygiene products on vaginal microbiome and genitourinary symptoms in pre- and postmenopausal women: a pilot randomized controlled trial. PLoS One 17:e0270242. doi: 10.1371/journal.pone.0270242, PMID: 36584204