Longitudinal Changes in the Composition of the Penile Microbiome Are Associated With Circumcision Status, HIV and HSV-2 Status, Sexual Practices, and Female Partner Microbiome Composition

Introduction

Randomized controlled trials demonstrated that voluntary medical male circumcision reduces penile anaerobic bacteria (Price et al., 2010; Liu et al., 2013), and this is a potential mechanism by which HIV risk is reduced (Liu et al., 2017). Subsequent studies demonstrate that the penile microbiome is associated with bacterial vaginosis (BV) (Eren et al., 2011; Zozaya et al., 2016; Mehta et al., 2020) and HSV-2 (Mehta et al., 2020) in female sex partners and with human papillomavirus and HIV (Onywera et al., 2020a) infections in men. Despite the potential contribution of the penile microbiome to sexual health outcomes in men and their female sex partners’ risk of BV and HSV-2, there is a limited description of the penile microbiome over time and factors associated with the longitudinal variation. Most studies have been cross-sectional or focused primarily on the role of circumcision status and have a short duration. We identified only one non-interventional, longitudinal study of the penile microbiome in the peer-reviewed literature. In a U.S. sample of 18 male adolescent participants aged 14–17 years (5 circumcised and 10 having had vaginal sex), Nelson et al. observed that several taxa with the highest mean relative abundance (RA) recovered from the coronal sulcus were stable across four measurements taken monthly over a 3-month period (Nelson et al., 2012).

In our analysis of 168 heterosexual couples in Kenya, we observed that baseline penile meatus microbiome composition predicted incident BV in female partners with high accuracy, up to 6 to 12 months after penile microbiome measure (Mehta et al., 2020). We posited that one explanation for the prediction’s durability over many months may be due to the temporal stability of the penile microbiome. The goals of this analysis were to 1) characterize penile microbiome composition over a 1-year period and 2) identify factors associated with penile microbiome composition over time. Given the interrelatedness of the penile microbiome and vaginal microbiome (VMB) in relation to several sexual and reproductive health outcomes, such longitudinal studies can yield value to identify epidemiological associations that may lead to improved sexual health outcomes for men and their female sex partners.

Materials and Methods

This study was approved by the Institutional Review Board of the University of Chicago, IL, USA (2013-0511) and the Ethical Review Committee of the Maseno University of Kisumu, Kenya (MSU/DRPC/MUERC/00054/13).

Study Design and Participants

This study used data and biological specimens from a prospective cohort study of community-recruited heterosexual couples in Kisumu, Kenya. The details of recruitment have been previously described (Mehta et al., 2018). Briefly, cisgender male (aged 18–35 years) and female (aged 16 years and older) couples who had been in a sexual relationship for at least 6 months were recruited. The study was powered to estimate the effect of the penile microbiome on female sex partners’ risk of BV, with parameters of sample size estimation detailed in (Mehta et al., 2018).

Trained clinicians and counselors obtained standardized medical history, physical examination, and personal interviews to obtain information on sociodemographics and sexual behavior. Laboratory methods have been described in detail previously. Briefly, serum was tested for HSV-2 antibody (Kalon HSV-2 IgG ELISA, Kalon Biological Limited, Aldershot, UK) using the manufacturer’s recommended cutoff, a clinician-collected vaginal swab was used to assess BV according to Nugent’s criteria (Nugent et al., 1991), and testing for HIV infection was conducted using a serial rapid test protocol that followed the Kenyan national guidelines (National and Control Programme, 2015).

Penile and Vaginal Microbiome Characterization

Clinicians twirled pre-moistened mini-tip flocked swabs (Copan Diagnostics, Inc., Corona, CA, USA) at the meatal opening for 3–5 rotations (Mehta et al., 2020). In uncircumcised men, the foreskin was retracted prior to sampling. The VMB was characterized in cervicovaginal lavage samples. Penile swabs and cervicovaginal lavage specimens were stored at −80°C until shipment. PCR amplification of the V3–V4 variable region of bacterial 16S rRNA genes was performed using a two-step PCR protocol with primers 341F and 806R, as described previously (Naqib et al., 2018). Amplicons were sequenced on an Illumina MiSeq instrument, implementing V3 chemistry (600 cycles). Forward and reverse reads were merged using the software package PEAR (Zhang et al., 2014). After standard read processing (i.e., quality trimming, primer removal, and chimera removal), annotation was conducted by the University of Maryland Institute for Genomic Science (Holm et al., 2019), as performed previously (Mehta et al., 2020). Subsequently, a biological observation matrix was generated at the lowest taxonomic level identifiable. Taxa identified as contaminants in negative controls [via decontam, (1.16.0) (Davis et al., 2018) in R] were removed, and samples with <5,000 reads (n = 19) were excluded. Data were filtered to retain taxa that contributed at least 0.1% of the total sequence reads, resulting in the retention of 54 penile taxa (from 2,138 taxa in the initial observation matrix). Raw sequence data files are available in the Sequence Read Archive (National Center for Biotechnology Information; BioProject identifier PRJNA 516684). Vaginal community state types (CSTs) were identified in a reference dataset using the nearest centroid classifier approach, VALENCIA (VAginaL community state typE Nearest CentroId clAssifier) (France et al., 2020). Of 252 men enrolled in the study, 249 men had microbiome data available. Participants with one visit were excluded from this analysis, resulting in 218 men contributing 740 observations: 47 with 2 visits, 58 with 3 visits, and 118 with 4 visits.

Statistical Analysis

Inferential Analyses

Explanatory variables that were a priori hypothesized to be associated with the penile microbiome composition and change over time were tested: circumcision status, HIV status, HSV-2 status, condom use, number of sex partners, days since last sex, age, socioeconomic status (employment, educational attainment), female sex partner BV status, vaginal CST, and HIV and HSV-2 status.

