Pathobionts in the Vaginal Microbiota: Individual Participant Data Meta-Analysis of Three Sequencing Studies

Sequencing studies have shown that optimal vaginal microbiota (VMB) are lactobacilli-dominated and that anaerobes associated with bacterial vaginosis (BV-anaerobes) are commonly present. However, they overlooked a less prevalent but more pathogenic group of vaginal bacteria: the pathobionts that cause maternal and neonatal infections and pelvic inflammatory disease. We conducted an individual participant data meta-analysis of three VMB sequencing studies that included diverse groups of women in Rwanda, South Africa, and the Netherlands (2,044 samples from 1,163 women in total). We identified 40 pathobiont taxa but only six were non-minority taxa (at least 1% relative abundance in at least one sample) in all studies: Streptococcus (54% of pathobionts reads), Staphylococcus, Enterococcus, Escherichia/Shigella, Haemophilus, and Campylobacter. When all pathobionts were combined into one bacterial group, the VMB of 17% of women contained a relative abundance of at least 1%. We found a significant negative correlation between relative abundances (ρ = −0.9234), but not estimated concentrations (r = 0.0031), of lactobacilli and BV-anaerobes; and a significant positive correlation between estimated concentrations of pathobionts and BV-anaerobes (r = 0.1938) but not between pathobionts and lactobacilli (r = 0.0436; although lactobacilli declined non-significantly with increasing pathobionts proportions). VMB sequencing data were also classified into mutually exclusive VMB types. The overall mean bacterial load of the ≥20% pathobionts VMB type (5.85 log₁₀ cells/μl) was similar to those of the three lactobacilli-dominated VMB types (means 5.13–5.83 log₁₀ cells/μl) but lower than those of the four anaerobic dysbiosis VMB types (means 6.11–6.87 log₁₀ cells/μl). These results suggest that pathobionts co-occur with both lactobacilli and BV-anaerobes and do not expand as much as BV-anaerobes do in a dysbiotic situation. Pathobionts detection/levels were increased in samples with a Nugent score of 4–6 in both studies that conducted Nugent-scoring. Having pathobionts was positively associated with young age, non-Dutch origin, hormonal contraceptive use, smoking, antibiotic use in the 14 days prior to sampling, HIV status, and the presence of sexually transmitted pathogens, in at least one but not all studies; inconsistently associated with sexual risk-taking and unusual vaginal discharge reporting; and not associated with vaginal yeasts detection by microscopy. We recommend that future VMB studies quantify common vaginal pathobiont genera.

Introduction

Understanding of the vaginal microbiota (VMB) has increased significantly since the turn of the century due to the increased availability of molecular laboratory techniques such as next-generation sequencing (van de Wijgert et al., 2014). Molecular studies have shown that most women have a VMB consisting of lactobacilli (most commonly Lactobacillus crispatus or L. iners), but that vaginal dysbiosis is highly prevalent worldwide (van de Wijgert and Jespers, 2017). The most common type of vaginal dysbiosis is anaerobic dysbiosis, which is characterized by a decrease of lactobacilli and an increase of fastidious anaerobes (van de Wijgert et al., 2014). Clinicians refer to symptomatic anaerobic dysbiosis as bacterial vaginosis (BV): patients typically have mild vaginal inflammation and a fishy-smelling vaginal discharge. It should be noted, however, that anaerobic dysbiosis is also frequently asymptomatic. The VMB of most women with anaerobic dysbiosis consists of a highly diverse mixture of fastidious anaerobes, usually including Gardnerella vaginalis. However, a substantial proportion of women with anaerobic dysbiosis are dominated by G. vaginalis, and this type of low diverse anaerobic dysbiosis is often overlooked. Recent studies have suggested that these women might be more difficult to treat, potentially due to the presence of a G. vaginalis-initiated vaginal mucosal biofilm (Verwijs et al., 2019a).

Another clinically relevant type of vaginal dysbiosis that has systematically been overlooked is the presence of bacterial pathobionts in the VMB (van de Wijgert and Jespers, 2017). Microbiologists define the term pathobiont as any potentially pathological organism which, under normal circumstances, lives as a non-harming symbiont. In the vaginal niche, this would include—among others—Streptococcus agalactiae (Group B streptococcus), Staphylococcus aureus, and species in the Enterobacteriaceae family. These bacteria have often been associated with maternal and neonatal infections (Cools et al., 2016; Black et al., 2018), as well as invasive infections in non-pregnant women such as pelvic inflammatory disease (Brunham et al., 2015). Some clinical researchers have hypothesized that a distinct type of vaginitis (desquamative inflammatory vaginitis), which is characterized by much more severe vaginal inflammation than BV and with desquamation of vaginal epithelial cells including parabasal cells (Sobel, 1994; Paavonen and Brunham, 2018), may be caused by pathobionts in the VMB (Donders et al., 2017). Two cases that appear to have been triggered by toxic shock syndrome toxin-1-producing Staphylococcus aureus strains have indeed been reported (Pereira et al., 2013). However, others believe that the condition is caused by estrogen deficiency or an immunologic disorder, and that vaginal dysbiosis develops secondarily (Sobel et al., 2011). A recent study found that most patients with vaginitis, parabasal cells, and lactobacilli-deficiency by microscopy did not have consistent VMB patterns by VMB sequencing (Oerlemans, 2019). We conclude that there is sufficient evidence to consider VMB pathobionts clinically relevant, but that the evidence-base related to both symptoms and complications is weak.