First, we used permutational ANOVA (PERMANOVA) to test for differences in penile bacterial community composition between study time points using pairwise comparisons with correction for multiple testing (RVAideMemoire (0.9-81-2) (Herve, [[NoYear]]), in R) of the Bray–Curtis resemblance matrix. Next, hierarchical clustering was conducted to identify penile microbiome community types (CTs). Following methods recommended for reliable clustering (Koren et al., 2013), partition was applied around the medoid using five distance measures (Euclidean, Bray–Curtis, Mahalanobis, Jaccard, and Kulczynski). The optimum number of clusters that maximized both silhouette index (silhouette function (Desgraupes, 2017) in cluster (2.1.2) R package (Maechler et al., 2022), a measure of similarity of samples within their own cluster compared to other clusters) and prediction strength (prediction.strength fpc (2.2-9) (Hennig, [[NoYear]]), a validation test to check if the clusters are replicable in training and test sets) was selected. The optimal number of clusters based on silhouette index and prediction score was 7, using the Euclidean distance of the natural log-transformed RAs. Multinomial regression to identify changes in each of the 7 penile CTs over time and associated factors were used, with a random effect for subject intercept. Third, a linear mixed-effects model was used to test whether the RA (outcome) of the twelve meatal taxa with the highest mean RA was impacted by time and other subject-specific characteristics. Following geometric Bayesian multiplicative prior imputation of zeros [zCompositions (1.3.4) (Palarea-Albaladejo and Martin-Fernandez, 2015)], separate models were fitted with each center log ratio-transformed taxon as an outcome with a random effect for individual, and covariates were treated as fixed effects [lme4 (1.1-27.1)] (Bates et al., 2015). Step-down Akaike information criterion (AIC) (Posada and Buckley, 2004) guided the final model selection. With the use of this same linear mixed model approach, factors associated with alpha diversity metrics over time were examined. Data were rarefied to 6,000 sequences prior to estimation of alpha diversity measures (vegan (2.5-7) (Oksanen et al., 2022)). Measures of Shannon, Simpson, and Evenness were highly correlated (0.94–0.96), and therefore, Shannon diversity index and richness were modeled due to their lower correlation (0.48–0.73) with the other alpha diversity measures. In all models, interactions with circumcision status and time, and with BV and time, were tested as an a priori hypothesis for association with stability and change. Fourth, group-based trajectories (GBTs) were estimated with zero-inflated Poisson regression for the four taxa with the highest RA (Corynebacterium, Streptococcus, Staphylococcus, and Anaerococcus) (Nagin and Odgers, 2010; Jones and Nagin, 2013). These GBTs shed light on longitudinal trends, by delineating stable and changing trajectories of RA. The number of groups and polynomial type for each trajectory was selected based on minimized AIC, although if an outlier group (<5% of observations) was produced, the number of groups was reduced. The final model selection for outcomes of CT, individual taxa, and alpha diversity measures was guided by minimizing AIC (Posada and Buckley, 2004). Lastly, the within-subject (i.e., paired) Bray–Curtis similarity was estimated from baseline to 6 months and from 6 months to 12 months, and multinomial mixed-effects modeling was used to assess factors associated with the Bray–Curtis similarity over time. The Bray–Curtis similarity was categorized as follows: 1) “less stable” included observations in the lowest 15th percentile and corresponded to <41% of Bray–Curtis similarity in paired measurements, 2) “stable” included observations between the 15th and 85th percentiles and corresponded to 41% to <70% similarity, and 3) “more stable” included observations in the highest 15th percentile of similarity and corresponded to >70% similarity. Final model selection for outcomes of CT, individual taxa, alpha diversity measures, and Bray–Curtis similarity was guided by minimizing AIC (Posada and Buckley, 2004).

Visualizations

Stacked bar charts, coefficient plots, and GBT plots were generated in Stata/SE 17. Change in the RA of individual taxa and alpha diversity measures across time by circumcision status and female partner BV status were visualized using line plots in R and visualization of transition of subjects across clusters with Sankey diagrams (ggalluvial (0.12.3) (Brunson and Read, 2020) in ggplot2 (3.3.6) (Wickham, 2016)). The relationship of global bacterial communities was visualized by study time point using metric multidimensional scaling of bootstrapped averages of centroids with 100 replicates per group (Primer-E, version 7, UK). A beta coefficient heatmap was created in R using ggplot function within ggplot2 (Wickham, 2016).

Results

At baseline, men had a median age of 26 years, 11.6% were living with HIV, 46.1% were HSV-2 seropositive, and 55.5% were circumcised (Table 1). Female sex partners have a median age of 22 years, 10.2% were living with HIV, 55.4% were HSV-2 seropositive, 20.2% were with BV, and the majority (89.2%) of CSTs were either diverse (CST-IV, 46.5%) or Lactobacillus iners dominant (CST-III, 43.7%).

TABLE 1

Table 1 Distribution of characteristics by study visit, N = 740 observations.

Penile Community Type

We identified seven penile CTs (Figure 1, Table 2): mixed (CT-1), Corynebacterium dominated (CT-2), Streptococcus dominated (CT-3), Sneathia sanguinegens dominated (CT-4), Finegoldia/Anaerococcus dominated (CT-5), L. iners dominated (CT-6), and Clostridiales Family XI/Prevotella timonensis dominated (CT-7). Although clustering was unsupervised, clusters were delineated by circumcision status, emphasizing the role of circumcision as a driver of penile microbiome composition: just 7.7% and 4.6% of men in CT-5 and CT-7 were circumcised, respectively, compared to much higher proportions (60%–85.7%) in other CTs. CT-7 had the highest proportion of female partners with BV (36.4%) and with HIV (18.2%), the highest proportion of never using condoms (86.4%), and the lowest proportion of men with at least some secondary educational attainment (36.4%). By comparison, men in CT-1 had greater educational attainment (72.3% secondary or higher), the lowest prevalence of HSV-2 (35.4%), and a moderate prevalence of BV among female partners (21.2%). Compared to other CTs, for CT-2, CT-3, and CT-6, alpha diversity was decreased, BV prevalence among female partners was decreased, and the proportion of female partners with CST-I (Lactobacillus crispatus dominated) was increased.

FIGURE 1

Figure 1 Stacked bar charts of the relative abundance of the penile taxa with greatest mean relative abundance for individuals by community type. Legend: observations from 740 samples are sorted by penile community type. The relative abundance of the taxa with the highest mean relative abundance is shown (0%–100%; y-axis), with individual subjects represented by individual bars, separated by community type (CT; x-axis). The bar at the top represents circumcision status, with circumcised individuals represented by red vertical lines.

TABLE 2

Table 2 Distribution of baseline male and female partner characteristics by penile community type (CT).

Factors Associated With Penile Community Type Over Time

There were statistically significant (p < 0.05) changes over time from baseline and 1-month visits to 6- and 12-month visits (Figure 2). Additionally, alluvial visualization indicates that individuals changed meatal CT over time (Figure 3), with an increasing proportion of individuals in CT-1. In multivariable-adjusted analyses (Table 3), relative to CT-2 (Corynebacterium dominated), CT-1 (mixed) increased over time relative to baseline, with transient increases in CT-5 (Finegoldia dominated) and CT-6 (L. iners dominated) at 6 months. Relative to CT-2, men were less likely to be circumcised if they were in CT-4 (S. sanguinegens dominated; adjusted odds ratio (aOR) = 0.45), CT-5 (aOR = 0.04), CT-6 (aOR = 0.28), and CT-7 (Clostridiales Family XI/P. timonensis; aOR = 0.03). Men who reported using condoms at the last sex act were less likely to have penile CT-4 than CT-2 (aOR = 0.43). Men who reported washing their penis more than 1 h after sex were more likely to have CT-5 (aOR = 2.39), CT-6 (aOR = 2.61), or CT-7 (aOR = 2.31) penile microbiome than CT-2. The association of vaginal CST-IV with penile CT-7 (aOR = 5.30) remained statistically significant in adjusted analyses. Univariate results were largely similar in terms of the magnitude of associations and significance (Supplementary Table 1).

FIGURE 2

Figure 2 Metric multidimensional scaling plot of penile meatal microbial community by study time point. Legend: the four different colors represent the four study time points: red = baseline; yellow = 1 month; green = 6 months; blue = 12 months. Each colored mark indicates one of 100 bootstraps of the dataset. The matching shaded area represents 95% coverage. The black symbol at the center of each colored shape represents the average centroid of the 100 bootstraps. The false discovery rate-corrected p-values of the pairwise comparisons by time point are shown below the figure.

FIGURE 3

Figure 3 Sankey diagram showing change in penile meatal community type (CT) over time. Legend: the percent of samples (y-axis) classified in seven penile meatal community types (CT) over time (x-axis) are represented by seven colored blocks, with key to the right.

TABLE 3

Table 3 Results of multivariable mixed effects multinomial modeling*: adjusted^ odd ratio and 95% CI of penile meatal community type.

Change in Individual Penile Taxa Over Time and Associated Factors

The twelve taxa with the highest mean RA accounted for 77.0% of the overall mean RA, with Corynebacterium (16.2%), Anaerococcus (9.4%), Streptococcus (8.5%), Finegoldia (8.0%), and L. iners (7.4%) having the greatest mean RA (Figure 4). Results of multivariable models are summarized visually (Figure 5; Supplementary Table 2 for tabular details). On average, changes over time were not statistically significant for most taxa; only the RAs of Corynebacterium and Anaerococcus changed, both declining over time. Men who were circumcised had higher RA of Corynebacterium, Streptococcus, and Staphylococcus and lower RAs of several anaerobic bacteria. Unsurprisingly, CST-III (L. iners dominant) in female sex partners had a strong positive association with the RA of penile L. iners. Compared to men whose sex partners had CST-I (L. crispatus dominant), those whose partners had CST-III and CST-IV (Gardnerella vaginalis dominant) had decreased RA of Staphylococcus and increased RA of anaerobes.

FIGURE 4

Figure 4 Mean relative abundance over time for the twelve penile meatal taxa with highest mean relative abundance, by circumcision status. Legend: the mean relative abundance of the twelve penile meatal taxa with highest mean relative abundance (RA; y-axis) are shown over time (x-axis) with the 95% CI shown in the shaded areas around the mean line. Mean relative abundances over time are stratified by circumcision status, as represented by legend to the right of the graphs.

FIGURE 5

Figure 5 Coefficient heatmap showing results of multivariable-adjusted random-effects regression: factors associated with 12 taxa with highest mean relative abundance over time. Legend: the columns represent the 12 taxa with the highest mean relative abundance. The rows represent the covariates associated with taxa in multivariable models that are time-adjusted, with final models selected according to minimized Akaike’s information criterion. A red circle represents a positive association between covariate and taxon, and a blue circle represents a negative (inverse) association. The size of the dot is proportional to the absolute magnitude of the coefficient. White space indicates that the association between covariate and taxon did not improve explanation and is not included in the final model.

Trajectories of Taxa Over Time and Associated Factors