An important reason why the evidence-base is weak is because pathobionts are often not assessed properly. For example, neonatal invasive infection studies have focused on only one pathobiont (S. agalactiae) by culture (Kwatra et al., 2016), and VMB sequencing studies have systematically under-reported pathobionts. Authors of such studies typically use bioinformatical methods, such as hierarchical clustering, to summarize the sequencing data into a few VMB types. The first set of VMB types were published by Ravel et al. (2011) based on a study in asymptomatic American women: they referred to these as community state types I (L. crispatus-dominated), II (L. gasseri-dominated), III (L. iners-dominated), IV (diverse group), and V (L. jensenii-dominated (Ravel et al., 2011). The only pathobiont that was mentioned in this publication was Streptococcus, as one of the taxa included in the “diverse group.” However, hierarchical clustering only takes relative abundances into account and not the pathogenic potential of individual bacteria. Pathobionts usually occur at lower levels than BV-anaerobes, but have a higher pathogenic potential: these lower levels may therefore be clinically relevant. Because pathobionts rarely dominate the VMB, samples that contain pathobionts are often classified based on the other bacteria that are also present in that sample. For example, a sample containing 70% L. iners and 30% S. agalactiae would be classified as community state type III (L. iners-dominated) in most studies.

We believe that this vaginal pathobionts knowledge gap is hampering clinical progress in the field. We therefore performed an individual participant data meta-analysis of three VMB sequencing studies that enrolled diverse groups of women in Rwanda, South Africa, and the Netherlands (2,044 samples from 1,163 women in total), with as main aim to describe the presence and levels of all pathobionts identified in the sequencing data, their correlations with lactobacilli and BV-anaerobes, and their associations with participant sociodemographic, behavioral, and clinical/laboratory characteristics.

Materials and Methods

Studies Included in the Meta-Analysis

We performed an individual participant data meta-analysis of three VMB sequencing studies that were conducted in three different countries to account for regional and ethnic differences in VMB composition: (1) a clinical trial of intermittent oral metronidazole or vaginal probiotic use in Kigali, Rwanda (referred to as the Rwanda VMB study); (2) the VMB sub-study of the South African HPV in Africa Research Partnership (HARP) study in Johannesburg, South Africa; and (3) the VMB sub-study of the Healthy Life in an Urban Setting (HELIUS) study in Amsterdam, the Netherlands.

The Rwanda VMB study screened HIV-negative, non-pregnant, pre-menopausal women at high risk of sexually transmitted infections (STIs) for BV (van de Wijgert et al., 2020a). Women with BV were treated with metronidazole for seven days, and when cured of BV and other urogenital infections, were randomized to no intervention, or intermittent use of oral metronidazole or two different lactobacilli-containing vaginal probiotics for 2 months. The lactobacilli contained in the vaginal probiotics did not include any naturally occurring vaginal lactobacilli. Women were sampled at screening (start of BV/urogenital infection treatment, if applicable), enrollment (start of the interventions), Day 7, Month 1, Month 2 (cessation of the interventions), and Month 6. The study found that all three interventions were safe and affected naturally occurring lactobacilli and BV-anaerobes (in favor of the lactobacilli, particularly L. iners), but not pathobionts. However, to avoid any bias in this meta-analysis due to exposure to interventions, we conducted analyses that included lactobacilli and BV-anaerobes levels on VMB data that were not influenced by the interventions (N = 366 of 629 samples): data from samples collected in all randomization groups at the screening visit prior to any treatments (if applicable) and at the Month 6 visit (4 months after cessation of the interventions), as well as samples collected in the no intervention group at the Month 1 and Month 2 visits.

The VMB sub-study of the HARP study was a nested case-control study within a prospective cohort study conducted in Johannesburg, South Africa (van de Wijgert et al., 2020b). The study enrolled HIV-positive women and investigated the associations of VMB composition with high-risk human papillomavirus (hrHPV) and cervical intraepithelial neoplasia (CIN) acquisition, clearance, and/or persistence. All but one participant were of sub-Saharan African origin. Samples for VMB analyses were collected at baseline (N = 445) and at endline (N = 414), a median of 16 months later. The study concluded that hrHPV infection (and/or increased sexual risk-taking) likely causes anaerobic vaginal dysbiosis, but that a bidirectional relationship is also possible. Furthermore, in this population, dysbiosis did not increase CIN2+ risk, but CIN2+ increased dysbiosis risk. Since the study did not include an intervention, we used all available VMB data for the analyses presented in this paper.

The HELIUS study is a large, multi-ethnic cohort study in Amsterdam, the Netherlands (Snijder et al., 2017). Sampling was stratified by ethnic group and included the six largest ethnic groups in the city (Dutch, African Surinamese, South-Asian Surinamese, Turkish, Moroccan, and Ghanaian). In a sub-sample, a cross-sectional study on the association of ethnicity with VMB composition was performed (Borgdorff et al., 2017). For this sub-study, vaginal samples of 546 pre-menopausal women were sequenced. The most prevalent VMB composition in ethnically Dutch women was a L. crispatus-dominated VMB, in African Surinamese and Ghanaian women a polybacterial G. vaginalis-containing VMB, and in the other ethnic groups a L. iners-dominated VMB. This study did not include an intervention either, and we therefore used all available VMB data for the analyses presented in the current paper.

Sequencing and Other Laboratory Methods

All three studies extracted DNA from vaginal swabs and conducted 16S rRNA gene sequencing of the V3–V4 variable regions on Illumina platforms (San Diego, CA, USA) as described by Fadrosh et al. (2014). Standard diagnostic tests were used to test for STIs (all three studies), and BV and vulvovaginal candidiasis (Rwanda VMB and HARP studies only). There were some differences in sequencing and diagnostic methods used, as outlined in Table 1. Because of these differences, we conducted all analyses on each study separately as well as the three studies combined.

TABLE 1

Table 1. Main study characteristics of the three studies.

The BactQuant 16S rRNA gene qPCR (Liu et al., 2012) was only done in the Rwanda VMB study (N = 379, of which 158 samples were not influenced by the interventions). The 16S rRNA gene concentration per sample was used to convert the relative abundances of that sample into estimated concentrations as previously described (van de Wijgert et al., 2020a). We used both relative abundances as well as estimated concentrations for the Rwanda VMB study, but only had relative abundances for the other two studies.

Sequencing Data Processing

The 16S rRNA gene sequencing and initial data processing yielded two-dimensional tables with samples and bacterial taxa on the axes, and relative abundances in the cells, for each of the three studies. DADA2 assigns amplicon sequence variants (ASVs) to taxa (in the Rwanda VMB study) and Swarm and USEARCH assign them to operational taxonomic units (OTUs; in the HARP and HELIUS studies, respectively). Details on quality control and cleaning of reads, taxonomic assignments, conversion of read counts into relative abundances, and rarefaction are summarized in Table 1 and explained in the original publications. The three study-specific relative abundance tables were combined into a single table, and all subsequent data processing steps were redone for the combined table (i.e., are slightly different from the original publications) to ensure that they were identical for the three studies.

Data reduction for biostatistical modeling was done in three different ways. First, the Simpson diversity index (1-D) was calculated for each sample, ranging from 0 (no diversity) to 1 (infinite diversity). Second, each ASV/OTU was assigned to one of four “bacterial groups” based on the published literature (Supplementary Material 1) as follows: (1) lactobacilli; (2) BV-anaerobes; (3) pathobionts; and (4) a rest group called “other bacteria” (which contained mostly skin and Bifidobacteria). Pathobionts were defined as all bacterial taxa that have been reported in the literature as having been associated with invasive disease, and are not typically associated with BV; we also included STI pathogens in this category because their mean relative abundances were too low to justify a separate bacterial group. For each sample, relative abundances of ASVs/OTUs belonging to the same bacterial group were summed. This resulted in four relative abundances (one for each bacterial group) per sample, which sum to one in total. For example, one sample could contain 0.5 (50%) lactobacilli reads, 0.4 (40%) BV-anaerobes reads, 0.08 (8%) pathobionts reads, and 0.02 (2%) other bacteria reads. Third, we classified samples into nine VMB types (with each sample assigned to only one VMB type): (1) Lactobacillus iners-dominated (Li; ≥75% relative abundance of lactobacilli of which L. iners was the most common); (2) L. crispatus-dominated (Lcr; also ≥75% lactobacilli of which L. crispatus was the most common); (3) dominated by other Lactobacillus species (Lo; also containing ≥75% lactobacilli); (4) lactobacilli and anaerobes (LA; ≥25% lactobacilli with the remainder BV-anaerobes); (5) high diversity BV-anaerobes with ≥10% G. vaginalis presence (BV_GV); (6) high diversity BV-anaerobes with <10% G. vaginalis presence (BV_noGV); (7) G. vaginalis-dominated (GV; G. vaginalis ≥50%); (8) substantial presence of pathobionts (PB; ≥20% pathobiont taxa); and (9) Bifidobacterium-dominated (BD; ≥50% Bifidobacteria).

Statistical Analyses and Figures

Statistical analyses were performed using Stata version 13 (StataCorp, College Station, TX, USA) and R version 3.2.3 (R foundation for Statistical Computing 2016, Vienna, Austria). All analyses were cross-sectional, sometimes including samples collected at baseline only (one sample per woman) and sometimes including all samples (in case of the Rwanda VMB and HARP studies, more than one sample per woman). Women in the Rwanda VMB study were exposed to antibiotic and/or probiotic interventions, and samples that could potentially have been influenced by these interventions were excluded from most analyses as described above and as indicated in the tables and text. Unadjusted differences between groups of interest were tested by Fisher's exact test for binary variables, Chi-squared test for categorical variables, and Kruskal-Wallis test for continuous variables. Pathobiont levels (relative abundances or estimated concentrations) were correlated with those of other bacterial groups or taxa by Spearman's rank correlation when all samples were included and by Pearson's correlation coefficient when samples with <1% pathobionts were excluded. To assess sociodemographic, behavioral, and clinical determinants of pathobionts detection (≥1% compared to <1%) and levels, we used unadjusted logistic regression models for analyses including one sample per woman, and Kruskal-Wallis tests for analyses of pathobionts levels that included all samples. The heatmap was made with the gplots package in R (Warnes et al., 2016), bar charts in Stata, and correlation matrices with the corrplot package in R (Taiyun, 2019).

Results

Participant Characteristics

The median age in the three studies combined was 30 years (inter-quartile range 26–34) and most women were non-pregnant by design (Table 2). The majority of women in the Rwanda VMB study were at high risk of HIV/STI by design: 93.2% reported two or more partners in the month prior to the baseline visit and 76.6% reported no or inconsistent condom use. Much lower proportions of women in the HARP study reported these current sexual risk behaviors but they were all HIV-positive by design. Women in the HELIUS study were not selected based on sexual risk or HIV-status. The proportions of women reporting sexual risk behaviors could be considered average for a young, urban, Dutch population, but with differences by ethnic group: proportions were highest in women of Dutch origin, followed by Dutch women of sub-Saharan African (African Surinamese or Ghanaian) origin, and Dutch women of Turkish, Moroccan, or South Asian Surinamese origin. Current hormonal contraceptive use varied substantially between studies and ethnic groups, as did current smoking habits. Almost half of the women (39.4–45.1%; not assessed in the HARP study) reported current urogenital symptoms, but none of them had sought care for them. Laboratory-confirmed viral and bacterial STI prevalences were high in the Rwanda VMB study and low in the HELIUS study, whereas viral STI prevalences were high and bacterial STI prevalences low in the HARP study (consistent with high past but low current sexual risk). Antibiotic use in the 2 weeks prior to baseline was rare.

TABLE 2

Table 2. Baseline characteristics by study, and by ethnic group within the HELIUS study.

Overall Vaginal Microbiota Characteristics

A heatmap of key taxa for all 2,044 samples from all three studies combined is shown in Supplementary Material 2, Table S1. In Table 3, Figure 1, VMB study samples were stratified by exposure to interventions, and HELIUS study samples by ethnic group. Mean Simpson diversity indexes and mean relative abundances of bacterial groups, and key taxa within these groups, differed significantly between the three studies and these pre-specified strata within studies (Table 3, Figure 1A). By far the most common bacterial groups in all studies and strata were the lactobacilli and BV-anaerobes, with mean relative abundances ranging from 0.46 to 0.73 and 0.25 to 0.49, respectively. The differences between studies and strata were as expected, with lower lactobacilli and higher BV-anaerobes proportions in women with higher sexual risk profiles and/or STI exposures and in women of sub-Saharan African ethnicities. In contrast, Rwanda VMB study participants who had recently been exposed to antibiotic/probiotic interventions had higher lactobacilli and lower BV-anaerobes proportions. The differences in L. crispatus mean relative abundance were especially striking, ranging from only 0.03 in the Rwanda VMB study samples that were not influenced by interventions to 0.38 in the HELIUS samples from women of Dutch origin. Mean relative abundances for pathobionts and the “other bacteria” group were low in all studies and strata, ranging from 0.01 to 0.07 and 0 to 0.05, respectively. The mean pathobionts relative abundance did not show a clear pattern between studies and strata but was lowest in the HELIUS women of Dutch origin. Estimated concentrations were only available for the Rwanda VMB study, and mean estimated concentrations in log₁₀ cells/μl were 5.12 for lactobacilli (mostly consisting of L. iners), 5.17 for BV-anaerobes, 2.18 for pathobionts, and 1.92 for other bacteria in samples not influenced by interventions (Table 3). The mean pathobionts estimated concentration was therefore 871 times lower than the mean lactobacilli estimated concentration, and 977 times lower than the mean BV-anaerobes estimated concentration.

TABLE 3

Table 3. Overview of VMB composition characteristics by study, stratified by intervention exposure (VMB Rwanda) and ethnic group (HELIUS).

FIGURE 1

Figure 1. Bar charts by study, intervention exposure (VMB Rwanda) and ethnic group (HELIUS). (A) Mean relative abundance of bacterial groups. (B) Cumulative percentage of women with a specific VMB type. BV, bacterial vaginosis; BD, Bifidobacterium-dominated; BV_GV, polybacterial Gardnerella vaginalis-containing; BV_noGV, polybacterial but low G. vaginalis; CI, confidence interval; GV, G. vaginalis-dominated; interv, (study product) interventions; LA, lactobacilli and anaerobes; Lcr, L. crispatus-dominated; Li, L. iners-dominated; Lo, other lactobacilli-dominated; NA, not applicable; PB, pathobionts-containing; SA, South-Asian; SSA, sub-Saharan African; VMB, vaginal microbiota; VMB-RW, Rwanda VMB study. Rwanda VMB study samples collected at the screening and Month 6 visits in all randomization groups, and at the Month 1 and Month 2 visits in the no-intervention group, were considered not influenced by the interventions.

The VMB types for all samples combined (N = 2,044) were distributed as follows: Li 31.6%, Lcr 10.5%, Lo 2.4%, LA 15.8%, BV_GV 20.3%, BV_noGV 5.5%, GV 8.3%, PB 5.2%, and BD 0.5%. The latter VMB type included only 10 samples and was therefore not included in subsequent comparisons. Consistent with the bacterial group findings, the VMB types characterized by lactobacilli-domination (Li, Lcr, and Lo; 44.5%) or by anaerobic dysbiosis (LA, BV_GV, BV_noGV, GV; 49.9%) were much more common than the VMB type characterized by ≥20% pathobionts (5.2%). VMB type distributions differed significantly between the studies and strata, following the same patterns as described above for the bacterial group findings (Table 3, Figure 1B).

Identification of Common Vaginal Pathobionts

We identified 40 different pathobiont taxa in all 2,044 samples combined (Supplementary Material 1; reported at species level if only one species was identified, genus level if multiple species and/or the genus was identified, and family or class level if only that level was identified). However, 20 of these were never a non-minority taxon (defined as present at a relative abundance of at least 1% in at least one sample) in any of the studies. Only six taxa were a non-minority taxon in all three studies: Streptococcus, Staphylococcus, Enterococcus, Escherichia/Shigella, Haemophilus, and Campylobacter. Chlamydia was consistently detected in all three studies but only as a non-minority taxon in HELIUS, and Neisseria and Treponema were detected in the two African studies only. The remaining 11 taxons varied in their detection (yes vs. no) and relative abundance (minority vs. non-minority) status between the three studies. More than half (54%) of all pathobiont sequencing reads were assigned to Streptococcus genus/species and 24% of the Streptococcus reads were assigned to S. agalactiae or S. agalactiae/pyogenes. Seventeen percent of all women (196/1,153) at baseline, and 19.5% (399/2,044) of all samples, had at least 1% pathobionts in their VMB; these proportions were 12.7% (147/1,153), and 14.9% (304/2,044) for at least 1% Streptoccoccus. Among samples with ≥20% pathobionts (N = 107; Supplementary Material 2, Table S1), 33 contained Streptococcus genus/species as the only pathobionts (relative abundances of 0.53–0.98) and an additional 52 contained multiple pathobionts including substantial relative abundances of Streptococcus genus/species (0.13–0.73). The other 22 samples contained other pathobionts (most commonly staphylococci, Escherichia/Shigella species, Haemophilus species, and/or enterococci), with <5% Streptococcus.

Of note, the total estimated bacterial concentration differed significantly per VMB type (Table 4; data available for the Rwanda VMB study only). The mean total bacterial estimated concentration of women with the PB VMB type (5.85 log₁₀ cells/μl) was comparable to those of women with Lactobacillus-dominated VMB types (5.13–5.83 log₁₀ cells/μl) but lower than those of women with VMB types associated with anaerobic dysbiosis (6.11–6.87 log₁₀ cells/μl). When samples were stratified by the proportion of pathobionts in the VMB (<1%, 1- <10%, 10%- <20%, 20- <50%, and ≥50%), the mean estimated concentration of pathobionts increased as expected, but reached a plateau at proportions of 10% or more. The mean estimated concentration of total bacteria remained stable but declined somewhat when the pathobionts proportion reached above 20%. Results were similar when only samples not influenced by interventions were included in these analyses (Table 4).

TABLE 4

Table 4. Estimated bacterial concentration per VMB type and by proportion of pathobionts (Rwanda VMB study only, multiple samples per woman).

Correlations Between Vaginal Pathobionts and Other VMB Characteristics

We next investigated correlations between pathobionts levels (relative abundances or estimated concentrations), lactobacilli levels, and BV-anaerobes levels for samples not influenced by interventions (Table 5). With increasing pathobionts proportion (from <1% to 1- <10% to 10%- <20%, etc.), the mean relative abundance of lactobacilli declined significantly (ρ = −0.1851; 95% confidence interval (CI) −0.2286 to −0.1416). The same applied to estimated concentrations, but this trend was not significant (ρ = −0.0132; 95% CI −0.1891 to 0.1627). We could not detect a pathobionts proportion threshold: the weak negative effect on lactobacilli was detectable even in the lowest pathobionts proportion categories. By contrast, the mean relative abundance of BV-anaerobes remained stable initially, and only declined when pathobionts made up 30% or more of the VMB. The mean BV-anaerobes estimated concentration significantly increased with increasing estimated pathobionts concentration and did not reach a plateau. In all pathobionts concentration categories, BV-anaerobes outnumbered pathobionts.

TABLE 5

Table 5. VMB correlates by pathobionts relative abundance category, all samples from all studies not influenced by interventions (N = 1,781)^a.

Correlation matrixes for samples not influenced by interventions confirmed that relative abundances of BV-anaerobes and lactobacilli were strongly negatively correlated (ρ = −0.9234; Figure 2A), but showed that their estimated concentrations were not correlated (r = 0.0031; Figure 2B; for correlation coefficients with 95% confidence intervals, see Supplementary Material 2, Table S2). Pathobionts and lactobacilli relative abundances were also negatively correlated, albeit less strongly (ρ = −0.2076), and their estimated concentrations were not (r = 0.0436). Pathobionts and BV-anaerobes relative abundances were not correlated (ρ = 0.0160) and their estimated concentrations were weakly positively correlated (r = 0.1938). Pathobionts also correlated positively with the “other bacteria” rest group (ρ = 0.3831 for relative abundances and r = 0.3388 for estimated concentrations). At individual genus level, the estimated concentrations of the six pathobionts that were a non-minority genus in all three studies correlated positively with one another except for Campylobacter with the other five taxa, and Haemophilus and Escherichia/Shigella.

FIGURE 2

Figure 2. Correlation matrices. (A) Relative abundances, all three studies (N = 1,781 samples)^a. (B) Estimated concentrations, VMB study only (N = 158 samples)^a. BV, bacterial vaginosis.^aOnly includes VMB samples not influenced by the interventions, collected at Screening, Month 6, and Month 1 and 2 in the no-intervention group. Squares indicate the Spearman's rank correlation between −1 and +1. Positive correlations are shown in white and negative correlations in black. The size of each square is loosely proportionate to the magnitude of the correlation. For the actual correlation coefficients, (see Supplementary Material, Tables S2C,D).

Correlates of Vaginal Pathobionts Detection

Finally, we investigated the correlates of pathobionts detection (≥1% vs. <1%), relative abundance, and estimated concentration for all studies combined (Table 6) and for each study separately (Supplementary Material 2, Tables S3A–C). The mean relative abundance and mean estimated concentration of pathobionts decreased with increasing age (except in the HELIUS study), and with ethnicities other than Dutch. The data consistently showed strong associations with Nugent score categories in both studies that assessed these (the Rwanda VMB and HARP studies): the likelihood of detection (OR = 5.29; 95% CI 2.82–9.90), mean relative abundance, and mean estimated concentration of pathobionts were highest for Nugent score category 4–6 (intermediate), followed by 7–10 (BV-positive), and 0–3 (BV-negative). Positive associations between pathobionts detection or levels and hormonal contraceptive use, smoking, antibiotic use in the 14 days prior to sampling, HIV status, and the presence of STI pathogens were found in at least one but not in all studies. Associations with sexual risk-taking and unusual vaginal discharge reporting were inconsistent between studies, and we did not find associations with detection of vaginal yeasts by microscopy.

TABLE 6

Table 6. Correlates of pathobionts detection, relative abundance, and estimated concentration.

Discussion

Seventeen percent of this highly diverse group of women from Africa and Europe had a VMB containing at least 1% pathobionts, and 5.2% had a VMB containing at least 20% pathobionts. Streptococcus was most common (54% of the pathobionts sequencing reads), but Staphylococcus, Enterococcus, Escherichia/Shigella, Haemophilus, and Campylobacter were also detected as non-minority genera in all three studies. Mean relative abundances and estimated concentrations of pathobionts were much lower than those of lactobacilli and BV-anaerobes, but the pathogenic potential may be higher, and these levels may therefore be clinically relevant.

The meta-analysis confirmed that the VMB of many women contains both lactobacilli and BV-anaerobes, but that the BV-anaerobes concentration is low in “healthy” women with lactobacilli-domination. Our relative abundance, estimated concentration, and correlation data may be best explained by the following hypothesis. BV-anaerobes are present or frequently introduced into the vagina of most women, and may start to expand in response to a trigger, such as recent sex or menses (Jespers et al., 2017). When the BV-anaerobes concentration increases, the lactobacilli concentration does not seem to decline much initially, but instead, the total bacterial concentration increases. The lactobacilli relative abundance therefore does decline. We cannot test this hypothesis directly because our analyses were cross-sectional, but the strong negative correlation between lactobacilli and BV-anaerobes relative abundances (ρ = −0.9234) but not estimated concentrations (r = 0.0031), and the higher overall bacterial load of the anaerobic dysbiotic VMB types (means 6.11–6.87 log₁₀ cells/μl) compared to the lactobacilli-dominated VMB types (means 5.13–5.83 log₁₀ cells/μl) fit this hypothesis.

By contrast, a much smaller proportion of women in our study carried pathobionts in their VMB (17% if a 1% relative abundance is used as a cut-off). Our relative abundance, estimated concentration, and correlation data may be best explained by the following hypothesis. Pathobionts are occasionally introduced into the vagina from the gut, urinary tract, and perineum, or from the male partner external genitalia, but are usually cleared or remain at low levels. If they do persist and expand, lactobacilli decline somewhat, and BV-anaerobes expand alongside the pathobionts. As before, we cannot test this hypothesis directly, but the modest positive correlations between estimated concentrations of pathobionts and BV-anaerobes, the declines in estimated concentration/relative abundance of lactobacilli with increasing pathobionts level, and an overall bacterial load of the ≥20% pathobionts VMB type (5.85 log₁₀ cells/μl) that is similar to that of the lactobacilli-dominated VMB types (means 5.13–5.83 log₁₀ cells/μl) fit this hypothesis. Pathobionts also correlated positively with the “other bacteria” group, which contains non-pathogenic skin bacteria such as Corynebacterium. The pathobionts and non-pathogenic skin bacteria may have been introduced into the vagina from the woman's perineum or the skin of the external genitalia of her male partner, but specimen contamination via the hands of specimen handlers cannot be ruled out.

Gram stain Nugent scoring is the current gold standard for BV diagnosis (Nugent et al., 1991). In this method, Gram stained slides are viewed under a microscopy, and three bacterial morphotypes are scored: Gram-positive rods (presumed to be lactobacilli), small Gram-variable rods (presumed to be G. vaginalis), and curved Gram-variable rods (presumed to be Mobiluncus). A Nugent score of 0–3 is considered BV-negative, 4–6 intermediate microbiota, and 7–10 BV-positive. In this meta-analysis, the likelihood of detection, mean relative abundance, and mean estimated concentration of pathobionts were consistently highest for Nugent score category 4–6, followed by 7–10, and 0–3. These findings also fit the above-mentioned hypotheses, and provide a partial explanation for what a Nugent score of 4–6 signifies. A Nugent score of 4–6 should, however, not be used to diagnose pathobionts presence because another significant proportion of these samples likely contain lactobacilli plus BV-anaerobes.

Positive associations between pathobionts detection and/or levels and young age, non-Dutch origin, hormonal contraceptive use, smoking, antibiotic use in the 14 days prior to sampling, HIV status, and the presence of STI pathogens were found in at least one study. All of these factors are also risk factors for anaerobic dysbiosis, except for hormonal contraceptive use. Hormonal contraception, and especially methods containing estrogen, protects women from anaerobic dysbiosis (van de Wijgert et al., 2013). Authors have hypothesized that estrogen increases vaginal glycogen, which is converted into lactic acid by lactobacilli. This keeps BV-anaerobes at bay but perhaps not pathobionts. Streptococci, for example, can tolerate low vaginal pH very well (Shabayek and Spellerberg, 2017). Sexual risk-taking is an important proven risk factor for anaerobic dysbiosis (van de Wijgert et al., 2014), but associations with pathobionts detection and/or levels were inconsistent in this meta-analysis. This could be due to the fact that women in the two African cohorts were recruited based on sexual risk or HIV-status, and women in the Dutch cohort were not, thereby introducing collinearity between study/ethnic group and sexual risk. Associations between pathobionts detection/levels and unusual vaginal discharge reporting were also inconsistent between studies. In our experience, vaginal symptom-reporting rarely correlates well with the actual presence of a vaginal infection or vaginal dysbiosis (Verwijs et al., 2019b). None of the women in the three studies had severe symptoms, such as those associated with desquamative inflammatory vaginitis, and we could therefore not test the association between pathobionts levels and such symptoms.

A limitation of our study is that each of the three studies used slightly different sequencing-related laboratory and initial data processing methods (Table 1). However, we took this into account by stratifying most analyses by study. Another limitation is that we did not detect and quantify individual pathobiont species and genera by quantitative PCR. Past S. agalactiae prevalence studies using selective culture have shown average rectovaginal detection rates of 22% in sub-Saharan Africa and 19% in Europe (Kwatra et al., 2016), and a recent study using quantitative PCR found a vaginal detection rate of 20% in Kenyan women and 23% in South African women (Cools et al., 2016). Detection of vaginal Streptococcus was lower in our meta-analysis when a relative abundance cut-off of 1% was used (12.7% of all women at baseline; about a quarter of those were S. agalactiae). It is currently not known how rectovaginal and vaginal selective cultures, PCR, and sequencing results relate to one another, but it is possible that pathobionts are not only under-detected in sequencing studies due to the bioinformatics used, but also due to DNA extraction, amplification, and other biases. For example, the detection rate in the HARP study was especially low, which may have been due to the fact that we did not use bead-beating during DNA extraction in that study (Gill et al., 2016). Third, some of the standard diagnostic tests that we used are known to have lower sensitivity than NAAT-based tests (e.g., culture for vulvovaginal candidiasis and T. vaginalis), and not all women in all three studies were screened for all STI pathogens.

We also report some limitations related to our statistical analyses. Correlating variables derived from relative abundance data is problematic because they are not independent (Knight et al., 2018); estimated concentrations of these same variables are independent and did indeed provide new insights as described above. However, we only had estimated concentration data for the Rwanda VMB study. Furthermore, our analyses were cross-sectional, and some of them had limited statistical power. Our findings are therefore hypothesis-generating and the hypotheses should be tested in well-powered longitudinal studies that assess the VMB quantitatively. We did not exclude all women who had recently used antibiotics, but reported antibiotic use in the last 2 weeks was rare. A strength of our study is the inclusion of women and samples from three world regions and multiple ethnic groups, and with different behaviors and STI pathogen exposures. The variability in VMB compositions that we observed reflects this wide variety of study participants.

Conclusion

While substantial presence of pathobionts in the VMB was less common than anaerobic dysbiosis, the pathogenic potential of pathobionts is higher than that of BV-anaerobes, and modest levels could therefore be clinically relevant. The most frequently used VMB types, and analyses limited to relative abundance, are inadequate. We recommend that future etiologic and intervention studies quantify the most common vaginal pathobiont genera, as well as lactobacilli and BV-anaerobes. Furthermore, the various detection and quantification methods (culture, PCR, and sequencing) should be rigorously compared to one another to facilitate interpretation of clinical study results.

Data Availability Statement

The three studies included in this meta-analyses were governed by different institutions and ethics committees. Data availability therefore differs for each study. The three original publications, which are referenced in this publication, include data availability statements for each respective study. Additional data unique to this publication are provided in Supplementary Material 1.

Ethics Statement

The Rwanda VMB study was sponsored by the University of Liverpool, approved by the Rwanda National Ethics Committee and the University of Liverpool Research Ethics Subcommittee for Physical Interventions, and registered on ClinicalTrials.gov (NCT02459665). The South African portion of the HARP study was approved by the ethics committees of the University of Witwatersrand in Johannesburg and the London School of Hygiene and Tropical Medicine. The HELIUS study was approved by the Medical Ethics Committee of the Academic Medical Center in Amsterdam (protocol number: 10/100; amendment 10/100# 10.17.1729; NL32251.018.10). The patients/participants provided their written informed consent to participate in this study.

Author Contributions

JW was the Principal Investigator of all three VMB sequencing studies. PM was the Principal Investigator of the HARP parent study. MV, AG, HB, and CV contributed to the VMB sequencing laboratory work, processed the initial sequencing data, and compiled the metadata datasets. MV and JW wrote the data analysis plan. MV compiled the combined dataset and conducted the analyses presented in this paper. JW and MV wrote the manuscript. All authors commented on and approved the manuscript.

Funding

The Rwanda VMB study was supported by the DFID/MRC/Wellcome Trust Joint Global Health Trials Scheme as a Development Project [grant reference MR/M017443/1; grant title: Preparing for a clinical trial of interventions to maintain normal vaginal microbiota for preventing adverse reproductive health outcomes in Africa] and the University of Liverpool [institutional funding for the 16S rRNA sequencing]. The HARP study was funded by the European Commission 7th Framework Programme [original HARP study; grant agreement number HEALTH-2010-F2-265396] and the University of Liverpool [institutional funding for the 16S rRNA sequencing]. The HELIUS study is conducted by the Academic Medical Center Amsterdam and the Public Health Service of Amsterdam. Both organizations provided core support for HELIUS [institutional funding]. The HELIUS study is also funded by the Dutch Heart Foundation [grant number 2010T084], the Netherlands Organization for Health Research and Development (ZonMw) [grant number 200500003], the European Union (FP-7) [grant number 278901], and the European Fund for the Integration of non-EU immigrants (EIF) [grant number 278901]. The HELIUS VMB Study reported here was additionally supported by the Aids Fonds Netherlands [grant number 201102] and the Netherlands Organization for Health Research and Development (ZonMw) [grant numbers 7115 0001 and 204005002]. Some HELIUS analyses were carried out on the Dutch national e-infrastructure with the support of SURF Foundation. The views expressed in this paper are those of the authors and do not necessarily represent the views of the authors' institutions or the funders.

Conflict of Interest

The 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.

Acknowledgments

The authors thank all study participants and members of the three study groups. These consisted of:

Rwanda VMB Study Group: S. K. Agaba, C. Bronowski, A. C. Darby, E. Lievens, L. Mwambarangwe, A. Nivoliez, J. Ravel, M. Uwineza, M. C. Verwijs, J. H. H. M. van de Wijgert, and the other staff of Rinda Ubuzima in Kigali, the National Reference Laboratory in Kigali, and the Center for Genomic Research at the University of Liverpool.

HARP Study Group: A. Chikandiwa, E. Cutler, S. Delany-Moretlwe, L. Diep, C. Hwinya, D. A. Lewis, M. P. Magooa, V. Maseko, P. Michelow, B. Muzah, T. Omar, A. Puren (Johannesburg, South Africa); F. Djigma, J. Drabo, O. Goumbri-Lompo, D. Karambiri, N. Meda, B. Sawadogo, J. Simporé, R. Sombie, A. Yonli, Z. Zamane, S. Zan (Ouagadougou, Burkina Faso); V. Costes, M. N. Didelot, S. Doutre, N. Leventoux, N. Nagot, J. Ngou, M. Segondy (Montpellier, France); and A. Devine, C. Gilham, L. Gibson, H. Kelly, R. Legood, P. Mayaud, H. A. Weiss (London, UK). The HARP-VMB Study Group consisted of: A. C. Darby, S. C. Francis, A. C. Gill, P. Mayaud, J. H. H. M. van de Wijgert, and other staff of the Center for Genomic Research at the University of Liverpool.

HELIUS VMB Study Group: C. J. Alberts, H. Borgdorff, S. M. Bruisten, S. E. Geerlings, R. van Houdt, M. Prins, M. F. Schim van der Loeff, M. B. Snijder, C. van der Veer, H. J. de Vries, and J. H. H. M. van de Wijgert.

Supplementary Material

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

Supplementary Material 1. List of taxa identified in the three studies and classification into bacterial groups.

References

Black, C. G., Tavares, L., Stachel, A., Ratner, A. J., and Randis, T. M. (2018). Distribution of late-onset neonatal sepsis pathogens differs in inpatient and outpatient settings. Amer J Perinatol. 36, 1136–1141. doi: 10.1055/s-0038-1675643

PubMed Abstract | CrossRef Full Text | Google Scholar

Borgdorff, H., Veer, C., van der Houdt, R., van Alberts, C. J., Vries, H. J., de Bruisten, S. M., et al. (2017). The association between ethnicity and vaginal microbiota composition in Amsterdam, the Netherlands. PLoS ONE 12:e0181135. doi: 10.1371/journal.pone.0181135

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunham, R. C., Gottlieb, S. L., and Paavonen, J. (2015). Pelvic inflammatory disease. N. Engl. J. Med. 372, 2039–2048. doi: 10.1056/NEJMra1411426

PubMed Abstract | CrossRef Full Text | Google Scholar

Cools, P., Jespers, V., Hardy, L., Crucitti, T., Delany-Moretlwe, S., Mwaura, M., et al. (2016). A multi-country cross-sectional study of vaginal carriage of Group B streptococci (GBS) and Escherichia coli in resource-poor settings: prevalences and risk factors. PLoS ONE 11:e0148052. doi: 10.1371/journal.pone.0148052

PubMed Abstract | CrossRef Full Text | Google Scholar

Donders, G. G. G., Bellen, G., Grinceviciene, S., Ruban, K., and Vieira-Baptista, P. (2017). Aerobic vaginitis: no longer a stranger. Res. Microbiol. 168, 845–858. doi: 10.1016/j.resmic.2017.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Fadrosh, D. W., Ma, B., Gajer, P., Sengamalay, N., Ott, S., Brotman, R. M., et al. (2014). An improved dual-indexing approach for multiplexed 16S rRNA gene sequencing on the Illumina MiSeq platform. Microbiome 2:6. doi: 10.1186/2049-2618-2-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Gill, C., van de Wijgert, J. H. H. M., Blow, F., and Darby, A. C. (2016). Evaluation of lysis methods for the extraction of bacterial DNA for analysis of the vaginal microbiota. PLoS ONE 11:e0163148. doi: 10.1371/journal.pone.0163148

PubMed Abstract | CrossRef Full Text | Google Scholar

Jespers, V., Kyongo, J., Joseph, S., Hardy, L., Cools, P., Crucitti, T., et al. (2017). A longitudinal analysis of the vaginal microbiota and vaginal immune mediators in women from sub-Saharan Africa. Sci. Rep. 7:11974. doi: 10.1038/s41598-017-12198-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Knight, R., Vrbanac, A., Taylor, B. C., Aksenov, A., Callewaert, C., Debelius, J., et al. (2018). Best practices for analysing microbiomes. Nat. Rev. Microbiol. 16, 410–422. doi: 10.1038/s41579-018-0029-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Kwatra, G., Cunnington, M. C., Merrall, E., Adrian, P. V., Ip, M., Klugman, K. P., et al. (2016). Prevalence of maternal colonisation with group B streptococcus: a systematic review and meta-analysis. Lancet Infect. Dis. 16, 1076–1084. doi: 10.1016/S1473-3099(16)30055-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, C. M., Aziz, M., Kachur, S., Hsueh, P.-R., Huang, Y.-T., Keim, P., et al. (2012). BactQuant: an enhanced broad-coverage bacterial quantitative real-time PCR assay. BMC Microbiol. 12:56. doi: 10.1186/1471-2180-12-56

PubMed Abstract | CrossRef Full Text | Google Scholar

Nugent, R. P., Krohn, M. A., and Hillier, S. L. (1991). Reliability of diagnosing bacterial vaginosis is improved by a standardized method of gram stain interpretation. J. Clin. Microbiol. 29, 297–301. doi: 10.1128/JCM.29.2.297-301.1991

PubMed Abstract | CrossRef Full Text | Google Scholar

Oerlemans, E. F. M. (2019). Topical Lactobacillus Applications for Modulation of the Vaginal and Skin Microbiota. Available online at: https://hdl.handle.net/10067/1583510151162165141 (accessed November 30, 2019).

Google Scholar

Paavonen, J., and Brunham, R. C. (2018). Bacterial vaginosis and desquamative inflammatory vaginitis. New Engl. J. Med. 379, 2246–2254. doi: 10.1056/NEJMra1808418

PubMed Abstract | CrossRef Full Text | Google Scholar

Pereira, N., Edlind, T. D., Schlievert, P. M., and Nyirjesy, P. (2013). Vaginal toxic shock reaction triggering desquamative inflammatory vaginitis. J. Low. Genit. Tract Dis. 17, 88–91. doi: 10.1097/LGT.0b013e3182656991

PubMed Abstract | CrossRef Full Text | Google Scholar

Ravel, J., Gajer, P., Abdo, Z., Schneider, G. M., Koenig, S. S. K., McCulle, S. L., et al. (2011). Vaginal microbiome of reproductive-age women. Proc. Natl. Acad. Sci. U.S.A. 108(Suppl. 1), 4680–4687. doi: 10.1073/pnas.1002611107

PubMed Abstract | CrossRef Full Text | Google Scholar

Shabayek, S., and Spellerberg, B. (2017). Acid stress response mechanisms of group B streptococci. Front. Cell. Infect. Microbiol. 7:395. doi: 10.3389/fcimb.2017.00395

PubMed Abstract | CrossRef Full Text | Google Scholar

Snijder, M. B., Galenkamp, H., Prins, M., Derks, E. M., Peters, R. J. G., Zwinderman, A. H., et al. (2017). Cohort profile: the Healthy Life in an Urban Setting (HELIUS) study in Amsterdam, the Netherlands. BMJ Open 7:e017873. doi: 10.1136/bmjopen-2017-017873

PubMed Abstract | CrossRef Full Text | Google Scholar

Sobel, J. D. (1994). Desquamative inflammatory vaginitis: a new subgroup of purulent vaginitis responsive to topical 2% clindamycin therapy. Am. J. Obstet. Gynecol. 171, 1215–1220. doi: 10.1016/0002-9378(94)90135-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Sobel, J. D., Reichman, O., Misra, D., and Yoo, W. (2011). Prognosis and treatment of desquamative inflammatory vaginitis. Obstet. Gynecol. 117, 850–855. doi: 10.1097/AOG.0b013e3182117c9e

PubMed Abstract | CrossRef Full Text | Google Scholar

Taiyun (2019). Package Corrplot is for Visualizing a Correlation Matrix: Taiyun/Corrplot. Available online at: https://github.com/taiyun/corrplot (accessed July 9, 2019).

Google Scholar

van de Wijgert, J. H. H. M., Borgdorff, H., Verhelst, R., Crucitti, T., Francis, S., Verstraelen, H., et al. (2014). The vaginal microbiota: what have we learned after a decade of molecular characterization? PLoS ONE 9:e105998. doi: 10.1371/journal.pone.0105998

PubMed Abstract | CrossRef Full Text | Google Scholar

van de Wijgert, J. H. H. M., and Jespers, V. (2017). The global health impact of vaginal dysbiosis. Res. Microbiol. 168, 859–864. doi: 10.1016/j.resmic.2017.02.003

PubMed Abstract | CrossRef Full Text | Google Scholar

van de Wijgert, J. H. H. M., Verwijs, M. C., Turner, A. N., and Morrison, C. S. (2013). Hormonal contraception decreases bacterial vaginosis but oral contraception may increase candidiasis: implications for HIV transmission. AIDS 27, 2141–2153. doi: 10.1097/QAD.0b013e32836290b6

PubMed Abstract | CrossRef Full Text | Google Scholar

van de Wijgert, J. H. H. M., Gill, A. C., Chikandiwa, A., Verwijs, M. C., Kelly, H. A., Omar, T., et al. (2020b). Human papillomavirus infection and cervical dysplasia in HIV-positive women: potential role of the vaginal microbiota. AIDS 34, 115–125. doi: 10.1097/QAD.0000000000002381

PubMed Abstract | CrossRef Full Text | Google Scholar

van de Wijgert, J. H. H. M., Verwijs, M. C., Agaba, S. K., Bronowski, C., Mwambarangwe, L., Uwineza, M., et al. (2020a). Intermittent lactobacilli-containing vaginal probiotic or metronidazole use to prevent bacterial vaginosis recurrence: a pilot study incorporating microscopy and sequencing. Sci. Rep. 3:3884. doi: 10.1038/s41598-020-60671-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Verwijs, M. C., Agaba, S. K., Darby, A. C., and van de Wijgert, J. H. H. M. (2019a). Impact of oral metronidazole treatment on the vaginal microbiota and correlates of treatment failure. Am J Obstet Gynecol. 222, 157.e1–157.e13. doi: 10.1016/j.ajog.2019.08.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Verwijs, M. C., Agaba, S. K., Sumanyi, J.-C., Umulisa, M. M., Mwambarangwe, L., Musengamana, V., et al. (2019b). Targeted point-of-care testing compared with syndromic management of urogenital infections in women (WISH): a cross-sectional screening and diagnostic accuracy study. Lancet Infect. Dis. 19, 658–669. doi: 10.1016/S1473-3099(18)30724-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Warnes, G. R., Bolker, B., Bonebakker, L., Gentleman, R., Liaw, W. H. A., Lumley, T., et al. (2016). gplots: Various R Programming Tools For Plotting Data. Available online at: https://CRAN.R-project.org/package=gplots (accessed September 12, 2018).

Google Scholar