We identified three trajectories for each of the four taxa with the highest mean RA (Figure 6). Within Corynebacterium, 40.3% of men exhibited a stable, low RA trajectory and 35.2% with a stable, medium RA trajectory. Approximately a quarter of men (24.5%) were classified into a high RA trajectory that decreased over time. Staphylococcus and Streptococcus trajectories were similar, with a majority of men (59.7% and 59.3%, respectively) having a stable, low RA trajectory; the proportions of men with medium and high RA trajectories were also similar, which is in keeping with the correlation between Staphylococcus and Streptococcus RA (Supplementary Figure 1). Although Anaerococcus also had “low,” “medium,” and “high” RA trajectory groups, the low RA trajectory (27.6% of men) showed substantially increasing RA over time, and the medium RA (47.3%) and high RA (25.1% of men) groups showed decreasing RA over time. In general, our models appear to identify factors associated with composition rather than changing trajectory. For example, men grouped into the high, decreasing RA trajectory of Corynebacterium were more likely to be circumcised (aOR = 20.4) and less likely to be HSV-2 seropositive (aOR = 0.42), which are two factors with few time-varying values. Aside from circumcision status, models failed to identify factors associated with Staphylococcus and Anaerococcus trajectories.

FIGURE 6

Figure 6 Group-based trajectories and results of multivariable multinomial modeling of factors associated with trajectory. Legend: the figures show the estimated group trajectories of relative abundances of the four taxa with highest relative abundance: (A) Corynebacterium, (B) Streptococcus, (C) Staphylococcus, and (D) Anaerococcus. The y-axis represents relative abundance (RA), and the x-axis represents time in months, with study visits indicated. Each trajectory is represented by a different color (blue, red, or green) with group membership shown in the legend below each figure. For example, for Corynebacterium, there is a stable trajectory with low relative abundance (dark blue), and 40.3% of observations are grouped in this trajectory. There is a stable, “medium” relative abundance (dark red), and 35.2% of observations are grouped in this trajectory. There is a (dark green) “high” relative abundance trajectory (with decline over time), and 24.5% of observations are grouped in this trajectory. The tables to the right of each figure (i–iv) show the results of multivariable multinomial modeling of baseline factors associated with group trajectories. For Corynebacterium, compared to men who are grouped in the stable, low RA group trajectory, men who are grouped in the stable, “medium” RA trajectory are more likely to be circumcised (adjusted odds ratio 7.17) and less likely to be HSV-2 seropositive (adjusted odds ratio = 0.37). Men who are grouped in the decreasing, “high” RA trajectory are more likely to be circumcised (adjusted odds ratio 20.4), less likely to be HSV-2 seropositive, more likely to be HIV positive, and more likely to report condom use at the last sexual encounter.

Change in Alpha Diversity Measures Over Time and Associated Factors

In univariate analysis, alpha diversity measures did not statistically significantly vary over time (Supplementary Figure 2, Supplementary Table 3). In adjusted mixed-effects modeling (Figure 7, Supplementary Table 3), variables showed similar patterns of association with Shannon diversity and richness. Circumcised men had decreased alpha diversity measures, as did those who reported condom use, and a longer time since last sexual intercourse. Employed men and those whose female partners had BV or a VMB classified as CST-III or CST-IV (vs. CST-I) had increased alpha diversity measures. Shannon diversity index did not change over time, and there were no statistically significant or meaningful interactions of circumcision status or female partner BV with time. A statistically significant increase in richness (coefficient = 1.96) was observed at the 12-month visit in the adjusted model: this was due to an interaction over time with the female partner’s BV status. Although overall penile microbiome richness was greater for men whose partners had BV, there was a reduction in penile richness at the 12-month visit for men whose partners had BV and an increase in richness for men whose partners did not have BV; a similar but non-significant trend was observed for Shannon diversity index (Supplementary Figure 3).

FIGURE 7

Figure 7 Results of multivariable adjusted mixed-effects linear regression: factors associated with penile meatal microbiome alpha diversity measures over time: (A) Shannon diversity index and (B) richness. Legend: the coefficient plots reflect the coefficient (white diamond) placed over the 95% CI (gray-shaded bars) for each factor (y-axis) in relation to the alpha diversity measures: (A) Shannon diversity index and (B) richness. p-Values are indicated under each variable. N = 654 observations among 221 individuals in both models.

Paired Measures of Bray–Curtis Similarity Over Time

The mean Bray–Curtis similarity was 56.8% for paired observations between 0 and 6 months and 56.3% between 6 and 12 months. In multivariable mixed-effects regression (Supplementary Table 4), men with less stable (<41%) Bray–Curtis similarity were more likely to have HSV-2 (aOR = 2.99) and less likely to be currently employed (aOR = 0.42). Only two variables were marginally associated with higher stability (>71%): increasing age was inversely associated (aOR = 0.92; 95% CI: 0.84–1.01), as was more frequent condom use (aOR = 0.50; 95% CI: 0.23–1.11). Circumcision status, female partner BV status and VMB CST, HIV status (male or female), and sexual practices were not associated with the Bray–Curtis similarity between visits.

Discussion

This study had two major findings, regarding 1) stability of the penile microbiome and 2) longitudinal factors shaping the penile microbiome composition. A substantial proportion of men had stable penile microbiome composition over 1 year as reflected by GBT estimation, individual taxa, and alpha diversity measures. In keeping with the GBT analyses, taxa with the highest mean RA were relatively stable over time for the study sample overall, except for modest decreases in Corynebacterium and Anaerococcus, concomitant with an increasing proportion of men with CT-1 and decreasing CT-2. Through a series of analyses, we identified that circumcision status, HSV-2, condom use, and female sex partner VMB composition are important factors influencing penile microbiome composition and alpha diversity over a 1-year period, which is the longest period of observation of the penile microbiome that we are aware of.

There can be substantial fluctuation in daily (Srinivasan et al., 2010; Song et al., 2020) and weekly (Srinivasan et al., 2010; Mirmonsef et al., 2015) VMB composition, although some studies have shown that community compositions can be quite stable over a period of months (Gajer et al., 2012). VMB composition is temporally influenced by the stage of the menstrual cycle, sexual activity, type of hormonal contraceptive use, diet, exercise, and antibiotic use. In our study, sexual practices (number of sex partners, condom use, and time since the last sex act) were associated with penile microbiome composition. This is intuitive, as these activities represent direct exposure to another microbiome. However, within this cohort of couples who had been sexually active together for at least 6 months (with a median relationship duration of 3 years as reported previously (Mehta et al., 2018)), there was little change in sexual practices (condom use and having multiple sex partners). Nelson et al. also observed a stable penile microbiome at monthly visits in their study of 18 male participants aged 14–17 years, with 61% sexually active (Nelson et al., 2012). Unlike the VMB, the penile microbiome may have fewer factors regulating homeostasis or may be less susceptible to perturbation in general.

Antibiotic use was not consistently associated with penile microbiome composition in our cohort (Supplementary Tables and Figure 5), although only a small number of men were given and/or reported taking antibiotics. Because of the evidence for the sexual exchangeability of the penile microbiome and VMB and their relation to BV, studies are examining the effects of antibiotics on the penile microbiome. Whether oral or topical antibiotic regimens similar to those used by women to treat BV can affect change in the penile microbiome composition persisting for 4 weeks was evaluated by Galiwango et al. among 125 HIV-negative Ugandan men (Galiwango et al., 2019). In a pilot study of 34 heterosexual couples followed up for 12 weeks, Plummer et al. evaluated the effect of standard antimicrobial treatments for BV in women (oral metronidazole and clindamycin cream) on the penile microbiome of men whose partners were treated for BV (Plummer et al., 2021). The presence and abundance of BV-associated bacteria were reduced, and Staphylococcus increased in men’s first-pass urine (representing the urethra), but this was not sustained by 12 weeks posttreatment, despite 81% of women having suppression of BV-associated bacteria over the 12-week period. These findings highlight the need to characterize the male genitourinary microbiome over time and associated factors to understand why effects may not be sustained and to better select and design therapeutic or interventional options related to microbiome manipulation. As evidence supports the sexual exchangeability of the penile microbiome and VMB, antibiotic treatment to modify the penile microbiome may need to be simultaneously administered to female partners. In our study, since 92% of instances of BV in female partners were treated with antibiotics, we were unable to assess the impact of female partner treatment for BV on the penile microbiome. However, this may explain the decline in alpha diversity among men whose partners had BV (Supplementary Figure 3).

Although standard penile CTs have not been established in the same way as vaginal CSTs, our findings share some similarities to dominant penile CTs identified in two cross-sectional studies. Among Ugandan men (Liu et al., 2017), 7 CTs were identified: Corynebacterium and Staphylococcus were more prevalent and abundant in CT-1 to CT-3, while men with the other four CTs had higher prevalence and abundance of Porphyromonas, Prevotella, and Clostridiales Family XI. Men with penile CTs with greater presence and abundance of anaerobes were more likely to have multiple sex partners and female partners with BV. Among South African men (Onywera et al., 2020a), 6 CTs were identified: the most common CT was Corynebacterium dominated, present in 53.4% of men. Other CTs were Lactobacillus dominated (2.5% of men), Prevotella dominated (18.5% of men), Gardnerella dominated (8.8% of men), mixed (9.2% of men), and Chryseobacterium dominated (7.6% of men). As in our study, alpha diversity was significantly lower for the Corynebacterium-dominated CT, and BV-associated bacteria were more predominant in non-Corynebacterium-dominated communities. Men with Corynebacterium-dominated CT were also less likely to harbor high-risk HPV. Differences across studies in penile CT identification and definition may arise from how the penis is sampled; the DNA extraction, sequencing, and annotation processes; and the statistical approach to clustering. For comparability across studies, we need a standardized definition of penile community type, similar to that which has emerged for the VMB (McKinnon et al., 2019; France et al., 2020). Nevertheless, our study and others indicate that a Corynebacterium-dominated, low-diversity penile microbiome community may have beneficial health associations for men and their female partners.

In the female genital tract, a diverse VMB is often associated with non-optimal VMB and poorer health outcomes (McKinnon et al., 2019). Our results contribute toward understanding whether this framework extends to the penile microbiome and whether diversity and the associated penile microbiome compositions influence susceptibility or resilience to poorer health outcomes in men. Penile microbiome alpha diversity was increased when female partners had BV and non-optimal CST, and men with CT-4 and CT-7 (highest Shannon diversity) had the highest prevalence of HSV-2. We observed decreased alpha diversity with circumcision and condom use, indicating that these factors may contribute to the stability of the penile microbiome. Among South African men, alpha diversity was increased among men dually infected with HIV and human papillomavirus (Onywera et al., 2020a). While more research is needed, these results suggest that greater diversity of the penile microbiome may be part of a non-optimal penile microbiome composition. Understanding what constitutes an optimal or non-optimal penile microbiome and identifying intrinsic (e.g., host interactions) and extrinsic (e.g., behavioral and environmental) mediating or modifying factors as has been done for the VMB (Murphy and Mitchell, 2016; McKinnon et al., 2019) will require longitudinal studies defining the penile microbiome and effects on men’s health. Building this knowledge base will have important implications for the health of men and their female partners.

Strengths and Limitations

As reviewed in (Onywera et al., 2020b), most studies of the penile microbiome have been cross-sectional or longitudinal with a small sample size and/or limited duration. In this large, well-characterized longitudinal cohort of community-recruited young Kenyan men and their female partners, we conducted repeated penile microbiota assessments and multifaceted statistical analyses. Although we were able to measure change over time, it was not possible to measure the penile VMB and VMB prior to sexual activity and determine conclusive directionality; even if we had asked couples to refrain from sexual intercourse for a period of time prior to sampling, and then obtain penile and vaginal samples within a short period of time after intercourse, this would not be a true “baseline” because couples were already exposed to each other’s microbiota. Moreover, such a study would have been unethical in this setting, as couples may feel coerced into unwanted sexual practices. This study had limitations inherent to microbiome studies, including measurement error in the assignment of taxa. Unmeasured behavioral, dietary, immunologic, genetic, and other health-related factors that may shape the penile microbiome or could confound the associations observed need to be studied.

Conclusions

The penile microbiome was stable over 1 year for at least half of men; and stable or varying composition was influenced by circumcision, HSV-2, female partner BV status and vaginal CST, and sexual practices. Identification of an optimal penile microbiome through association with health outcomes, a meta-transcriptome study to understand the activity of key taxa, and a comprehensive measure of factors shaping the penile microbiome are necessary to consider interventions to prevent adverse outcomes related to the non-optimal penile microbiome in men and associated conditions in female sex partners.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

This study was approved by the Institutional Review Board of University of Chicago, Illinois USA (2013-0511) and Ethical Review Committee of Maseno University of Kisumu, Kenya (MSU/DRPC/MUERC/00054/13). Written informed consent from the participants’ legal guardian/next of kin was not required to participate in this study in accordance with the national legislation and the institutional requirements.

Author Contributions

SM: obtained funding, study conceptualization and design, implementation of statistical analysis approaches (mixed-effects modeling of alpha diversity, Bray–Curtis similarity), visualization (stacked bar chart, non-metric multidimensional scaling (NMDS), and alluvial diagram), and drafting of the manuscript. DN: design and implementation of statistical analysis approach (specifically hierarchical clustering, mixed-effects modeling of individual taxa, and GBTs), visualization (correlation matrix, coefficient heatmap, taxa, and alpha diversity over time), and critical review and revision of the manuscript. SG: development and oversight of protocols for amplicon sequencing, microbiologic analyses and interpretation, and critical review and revision of the manuscript. WA: development, implementation, and oversight of laboratory protocols in Kenya; acquisition of data; microbiologic analyses and interpretation; and critical review and revision of the manuscript. FO: study oversight and management to ensure integrity to protocols and critical review and revision of the manuscript. DB: conceptualization and design of statistical analysis approach (specifically evaluating overall approach and results) and critical review and revision of the manuscript. RB: study oversight and management to ensure integrity to protocols and critical review and revision of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by grant number R01-AI110369 (PI: Mehta) from the National Institutes of Health, National Institute of Allergy and Infectious Diseases, Division of Microbiology.

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.

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/fcimb.2022.916437/full#supplementary-material

Supplementary Figure 1 | Correlation heatmap showing the Spearman correlations among twelve taxa with highest relative abundance and selected covariates of importance. Spearman rank correlation heatmap of the top 12 penile taxa with greatest relative abundance, and selected covariates. Relative abundance is center log ratio transformed prior to correlation estimation. The correlation heat map represents the direction and magnitude of the Spearman rank correlation between the taxa and covariates. Negative correlations are shaded in blue and positive correlations in red, with deeper intensity representing the magnitude of the correlation. Taxa are represented at the genus level, except where species is noted, for the following: Prevotella timonensis, Sneathia sanguinegens, Gardnerella vaginalis, and Lactobacillus iners.

Supplementary Figure 2 | Mean alpha diversity over time, by circumcision status. The mean alpha diversity measures are shown over time (black lines) with the 95% confidence interval shown in the shaded areas around the mean line. Alpha diversity measures over time are stratified by circumcision status, in which pink represents men who are uncircumcised and blue represents men who are circumcised.

Supplementary Figure 3 | Fractional polynomial fitted plot with 95% confidence intervals for alpha diversity measures over time, stratified by Nugent Score Category. The fitted line plots show the change in (A) Shannon Diversity Index and (B) Richness over time in months (x-axis), stratified by circumcision status and female partner BV status.

References

Bates, D., Mächler, M., Bolker, B., Walker, S. (2015). Fitting Linear Mixed-Effects Models Using Lme4. J. Stat. Software 67, 1–48. doi: 10.18637/jss.v067.i01

CrossRef Full Text | Google Scholar

Brunson, J. C., Read, Q. D. (2020) “Ggalluvial: Alluvial Plots in 'Ggplot2'.” R Package Version 0.12.3. Available at: http://corybrunson.github.io/ggalluvial/.

Google Scholar

Davis, N. M., Proctor, D. M., Holmes, S. P., Relman, D. A., Callahan, B. J. (2018). Simple Statistical Identification and Removal of Contaminant Sequences in Marker-Gene and Metagenomics Data. Microbiome 6, 226. doi: 10.1186/s40168-018-0605-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Desgraupes, B. (2017) Clustering Indices (University Paris Ouest Lab Modal’X). Available at: https://cran.r-project.org/web/packages/clusterCrit/vignettes/clusterCrit.pdf (Accessed November 3, 2021).

Google Scholar

Eren, A. M., Zozaya, M., Taylor, C. M., Dowd, S. E., Martin, D. H., Ferris, M. J. (2011). Exploring the Diversity of Gardnerella Vaginalis in the Genitourinary Tract Microbiota of Monogamous Couples Through Subtle Nucleotide Variation. PloS One 6, e26732. doi: 10.1371/journal.pone.0026732

PubMed Abstract | CrossRef Full Text | Google Scholar

France, M. T., Ma, B., Gajer, P., Brown, S., Humphrys, M. S., Holm, J. B., et al. (2020). VALENCIA: A Nearest Centroid Classification Method for Vaginal Microbial Communities Based on Composition. Microbiome 8, 166. doi: 10.1186/s40168-020-00934-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Gajer, P., Brotman, R. M., Bai, G., Sakamoto, J., Schutte, U. M. E., Zhong, X., et al. (2012). Temporal Dynamics of the Human Vaginal Microbiota. Sci. Transl. Med. 4, 132ra52. doi: 10.1126/scitranslmed.3003605

PubMed Abstract | CrossRef Full Text | Google Scholar

Galiwango, R. M., Bagaya, B., Mpendo, J., Joag, V., Okech, B., Nanvubya, A., et al. (2019). Protocol for a Randomized Clinical Trial Exploring the Effect of Antimicrobial Agents on the Penile Microbiota, Immunology and HIV Susceptibility of Ugandan Men. Trials 20, 443. doi: 10.1186/s13063-019-3545-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hennig, C. Flexible Procedures for Clustering. Version 2.2-9. Available at: https://cran.r-project.org/web/packages/fpc/fpc.pdf.

Google Scholar

Herve, M. Package ‘Rvaidememoire’. Version 0.9-81-2. Available at: https://cran.r-project.org/web/packages/RVAideMemoire/index.html.

Google Scholar

Holm, J. B., Humphrys, M. S., Robinson, C. K., Settles, M. L., Ott, S., Fu, L., et al. (2019). Ultrahigh-Throughput Multiplexing and Sequencing of >500-Base-Pair Amplicon Regions on the Illumina HiSeq 2500 Platform. mSystems 4, e00029–e00019. doi: 10.1128/mSystems.00029-19

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, B. L., Nagin, D. S. (2013). A Note on a Stata Plugin for Estimating Group-Based Trajectory Models. Sociological Methods Res. 42, 608–613. doi: 10.1177/0049124113503141

CrossRef Full Text | Google Scholar

Koren, O., Knights, D., Gonzalez, A., Waldron, L., Segata, N., Knight, R., et al. (2013). A Guide to Enterotypes Across the Human Body: Meta-Analysis of Microbial Community Structures in Human Microbiome Dataset. PloS Comput. Biol. 9, 59–65. doi: 10.1371/journal.pcbi.1002863

CrossRef Full Text | Google Scholar

Liu, C. M., Hungate, B. A., Tobian, A. A. R., Serwadda, D., Ravel, J., Lester, R., et al. (2013). Male Circumcision Significantly Reduces Prevalence and Load of Genital Anaerobic Bacteria. mBio 4, e00076. doi: 10.1128/mBio.00076-13

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, C. M., Prodger, J. L., Tobian, A. A. R., Abraham, A. G., Kigozi, G., Hungate, B. A., et al. (2017). Penile Anaerobic Dysbiosis as a Risk Factor for HIV Infection. mBio 8, e00996–e00917. doi: 10.1128/mBio.00996-17

PubMed Abstract | CrossRef Full Text | Google Scholar

Maechler, M., Rousseeuw, P., Struyf, A., Hubert, M., Hornik, K. (2022) R Package Version 2.1.3. Available at: https://CRAN.R-project.org/package=cluster.

Google Scholar

McKinnon, L. R., Achilles, S. L., Bradshaw, C. S., Burgener, A., Crucitti, T., Fredricks, D. N., et al. (2019). The Evolving Facets of Bacterial Vaginosis: Implications for HIV Transmission. AIDS Res. Hum. Retroviruses 35, 219–228. doi: 10.1089/AID.2018.0304

PubMed Abstract | CrossRef Full Text | Google Scholar

Mehta, S. D., Nandi, D., Agingu, W., Green, S. J., Bhaumik, D. K., Bailey, R. C., et al. (2020). Vaginal and Penile Microbiome Associations With HSV-2 in Women and Their Male Sex Partners. J. Infect. Dis. 21, jiaa529. doi: 10.1093/infdis/jiaa529

CrossRef Full Text | Google Scholar

Mehta, S. D., Nordgren, R. K., Agingu, W., Otieno, F., Ochieng, W., Odhiambo, F., et al. (2018). Sexual Quality of Life and Association With HIV and STIs Among a Cohort of Heterosexual Couples in Kenya. J. Sexual Med. 15, 1446–1455. doi: 10.1016/j.jsxm.2018.08.007

CrossRef Full Text | Google Scholar

Mehta, S. D., Zhao, D., Green, S. J., Agingu, W., Otieno, F., Bhaumik, R., et al. (2020). The Microbiome Composition of a Man’s Penis Predicts Incident Bacterial Vaginosis in His Female Sex Partner With High Accuracy. Front. Cell Infect. Microbiol. 10, 433. doi: 10.3389/fcimb.2020.00433

PubMed Abstract | CrossRef Full Text | Google Scholar

Mirmonsef, P., Modur, S., Burgad, D., Gilbert, D., Golub, E. T., French, A. L., et al. (2015). Exploratory Comparison of Vaginal Glycogen and Lactobacillus Levels in Premenopausal and Postmenopausal Women. Menopause 22, 702–709. doi: 10.1097/GME.0000000000000397

PubMed Abstract | CrossRef Full Text | Google Scholar

Murphy, K., Mitchell, C. M. (2016). The Interplay of Host Immunity, Environment and the Risk of Bacterial Vaginosis and Associated Reproductive Health Outcomes. J. Infect. Dis. 214, S29–S35. doi: 10.1093/infdis/jiw140

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagin, D. S., Odgers, C. L. (2010). Group-Based Trajectory Modeling in Clinical Research. Ann. Rev. Clin. Psychol. 6, 109–138. doi: 10.1146/annurev.clinpsy.121208.131413

CrossRef Full Text | Google Scholar

Naqib, A., Poggi, S., Wang, W., Hyde, M., Kunstman, K., Green, S. J. (2018). Making and Sequencing Heavily Multiplexed, High-Throughput 16s Ribosomal RNA Gene Amplicon Libraries Using a Flexible, Two-Stage PCR Protocol. Methods Mol. Biol. 1783, 149–169. doi: 10.1007/978-1-4939-7834-2_7

PubMed Abstract | CrossRef Full Text | Google Scholar

National, A. I. D. S., Control Programme, S. T. I. (2015) Ministry of Health, Kenya. Guidelines for HIV Testing Services in Kenya (Nairobi: NASCOP). Available at: https://www.fast-trackcities.org/sites/default/files/Kenya%20HIV%20Testing%20Services%20Guidelines%20%282015%29.pdf (Accessed March 1, 2022).

Google Scholar

Nelson, D. E., Dong, Q., van der Pol, B., Toh, E., Fan, B., Katz, B. P., et al. (2012). Bacterial Communities of the Coronal Sulcus and Distal Urethral of Adolescent Males. PloS One 7, e36298. doi: 10.1371/journal.pone.0036298

PubMed Abstract | CrossRef Full Text | Google Scholar

Nugent, R. P., Krohn, M. A., 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

Oksanen, J., Simpson, G. L., Blanchet, F. G., Roeland, K., Legendre, P., Minchin, P. R., et al. (2022) Vegan: Community Ecology Package. R Package Version 2.6-2. Available at: https://CRAN.R-project.org/package=vegan.

Google Scholar

Onywera, H., Williamson, A.-L., Cozzuto, L., Bonnin, S., Mbulawa, Z. Z. A., Coetzee, D., et al. (2020a). The Penile Microbiota of Black South African Men: Relationship With Human Papillomavirus and HIV Infection. BMC Microbiol. 20, 78. doi: 10.1186/s12866-020-01759-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Onywera, H., Williamson, A.-L., Cozzuto, L., Bonnin, S., Mbulawa, Z. Z. A., Coetzee, D., et al. (2020b). The Penile Microbiota of Black South African Men: Relationship With Human Papillomavirus and HIV Infection. BMC Microbiol. 20, 78. doi: 10.1186/s12866-020-01759-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Palarea-Albaladejo, J., Martin-Fernandez, J. A. (2015). Zcompositions – R Package for Multivariate Imputation of Left-Censored Data Under a Compositional Approach. Chemometr Intell. Lab. Syst. 143, 85–96. doi: 10.1016/j.chemolab.2015.02.019

CrossRef Full Text | Google Scholar

Plummer, E. L., Vodstrcil, L. A., Doyle, M., Danielewski, J. A., Murray, G. L., Fehler, G., et al. (2021). A Prospective, Open-Label Pilot Study of Concurrent Male Partner Treatment for Bacterial Vaginosis. mBio 12, e0232321. doi: 10.1128/mBio.02323-21

PubMed Abstract | CrossRef Full Text | Google Scholar

Posada, D., Buckley, T. R. (2004). Model Selection and Model Averaging in Phylogenetics: Advantages of Akaike Information Criterion and Bayesian Approaches Over Likelihood Ratio Tests. Syst Biol. 53, 793–808. doi: 10.1080/10635150490522304

PubMed Abstract | CrossRef Full Text | Google Scholar

Price, L. B., Liu, C. M., Johnson, K. E., Aziz, M., Lauk, M. K., Bowers, J., et al. (2010). The Effects of Circumcision on the Penis Microbiome. PloS One 5, e8422. doi: 10.1371/journal.pone.0008422

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, S. D., Acharya, K. D., Zhu, J. E., Deveney, C. M., Walther-Antonio, M. R. S., Tetel, M. J., et al. (2020). Daily Vaginal Microbiota Fluctuations Associated With Natural Hormonal Cycle, Contraceptives, Diet, and Exercise. mSphere 5, e00593–e00520. doi: 10.1128/mSphere.00593-20

PubMed Abstract | CrossRef Full Text | Google Scholar

Srinivasan, S., Liu, C., Mitchell, C. M., Fiedle, T. L., Thomas, K. K., Agnew, K. J., et al. (2010). Temporal Variability of Human Vaginal Bacteria and Relationship With Bacterial Vaginosis. PloS One 5, e10197. doi: 10.1371/journal.pone.0010197

PubMed Abstract | CrossRef Full Text | Google Scholar

Wickham, H. (2016) (New York).

Google Scholar

Zhang, J., Kobert, K., Flouri, T., Stamatakis, A. (2014). PEAR: A Fast and Accurate Illumina Paired-End reAd Merger. Bioinformatics 30, 614–620. doi: 10.1093/bioinformatics/btt593

PubMed Abstract | CrossRef Full Text | Google Scholar

Zozaya, M., Ferris, M. J., Siren, J. D., Lillis, R., Myers, L., Nsuami, M. J., et al. (2016). Bacterial Communities in Penile Skin, Male Urethra, and Vaginas of Heterosexual Couples With and Without Bacterial Vaginosis. Microbiome 4, 16. doi: 10.1186/s40168-016-0161-6

PubMed Abstract | CrossRef Full Text | Google Scholar