Longitudinal cervicovaginal microbiome and virome alterations during ART and discordant shedding in women living with HIV

doi:10.21203/rs.3.rs-4078561/v1

Download PDF

Article

Longitudinal cervicovaginal microbiome and virome alterations during ART and discordant shedding in women living with HIV

https://doi.org/10.21203/rs.3.rs-4078561/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Despite successful suppression of plasma HIV replication by antiretroviral therapy (ART), some women living with HIV (WLHIV) can still experience genital HIV shedding (discordant shedding). Female genital tract (FGT) microbiome and virome dynamics during long-term ART in WLHIV are poorly understood but might contribute to discordant HIV shedding, as the microbiome and virome are known to influence FGT health. To understand FGT microbial communities over time during ART usage and discordant shedding, we characterized the microbiome and virome in 125 cervicovaginal specimens collected over two years in 31 WLHIV in Lima, Peru. Intrapersonal bacterial microbiome variation was higher in HIV shedders compared to non-shedders. Cervicovaginal virome composition changed over time, particularly in non-shedders. Specifically, anellovirus relative abundance was inversely associated with ART duration and CD4 counts. Our results suggest that discordant HIV shedding is associated with FGT microbiome instability, and immune recovery during ART influences FGT virome composition.

Biological sciences/Microbiology/Microbial communities/Microbiome

Health sciences/Diseases/Infectious diseases/HIV infections

With the development of effective antiretroviral therapies (ART), life expectancy for people living with HIV has increased significantly and HIV infection can now be treated as a chronic disease^1,2. Of the estimated 20 million women aged 15 and older living with HIV worldwide in 2022, an estimated 82% were taking ART³. While ART is effective in suppressing HIV replication, some women living with HIV (WLHIV) with undetectable plasma viral loads can still experience genital HIV shedding, described as discordant shedding^4-10. The microbiome has been implicated in female genital tract (FGT) health and disease¹¹. However, the relationship between discordant shedding and the FGT microbiome and virome during long-term ART in WLHIV remains poorly understood.

The vaginal microbiomes of reproductive-aged women typically consist mainly of Lactobacillus or diverse anaerobes including Gardnerella, Prevotella, Sneathia, Fannyhessea, Megasphaera, Dialister, Mobiluncus, Peptonophilus, Shuttleworthia, Parvimonas and other genera^11-15. Lactobacilli contribute to vaginal health by producing lactic acid, which lowers vaginal pH and creates a less favorable environment for pathogenic bacteria, and by competing with pathogens for space¹¹. Cervicovaginal microbiome alterations have been associated with inflammation¹⁶, preterm birth¹⁷, changes to the mucosal epithelial barrier¹⁸, cervical cancer^19,20, and sexually transmitted infections^14,20-22. Although higher Lactobacillus relative abundance has been associated with vaginal health, many asymptomatic women have diverse vaginal microbiomes dominated by bacteria other than lactobacilli^12-15. Longitudinal studies have demonstrated that the vaginal microbiome is stable over time in some women, while in others vaginal microbiome composition fluctuates more frequently^13,16.

The FGT microbiome is particularly important in HIV infection. Specific microbiome community states and bacterial species have been associated with increased HIV prevalence²² and acquisition¹⁴. In WLHIV, cervicovaginal microbiota have been associated with genital HIV shedding^5,22,23. While pre-exposure prophylaxis (PrEP) has not been shown to significantly influence the vaginal microbiome in women without HIV^24-26, multiple lines of evidence suggest that vaginal bacteria influence metabolism and concentrations of antiretroviral drugs in the FGT^24,27-32. Whether long-term ART influences the vaginal microbiome in WLHIV is poorly understood. One study found no microbiome changes after six months of ART³³, another reported higher Lactobacillus relative abundance in women with viral suppression after a year of ART³⁴, and a third reported an association between bacterial vaginosis and ART with protease inhibitors compared to regimens without protease inhibitors³⁵. Further study is needed to clarify the relationship between ART and the FGT microbiome over years or decades. The relationship between the microbiome and discordant HIV shedding is also unclear. Discordant HIV shedding is incompletely understood, but has been associated with genital inflammation³⁶, coinfections^36,37, trauma³⁸ and ulcers¹⁰, and has also been attributed to expression of integrated provirus from proliferating cells⁴. The importance of the microbiome in FGT health led us to consider the microbiome as a potential factor promoting discordant HIV shedding.

Compared to the bacterial microbiome, less is known about the role of the virome in FGT health. Viruses found in the FGT include papillomaviruses, anelloviruses, herpesviruses and bacteriophages^39-47. Vaginal virome community structure has been linked with the bacterial microbiome^40-42. Vaginal virome alterations have been associated with inflammation⁴¹, cervical lesions^43,44, and preterm birth⁴⁵. Human papillomaviruses are the major cause of cervical cancer, the fourth most-diagnosed cancer in women worldwide in 2020⁴⁸. Anelloviruses, a recently discovered, prevalent viral family⁴⁹, have been associated with FGT inflammation in women without HIV⁴¹ and CD4 T cell levels in WLHIV⁴⁴. Little is known about the relationship between the FGT virome and long-term ART or discordant HIV shedding.

To investigate the relationship between the cervicovaginal microbiome and virome and discordant shedding over time in WLHIV taking ART, we analyzed the microbiome and virome in 125 cervicovaginal specimens collected over a two-year period from 31 WLHIV in Lima, Peru, nine of whom experienced discordant HIV shedding. Intrapersonal bacterial microbiome stability differed significantly between HIV shedders and non-shedders. Significant virome changes were associated with ART duration, specifically decreases in anellovirus abundance, especially among non-shedders. Our results suggest that cervicovaginal microbiome instability is associated with discordant HIV shedding, while immune recovery during ART may influence the cervicovaginal virome.

Shedders’ bacterial microbiomes are less stable over time than non-shedders’ microbiomes.

To investigate the roles of both the microbiome and virome in discordant HIV shedding, we analyzed 125 cervicovaginal specimens collected longitudinally from 31 WLHIV in Lima, Peru (Table 1, Fig. 1A, Supplementary Table 1). Specimens were collected quarterly from participants over a period of two years as part of a larger study characterizing discordant shedding in people living with HIV⁴. In this study, we analyzed an average of 4 specimens per participant. Most specimens (n = 120) were cervicovaginal lavages (CVL). When CVL samples were not available, vaginal swabs (n = 2), cytobrushes (n = 2), or TearFlo paper strips (n = 1) were used to sample the cervicovaginal environment. Nine women experienced discordant HIV shedding during the study period (shedders) and 22 did not (non-shedders). Non-shedders were older than shedders at study enrollment (non-shedder median age = 34 years, shedder median age = 28 years; Mann-Whitney test, p = 0.047) (Extended Data Fig. 1A). Shedders and non-shedders had similar CD4 counts at enrollment and similar CD4 recoveries when values from enrollment and after 2 years of ART were compared (Mann-Whitney test; p = 0.94 at enrollment, p = 0.94 at final visit, p = 0.94 for change from enrollment to final visit) (Extended Data Fig. 1B). Using metagenomic next-generation sequencing (NGS) after physically enriching for bacteria, we obtained a median of 15.6x10⁶ (IQR 13.7x10⁶ – 17.8x10⁶) raw read pairs and 1.95x10⁶ (IQR 8.7x10⁵ – 4.7x10⁶) quality-filtered, deduplicated individual reads per sample. Most samples clustered into groups dominated by Lactobacillus crispatus (n = 29), Lactobacillus iners (n = 21), Gardnerella vaginalis (n = 48), or a mixed community (n = 21), with a few samples dominated by Lactobacillus jensenii (n = 2) or Escherichia coli (n = 4) (Fig. 1B). There were no significant differences in bacterial alpha diversity, beta diversity, or cluster assignment between different sample types (CVL, vaginal swabs, cytobrushes or TearFlo strips) (Extended Data Fig. 1C-1F, Supplementary Table 2).

Table 1

Cohort characteristics
	Shedders	Non-shedders	Statistical significance
Number of participants	9	22	N/A
Number of specimens analyzed per participant: median (interquartile range)	4 (4–4)	4 (4–4)	p = 0.52
Age at enrollment: median (interquartile range)	28 (24–30)	34 (30–39)	p = 0.047
CD4 at enrollment: median (interquartile range)	132 (89–223)	141 (51.25–205.5)	p = 0.94
CD4 after 2 years of ART: median (interquartile range)	454 (307.5–510)	356.5 (295.5–415.5)	p = 0.94
CD4 change from initial to final visit: median (interquartile range)	209.5 (140.25–304.5)	208 (152–329.75)	p = 0.94
Statistical significance was assessed by the Mann-Whitney test, with the Benjamini-Hochberg method used to correct for multiple comparisons where appropriate. Number of specimens, age, and CD4 counts are expressed as median and interquartile range.

Intrapersonal microbiome variation (beta diversity) was lower than interpersonal variation, suggesting intraindividual stability over time (Fig. 2A; Bray-Curtis dissimilarity; Mann-Whitney test, p = 1.9x10^− 19 for non-shedders, p = 0.0019 for shedders). Shedders had higher intrapersonal microbiome beta diversity than non-shedders (Fig. 2A; Bray-Curtis dissimilarity, Mann-Whitney test, p = 0.0019). We tested the hypothesis that perturbation was most pronounced during discordant shedding by comparing diversity between discordant and non-discordant timepoints within individuals (Fig. 2B). Beta diversity was significantly higher between samples when one timepoint was discordant and the other was non-discordant, compared to pairs of non-discordant timepoints (Bray-Curtis dissimilarity; Mann-Whitney test, p = 0.014), suggesting that higher intrapersonal variation in shedders reflects dysbiosis during discordant shedding rather than prolonged fluctuations. Principal coordinates analysis (PCoA) using Bray-Curtis dissimilarity showed no association between microbiome composition and ART duration (PERMANOVA, p = 0.33) or shedder/non-shedder status (PERMANOVA, p = 0.34) (Fig. 2C), although discriminant analysis with MaAsLin2⁵⁰ showed higher Mycoplasma hominis relative abundance in shedders compared to non-shedders (Fig. 2D, q = 0.048). PCoA of shedder samples showed no difference between discordant and non-discordant timepoints (Fig. 2E; PERMANOVA, p = 0.85). Microbiome community clusters were not associated with shedder/non-shedder status, ART duration, or discordant shedding timepoints (Fig. 1B, Supplementary Table 2). These results suggest that no specific microbiome composition is associated with shedder/non-shedder status, ART duration, or discordant shedding events, implying that higher variation during discordant shedding reflects individual-specific, not shared, dysbioses. To determine whether this could be explained by common functional characteristics encoded by bacterial genomes that might not be reflected in taxonomic profiles, we performed functional profiling with HUMAnN3⁵¹ (Extended Data Fig. 2A). We observed similar results for bacterial metabolic pathway and gene family beta diversity as for taxonomic beta diversity (Extended Data Fig. 2B-2I). 172 metabolic pathways were associated with microbiome community clusters but not with shedder/non-shedder status, discordant shedding events, or ART duration (MaAsLin2, Supplementary Table 3). Results were similar for gene families (Supplementary Table 4). This suggests that functional profiles mirror taxonomic profiles, with similar beta diversity dynamics during discordant HIV shedding.

Non-shedders’ viromes are less stable than shedders’ viromes, with anellovirus and papillomavirus relative abundance changing over time.

Given the temporal dynamics of the bacterial microbiome, we examined the virome for similar patterns. After enrichment for virus-like particles (VLPs), we visualized VLPs in specimens using epifluorescence microscopy to corroborate VLP enrichment (Extended Data Fig. 3A). VLP concentrations ranged from 2.5x10⁵ – 2.2x10⁸ particles per mL of VLP-enriched supernatant, with a median of 1.1x10⁷ particles per mL (Extended Data Fig. 3B). Using metagenomic NGS of nucleic acids extracted from VLP-enriched supernatants and amplified by multiple displacement amplification (MDA), we obtained a median of 1.6x10⁷ (IQR 4.5x10⁶ – 2.2x10⁷) raw read pairs and 2.4x10⁶ (IQR 1.1x10⁶ – 3.5x10⁶) quality-filtered, deduplicated individual reads per sample. 12,574 unique contigs > = 500 nt were assembled from the samples and sequencing controls, and 1,596 were identified as viral. After mapping quality-filtered reads to the viral contigs, 116 samples had > 1,000 viral reads and were retained for further analysis. Samples clustered into groups dominated by Anelloviridae (n = 46), Papillomaviridae (n = 54), or bacteriophages including Microviridae and Caudoviricetes (n = 16) (Fig. 3A). Virome alpha diversity was inversely associated with ART duration (linear mixed modeling (LME), p = 7.7x10^− 10 for viral contig richness, p = 8.3x10^− 8 for Shannon index), with no significant difference between shedders and non-shedders (LME, p = 0.19 for richness, p = 0.90 for Shannon index) (Fig. 3B). Because MDA only amplifies DNA, we performed sequence-independent amplification (SIA) on the viral nucleic acids to amplify RNA and DNA. We found plant-infecting RNA viruses in a few samples, but no RNA viruses known to infect humans or bacteria, and most samples had no RNA viruses identified (Extended Data Fig. 4A). Therefore, we focused our analyses on the DNA virome using the MDA data.

Like bacterial microbiome variation, intrapersonal virome variation (beta diversity) was significantly lower than interpersonal variation (Fig. 3C, Bray-Curtis dissimilarity (weighted); Mann-Whitney, p = 1.4x10^− 36 for non-shedders, p = 1.2x10^− 12 for shedders; Fig. 3D, Sorensen dissimilarity (unweighted); Mann-Whitney, p = 2.2x10^− 13 for non-shedders, p = 9.3x10^− 5 for shedders). However, weighted intrapersonal variation was not significantly different between shedders and non-shedders (Fig. 3C, Bray-Curtis dissimilarity; Mann-Whitney, p = 0.71). Further, unweighted intrapersonal variation was significantly higher in non-shedders than in shedders (Sorensen dissimilarity; Mann-Whitney test, p = 0.0020) (Fig. 3D). In weighted (Bray-Curtis) PCoA, virome composition was associated with ART duration (PERMANOVA, p = 0.001) and shedder/non-shedder status (PERMANOVA, p = 0.016) (Fig. 3E). Because ART duration trended toward interaction with shedder/non-shedder status (PERMANOVA, p = 0.069), we carried out PCoA on shedders and non-shedders separately (Extended Data Fig. 5A). Virome composition was associated with ART duration in non-shedders (PERMANOVA, p = 0.001), and trended toward association with ART duration in shedders but was not significant (PERMANOVA, p = 0.072). In unweighted (Sorensen) PCoA of combined shedders and non-shedders, virus presence-absence was associated with ART duration (PERMANOVA, p = 0.001) but not shedder/non-shedder status (PERMANOVA, p = 0.63), with no interaction between ART duration and shedder/non-shedder status (PERMANOVA, p = 0.15) (Extended Data Fig. 5B). These results suggest that virome composition changes over time in WLHIV taking ART, with more pronounced changes among non-shedders.

We used differential abundance analysis (MaAsLin2⁵⁰) to identify viruses responsible for changes in virome beta diversity over time. Anelloviridae relative abundance was inversely associated with ART duration (q = 8.8x10^− 5), while Papillomaviridae relative abundance increased over time (q = 0.037) (Fig. 4F). ART duration trended toward interaction with shedder/non-shedder status for several viral families, but these trends were not significant (Supplementary Table 5). When non-shedders were analyzed separately, Anelloviridae relative abundance decreased over time after ART initiation and Papillomaviridae relative abundance increased (Extended Data Fig. 5C). No taxa changed significantly over time in shedders. These results suggest that changes in virome composition after ART initiation are driven by Anelloviridae and Papillomaviridae, especially among non-shedders.

Anellovirus communities change over time, with decreased alphatorquevirus abundance and increased betatorquevirus proportions.

Because the Anelloviridae family contributed significantly to the temporal virome signature (Fig. 3), we performed phylogenetic analysis of the anellovirus contig ORF1 gene sequences to better define anellovirus community dynamics. We obtained predicted ORF1 sequences from 751/799 anellovirus contigs. 21 contigs had more than one predicted ORF1 sequence, for a total of 774 sequences. At the genus level, most anelloviruses were classified as alphatorqueviruses (n = 339), betatorqueviruses (n = 339), or gammatorqueviruses (n = 85), with a few hetorqueviruses (n = 4) and sequences falling outside known anellovirus genera (n = 7) (Fig. 4A). Compared to all viral genera, relative abundance of alphatorqueviruses and incomplete anellovirus genomes (i.e., without ORF1) decreased over time after ART initiation (MaAsLin2, q = 0.0019 and q = 0.0019) (Extended Data Fig. 6A). There were no significant changes in the relative abundance of other anellovirus genera. This suggests that decreased Anelloviridae relative abundance reflects decreased Alphatorquevirus relative abundance.

We used the anellovirus phylogeny to calculate UniFrac distance, a beta diversity metric that accounts for phylogenetic distance. In shedders, intrapersonal anellovirus beta diversity was significantly lower than interpersonal diversity (weighted UniFrac, Mann-Whitney, p = 2.7x10^− 4), while for non-shedders there was no significant difference between intrapersonal and interpersonal anellovirus beta diversity (weighted UniFrac, Mann-Whitney, p = 0.33) (Fig. 4B). Intrapersonal variation was also significantly lower in shedders than non-shedders (Fig. 4B; Mann-Whitney, p = 3.5x10^− 5). Results were similar for unweighted UniFrac distance (Extended Data Fig. 6B). These results suggest that anellovirus communities are less stable over time in non-shedders compared to shedders. In PCoA using weighted anellovirus UniFrac distance, anellovirus community compositions changed significantly by ART duration (PERMANOVA, p = 0.005), with no significant difference between shedders and non-shedders (PERMANOVA, p = 0.81) (Fig. 4C). Because ART duration trended toward interaction with shedder/non-shedder status (PERMANOVA, p = 0.078), we carried out separate PCoA analyses on non-shedder and shedder samples. In these separate analyses, ART duration was associated with anellovirus community structure in non-shedders (PERMANOVA, p = 0.006), but not in shedders (PERMANOVA, p = 0.12) (Extended Data Fig. 6C). In unweighted anellovirus UniFrac PCoA of combined shedders and non-shedders, anellovirus presence-absence was associated with ART duration (PERMANOVA, p = 0.001) but not with shedder/non-shedder status (PERMANOVA, p = 0.59), with no interaction between ART duration and shedder/non-shedder status (PERMANOVA, p = 0.31) (Extended Data Fig. 6D). We performed differential abundance analyses (MaAsLin2) to identify specific anelloviruses responsible for changing anellovirus diversity over time. When shedders and non-shedders were analyzed together, Betatorquevirus relative abundance increased and relative abundance of partial anellovirus genomes decreased relative to other anellovirus genera (Extended Data Fig. 6E; MaAsLin2, q = 0.036 and q = 0.036). When non-shedders were considered separately, Betatorquevirus relative abundance increased (MaAsLin2, q = 0.028) (Extended Fig. 6F) and Alphatorquevirus relative abundance tended to decrease but did not reach significance (MaAsLin2, q = 0.060). No anellovirus genera changed significantly over time in shedders. Together with the whole-virome analyses, where Alphatorquevirus relative abundance decreased over time in the virome as a whole, these results suggest that betatorqueviruses consequently make up a larger proportion of the remaining anellovirus community, and these changes are more pronounced among non-shedders.

We corroborated the NGS data with quantitative PCR (qPCR) specific to the untranslated region of the alphatorquevirus genome (Extended Data Fig. 7Α). Alphatorquevirus genome copy numbers were correlated with Alphatorquevirus NGS relative abundance (linear regression, p = 0.041) and trended toward association with Anelloviridae relative abundance (linear regression, p = 0.085), but not Betatorquevirus or Gammatorquevirus relative abundance (linear regression, p = 0.96 and p = 0.84, respectively) (Extended Data Fig. 7Β). Alphatorquevirus copy numbers tended toward inverse association with ART duration, but this was not significant (LME, p = 0.090) (Extended Data Fig. 7C), and there was no significant difference between shedders and non-shedders (LME, p = 0.51).

Because anellovirus abundance has previously been associated with immune status^44,52–59, we considered whether increasing CD4 counts could explain decreasing anellovirus relative abundance. CD4 counts were inversely associated with Anelloviridae relative abundance (LME, p = 1.1x10^− 4) (Fig. 4D), Alphatorquevirus relative abundance (Extended Data Fig. 7D, LME, p = 1.7x10^− 4) and alphatorquevirus qPCR copy numbers (Fig. 4E; LME, p = 0.04). There was no association between CD4 counts and Betatorquevirus or Gammatorquevirus relative abundance (LME, p = 0.77 and p = 0.36, respectively) (Extended Data Fig. 7D). When we analyzed non-shedders separately, Anelloviridae NGS relative abundance, Alphatorquevirus NGS relative abundance and alphatorquevirus qPCR copy numbers were inversely associated with CD4 counts (LME; p = 2.1x10^− 4, p = 1.3x10^− 4, and p = 0.04, respectively), but there were no significant associations in shedders (Extended Data Fig. 7E). These results suggest a relationship between immune status and anellovirus abundance, specifically alphatorqueviruses, in the FGT.

Papillomavirus populations are stable over time in shedders and non-shedders.

Because the Papillomaviridae family also contributed to the temporal virome signature (Fig. 3), we conducted phylogenetic analysis of the concatenated papillomavirus contig E1, E2, L2, and L1 genes to characterize papillomavirus diversity over time. We obtained E1-E2-L2-L1 sequences from 169/190 papillomavirus contigs. Most were alphapapillomaviruses (n = 128), with some betapapillomaviruses (n = 7) and gammapapillomaviruses (n = 34) (Fig. 5A). Alphapapillomavirus relative abundance increased over time (Extended Data Fig. 8A, MaAsLin2, q = 0.044), while other papillomavirus genera did not change significantly.

Unlike anellovirus beta diversity, both non-shedders and shedders had lower intrapersonal papillomavirus beta diversity compared to interpersonal diversity, with no significant difference between shedders and non-shedders (Fig. 5B; weighted papillomavirus UniFrac distance, Mann-Whitney, p = 8.8x10^− 9, p = 0.0083, and p = 0.95, respectively). Unweighted papillomavirus beta diversity results were similar (Extended Data Fig. 8B). In PCoA using papillomavirus UniFrac distance, there were no significant differences by ART duration or shedder/non-shedder status (Fig. 5C, Extended Data Fig. 8C; PERMANOVA; weighted, p = 0.33 and p = 0.28, respectively; unweighted, p = 0.18 and p = 0.10, respectively). These results suggest that papillomavirus diversity is stable over time in shedders and non-shedders.

We identified intraindividual and longitudinal patterns in the cervicovaginal microbiome and virome of WLHIV with and without discordant HIV shedding. Specifically, the bacterial microbiome was less stable over time in shedders. Virome composition changed significantly in association with ART duration, with more pronounced changes among non-shedders.

Higher intraindividual bacterial beta diversity in shedders suggests more frequent changes in microbiome composition compared to non-shedders. This was driven by shedding timepoints, suggesting that microbiome alterations may be associated with shedding events. By contrast, lower intraindividual diversity in non-shedders suggests stability over time. No specific microbiome composition was associated with shedders or shedding events, suggesting that stability over time is more important than specific taxa in the context of discordant shedding. This is interesting because previous studies have often focused on specific microbiome community states or bacterial taxa as markers of FGT health¹¹. Previous analysis of this cohort linked discordant shedding with proliferation of infected cells, rather than HIV replication⁴. Other studies have identified associations between cervicovaginal microbiota and T cell populations in the FGT¹⁴. Thus, the microbiome may contribute to discordant shedding by influencing the number of infected cells in the FGT, that may proliferate in response to disruption of the FGT. Alternatively, the microbiome may be influenced by other factors (trauma through intercourse, inflammation, etc.) that also contribute to discordant shedding.

We identified significant longitudinal changes in the cervicovaginal virome. Specifically, Anelloviridae and Alphatorquevirus relative abundance decreased and overall anellovirus community diversity changed during two years of ART. Anellovirus relative abundance and alphatorquevirus absolute abundance measured by qPCR were inversely associated with plasma CD4 counts. Anelloviruses are prevalent in healthy adult humans and are not known to be pathogenic⁴⁹. Interestingly, numerous studies have associated anellovirus loads with immune suppression, including in people living with HIV^44,52–58 and in organ transplant recipients^49,59. Thus, the anellovirus dynamics that we observed may reflect participants’ immune recovery during two years of ART. The fact that the associations were more pronounced among non-shedders may be a consequence of sampling size, or may suggest underlying immune dysregulation in shedders.

While Papillomaviridae relative abundance increased over time, papillomavirus diversity did not change significantly. Because relative abundance data is compositional, papillomavirus relative abundance may have been inflated when anellovirus relative abundance decreased, with no underlying change in papillomavirus viral load. Although Phi29 polymerase may preferentially amplify small circular single-stranded DNA⁶⁰, our alphatorquevirus qPCR absolute copy numbers were consistent with NGS relative abundance, suggesting that our results represent an actual biological signature. Besides plant virus reads identified in a few samples, we did not find any RNA viruses by sequence-independent amplification (SIA). This could mean that no RNA viruses were present in the specimens, or the RNA could have degraded during storage.

Our cohort consisted of ART-naïve women from a clinical setting who were initiating government-sponsored ART in Lima, Peru. Because all the study participants received ART, we cannot directly compare our results with data from WLHIV not receiving ART. While this comparison would be biologically informative, such a study would be incompatible with current treatment standards. To our knowledge, only a few studies have analyzed the cervicovaginal microbiome or virome in women from Peru^61–63. Therefore, our data contributes valuable information about the FGT microbiome and virome in women from this region. Prevalence of vaginal microbiome community states varies between cohorts, with some studies reporting differences between women from different geographic regions or different racial and ethnic backgrounds^{11–15,19,64–66}. Therefore, studying diverse groups of WLHIV worldwide would indicate whether our results are generalizable to other populations.

This longitudinal analysis of the cervicovaginal virome and microbiome in WLHIV identifies intrapersonal and temporal patterns associated with discordant HIV shedding and ART duration. These findings support further investigation of the role of bacterial communities in discordant HIV shedding in women with suppressed plasma viral loads, as well as the role of anelloviruses during immune suppression and recovery.

Specimens: All study participants provided informed consent. This study was approved by the Institutional Review Boards of the Hospital Dos de Mayo (Lima, Peru), Seattle Children's Hospital (Seattle, Washington, USA), and Arizona State University (Tempe, Arizona, USA). Cervicovaginal specimens were collected as part of a 18-24-month study evaluating low-level viremias (aka “blips”) and discordant genital shedding in ART naïve men and women from Lima, Peru^4,67. Specimens were collected every three months over a two-year period for each participant, and plasma CD4 counts were measured every six months. Cervicovaginal lavages were collected by washing the uterine cervix with sterile 1X phosphate buffered saline (PBS) and collecting the 1X PBS from the vaginal fornix as described⁶⁸. Plasma and CVL HIV RNA were quantified by qRT-PCR as described⁶⁹. Discordant shedding was assessed in participants with ART suppression, defined as median plasma HIV RNA below the limit-of-quantification (1.48 log₁₀ copies(c)/mL), without virologic failure (HIV RNA > 3.0 log₁₀c/mL), who had HIV RNA detected in the CVL (1.48 log₁₀ copies(c)/mL) for the women. Nine women experienced discordant HIV shedding (shedders) and 22 women did not experience discordant shedding (non-shedders). To select timepoints for microbiome and virome analysis, non-shedders were matched with shedders based on average CD4 counts over the duration of ART (2–4 non-shedders for each shedder). We analyzed the microbiome and virome at enrollment (month 3 for one shedder who did not have an enrollment specimen available); for shedders, discordant shedding timepoints and immediate pre- or post- shedding timepoints, as available; matched visits for non-shedders; and at 24 months (month 21 for two shedders who did not have month 24 specimens available) (Fig. 1A). When CVL specimens were not available, we analyzed other available specimen types: vaginal swabs (n = 2), cytobrushes (n = 2), or TearFlo paper strips (n = 1).

Sample processing

Specimens were thawed on ice and mixed by pipetting. Up to 800 µl of CVL were centrifuged for 10 minutes at 7,000xg at 4° C. The supernatants were used immediately for virus-like particle enrichment and virome sequencing, and the pellets were frozen at 80° C until used for bacterial microbiome sequencing. PBS negative controls (800 µl) were processed alongside the samples through microbiome and virome sequencing to assess contamination. In the case of vaginal swabs, cytobrushes, or TearFlo strips, 1 mL of 1X PBS was added to the specimen and vortexed for 30 seconds, and 800 µl was then processed in the same way as the CVLs.

Bacterial microbiome sequencing

Frozen bacterial pellets (CVL samples and 1X PBS controls) were thawed on ice. DNA was extracted with the DNeasy PowerSoil Pro kit (Qiagen; Germantown, MD), following the manufacturer’s instructions modified to include 775 µl of Solution CD1, 20 minutes bead beating in a TissueLyser II, 3 minutes of centrifugation at 16,000 g to dry the spin column membrane, and elution in 60 µl of Solution C6. DNA concentrations were measured on a NanoDrop spectrophotometer, and 5–25 µl of DNA (approximately 100–500 ng DNA) per sample were used as input for the Illumina DNA Prep library kit, following the manufacturer’s protocol with 6 cycles for the indexing PCR. Two µl of ZymoBIOMICS Microbial Community DNA Standard was used to build a positive control library. Libraries were pooled, quantified with a Qubit fluorometer and sequenced (paired-end, 2x150 bp) on a NextSeq 1000/2000 sequencer (Illumina; San Diego, CA). Depending on the number of reads obtained per library, some libraries were re-sequenced and the results concatenated to obtain a minimum of 10 million read pairs per sample.

Virus-like particle enrichment and virome sequencing: CVL supernatants and PBS controls were passed through a 0.2-µm pore filter to exclude bacterial and eukaryotic cells. Filtrates were treated with benzonase and Baseline-ZERO DNase to digest non-encapsidated nucleic acids. Total nucleic acids (TNA) were then extracted from 500 µl of the treated filtrates on the bioMérieux eMAG, with PBS added as needed if filtrates were below the required volume. Viral DNA was amplified by multiple displacement amplification (MDA) using Phi29 DNA polymerase (GenomiPhi V2 (Cytiva; Marlborough, MA)), with the kit’s lambda phage DNA for a positive control. The amplified DNA was variably diluted (up to 1:10 depending on band brightness after gel electrophoresis) and used as input for the Illumina DNA Prep library kit, with 6 cycles for the indexing PCR. Libraries were individually quantified with a Qubit fluorometer, pooled and sequenced on an Illumina NextSeq 2000 sequencer (paired-end, 2x150 bp). Separately, we performed sequence-independent amplification (SIA) of viral RNA and DNA by first-strand cDNA synthesis with SuperScript IV reverse transcriptase (ThermoFisher; Waltham, MA) and primers 5'-GTTTCCCAGTCACGATCNNNNNNNNN-3' and 5'-GTTTCCCAGTCACGATC-3', second strand DNA synthesis with DNA Polymerase I Klenow fragment (New England Biolabs; Ipswich, MA), and amplification with AccuPrime Taq high fidelity polymerase (ThermoFisher; Waltham, MA) and primer 5'-GTTTCCCAGTCACGATC-3'. We used cell-culture derived rhinovirus B14 RNA as a positive control. The SIA products were used undiluted as input for the DNA Prep library kit (Illumina; San Diego, CA) and sequenced in the same way as the MDA products.

Virus-like particle enumeration by fluorescence microscopy: Unextracted, nuclease-treated filtrates remaining from the virus-like particle (VLP) enrichment were frozen at 80° C until processed for VLP enumeration by fluorescence microscopy, using an adaption of previously described methods^70,71. Fifteen-100 µl of CVL filtrate and PBS negative controls were mixed with 10 mL sterile SM buffer and filtered onto a 0.02-µm pore, 25 mm diameter Anodisc aluminum oxide filter (Cytiva; Marlborough, MA). De-identified, remnant stool specimens (ASU IRB record STUDY00011967) processed similarly to the CVL samples (resuspension in SM buffer, 0.45-µm and 0.2-µm filtration, and nuclease treatment) and myxoma virus with a known titer were used as positive controls. Filters were dried, stained with 25X SYBR Gold (ThermoFisher; Waltham, MA), and mounted on slides in 15 µl of glycerol/PBS buffer with 0.1% ascorbic acid (working stock: 1 mL glycerol, 980 µl PBS and 20 µl 10% ascorbic acid). Slides were visualized at 60X with oil immersion on an AX R confocal microscope (Nikon; Tokyo, Japan) with a 488-nm wavelength laser. Five fields were randomly selected per sample, and NIS-Elements software (Nikon; Tokyo, Japan) was used to count particles between 50 nm and 500 nm in size. VLP counts for each sample were averaged across the number of fields viewed and used to calculate VLP concentrations in the CVL filtrates and controls.

Bacterial microbiome analysis

Bacterial sequencing reads were quality-filtered, sequencing adapters, PhiX and human sequences removed, reads deduplicated and paired reads merged using BBTools (v. 38.38)⁷². Taxonomy was assigned to the clean reads with KrakenUniq (v. 0.7.3)⁷³ using a custom database (RefSeq archaeal, bacterial, viral, plasmid, and fungal genomes, plus the GRCh38 human genome; downloaded in January 2022). To reduce false positives, we masked taxa that were assigned fewer than 1,800 unique kmers per million reads. For downstream analyses, we parsed the read counts at the species level, i.e., we included reads that could be classified at the species level or lower and excluded reads that could only be classified at the genus level or higher. We excluded probable false positive taxa (Cyanobacteria/Melainabacteria and Rickettsiales), as well as eukaryotes, archaea, viruses, and plasmids. Contaminant bacterial taxa were identified using decontam (R package v. 1.18.0)⁷⁴ and removed at a threshold of 0.1 (Acidovorax sp. KKS102, Streptococcus equi, and Cutibacterium acnes). To normalize for different sequencing depths between samples, species read counts were expressed as number of reads per 100,000 clean reads in a sample, or as relative abundance where appropriate.

Bacterial alpha diversity was measured by species richness, calculated as the number of taxa present in each sample, and Shannon index, calculated using vegan (R package, v. 2.6-4)⁷⁵. Bacterial beta diversity was measured by Bray-Curtis dissimilarity (weighted) and Sorensen dissimilarity (unweighted), calculated with vegan using species relative abundance and species presence-absence, respectively. PCoA analyses were conducted using phyloseq (R package, v. 1.42.0)⁷⁶. Microbiome community groups were defined by kmeans clustering using the R stats package (v. 4.2.2)⁷⁷ on species relative abundance with k = 6. MaAsLin2⁵⁰ was used on species relative abundance to identify species differentially associated either with ART duration or shedder/non-shedder status.

For functional profiling, we assigned reads to bacterial gene families and metabolic pathways using HUMAnN3⁷⁸. Contaminant identification and beta diversity analyses were performed as described above, using normalized gene family and metabolic pathway read counts or relative abundance for input as appropriate. We used MaAsLin2 to identify gene families and pathways associated with sample metadata and microbiome community clusters.

Virome analysis

MDA and SIA virome data were analyzed separately. For the MDA data, sequencing adapters were trimmed with Cutadapt (v. 4.0)⁷⁹. Reads were quality-filtered, PhiX and human reads removed, reads deduplicated, and paired reads merged with BBTools. Contigs were assembled with metaSPAdes (v. 3.15.4)⁸⁰ and deduplicated with cd-hit-est (v. 4.8.1)⁸¹ at 95% identity across 95% of contig length. We retained contigs ≥500 bases in length for further analysis.

Identification of viral contigs

We used Cenote-Taker2 (v. 2.1.5)⁸² with the standard database and VirSorter2 (v. 2.2.4)⁸³ to identify candidate viral contigs. We queried the candidate viral contigs against the NCBI NT database (downloaded June 2020) with megablast (BLAST + v. 2.13.0)⁸⁴ and removed contigs with ≥95% nucleotide identity and query coverage to human sequences. We used CheckV (v. 1.0.1)⁸⁵ to check for false positive, non-viral sequences. We considered contigs to be viral if they were classified as complete, high quality or medium quality by CheckV; or if they were low-quality with at least three viral genes and no more than one host gene; or if they were classified as proviruses. For contigs classified as proviruses, we kept only the viral regions and discarded flanking host sequences. We discarded contigs flagged by CheckV as having high kmer frequency, or which had been classified as bacteria by Cenote-Taker2 during candidate contig identification. To assign taxonomy, we queried the viral contigs against a custom viral protein database (NCBI RefSeq and neighbor viral sequences, downloaded January 2023) using blastx (e-value 1x10^− 3) and used the taxonomy of the best hit in the database. We discarded contigs classified as Mimiviridae as likely false positives. We used Bowtie 2 (v. 2.5.1)⁸⁶ to map the quality-filtered sample reads against the viral contigs. Contaminant contigs were identified using decontam (threshold = 0.1) and removed (six contigs). For downstream analyses, we kept samples that had at least 1,000 viral reads (116 of 125 samples (93%)).

Ecological analyses

Virome alpha (Shannon index and richness) and beta diversity metrics (Bray-Curtis and Sorensen dissimilarity) were calculated and PCoA was conducted as described for the bacterial microbiome, using viral contig relative abundance for weighted analyses and presence-absence for unweighted analyses. Kmeans clustering was conducted on virus family relative abundance in R as described for the bacterial microbiome, with k = 3. MaAsLin2 was used at different taxonomic levels (family, genus, and contig) to identify viral taxa associated with ART duration or shedder/non-shedder status.

Anellovirus phylogenetic analysis: We used Geneious Prime (v. 2023.0.4, http://www.geneious.com) to predict open reading frames (ORFs) from the anellovirus contigs called by blastx during taxonomy assignment. ORFs ≥1,000 nucleotides in length were extracted, manually curated, translated to amino acids and aligned with reference Alphatorquevirus, Betatorquevirus, Gammatorquevirus, Hetorquevirus, Omegatorquevirus, Epsilontorquevirus and Zetatorquevirus sequences retrieved from GenBank (accessions in Supplementary Table 6) using MAFFT (v. 7.505)⁸⁷ with auto algorithm selection. The alignment was trimmed with trimAl (v. 1.4.rev15)⁸⁸ using the -gappyout option and a maximum likelihood phylogeny was constructed with IQ-TREE (v. 2.0.3)⁸⁹ with the model LG + F + G4 and 1,000 bootstrap replicates. The phylogeny was visualized and rooted in FigTree (v. 1.4.4)⁹⁰ using Epsilontorquevirus/Zetatorquevirus as the outgroup, and contigs were assigned to genera based on clade location. Phylogeny-based genus assignments, rather than blastx-based genus assignments, were used for anellovirus contigs in genus-level virome analyses (e.g., MaAsLin2 analyses). For anellovirus UniFrac distance, a separate phylogeny was built with the contig ORF1 sequences and a single Zetatorquevirus reference sequence, using the method described above. This phylogeny was re-rooted in R using ape (v. 5.6-2)⁹¹, the reference sequence was pruned using phyloseq (v. 1.42.0)⁷⁶, and weighted and unweighted UniFrac distance were calculated with phyloseq using the pruned phylogeny and anellovirus contig relative abundance. PCoA analyses were performed using weighted and unweighted UniFrac distance based on the pruned phylogeny and anellovirus contig relative abundance or presence-absence.

Papillomavirus phylogenetic analysis

We used Geneious Prime to predict ORFs from the papillomavirus contigs called by blastx during taxonomy assignment. We mapped contig ORFs ≥1,000 nucleotides in length against reference Alphapapillomavirus, Betapapillomavirus, Gammapapillomavirus, Mupapillomavirus and Nupapillomavirus genomes and annotated the ORFs based on mapping to the reference genes. After manual curation, concatenated contig E1-E2-L2-L1 nucleotide sequences and reference Alphapapillomavirus, Betapapillomavirus, Gammapapillomavirus, Mupapillomavirus and Nupapillomavirus sequences retrieved from RefSeq (accessions in Supplementary Table 7) were aligned with MAFFT with auto algorithm selection, trimmed with trimAl with the -gappyout option, and a maximum likelihood phylogeny constructed with IQ-TREE with the model GTR + G and 1,000 bootstrap replicates. The phylogeny was visualized in FigTree and rooted using Mupapillomavirus/ Nupapillomavirus as the outgroup. Contigs were assigned to genera based on their clade, and phylogeny-based genus assignments were used for papillomavirus contigs in genus-level virome analyses. Papillomavirus beta diversity analyses were performed in the same way as for anelloviruses, with a separate phylogeny built from the contig E1-E2-L2-L1 sequences and a single Mupapillomavirus reference.

SIA virome analysis

SIA virome sequencing reads were processed through the same quality control pipeline as described for the MDA virome reads. SIA reads were then queried against the viral RefSeq and neighbor sequences database using blastx (e-value 1x10^− 3) and parsed with MEGAN (community edition, v. 6.24.21)⁹².

Alphatorquevirus qPCR assay

Forward primer 5’-GTGCCGIAGGTGAGTTTA-3’, reverse primer 5’-AGCCCGGCCAGTCC-3’ and probe 5’-[6-FAM]-TCAAGGGGCAATTCGGGCT-[TAMRA]-3’ were used with Applied Biosystems 2X TaqMan Fast Universal PCR Master Mix and 5 µl of sample total nucleic acid (TNA) per 25-µl reaction. Reactions were performed on an Applied Biosystems QuantStudio 7 Flex Real-Time PCR system (ThermoFisher; Waltham, MA), with reaction conditions of 95° C for 20 seconds and 40 cycles of 95° C for 1 second and 60° C for 20 seconds.

Statistical analysis

All statistical analyses were conducted in R (v. 4.2.2)⁷⁷. P values < 0.05 were considered statistically significant. Where appropriate, corrections were made for multiple comparisons using the Benjamini-Hochberg method⁹³. We corrected for repeated sampling using participant ID as a random effect, where appropriate. PERMANOVA was performed using the adonis2 function in the vegan package (v. 2.6-4)⁷⁵. Multinomial logistic regression was done with the mblogit function in the mclogit package (v. 0.9.6)⁹⁴. Differential abundance analysis was conducted using MaAsLin2 (v. 1.12.0)⁵⁰. Linear mixed modeling was done using the nlme package (v. 3.1–160)⁹⁵. Details of statistical models can be found in the GitHub code repository for this study. Plots were generated using R, GraphPad Prism (v. 9.3.1, www.graphpad.com) and FigTree.

Data availability

De-identified sequencing reads with human sequences removed have been deposited to the NCBI Sequence Read Archive under PRJNA1077994.

Code availability

Code used in this study can be found at https://github.com/ASU-Lim-Lab/CVL-Microbiome-HIV-Discordant-Shedding.

ACKNOWLEDGEMENTS

The authors would like to acknowledge the men and women who participated in this study and the work of the clinical study personnel (Carmela Ganoza, Marcela Rodriguez, Angie Roldan) who made this study possible. Additionally, we would like to acknowledge the technical support from Kelli Kraft, Jillian Legard, Jennifer McKernan-Mullin, and LaRinda Holland. Mara Fredenberg performed the VLP microscopy. We thank Dr. Bertram Jacobs (ASU) for providing the rhinovirus B14 stock sample and Dr. Masmudur Rahman (ASU) for providing the myxoma virus stock sample.

Funding: This work was supported by the National Institute of Allergy and Infectious Disease at the National Institutes of Health. R01 AI071212 and R01 AI091550 awarded to L.M.F.; UM1-AI-106701 awarded to R.W.C.; UW/FHCRC CFAR P30-AI-027757; and R00DK107923 awarded to E.S.L.

AUTHOR CONTRIBUTIONS

E.A.K: data curation, investigation, formal analysis, writing – original draft, writing – review & editing.

C.M: Clinical set-up for the trial, medical expertise as an OB GYN.

J.S: Clinical set-up for the trial, medical expertise with experience in this particular community.

A.L.R: Clinical set-up for the trial, medical expertise with experience in this particular community.

E.T: Clinical set-up for the trial, medical expertise with experience in this particular community.

R.W. C: CLIA certified quantification of HIV in the plasma and CVL.

L.M.F: Design and funding of the parent clinical trial, input on the current study design.

M.E.B: conceptualization, resources, data curation, writing – review & editing.

E.S.L: conceptualization, funding, supervision, writing – original draft, writing – review & editing.

COMPETING INTERESTS

All authors report no competing interests for this study.

Ghosn, J., Taiwo, B., Seedat, S., Autran, B. & Katlama, C. HIV. Lancet392, 685-697, doi:10.1016/s0140-6736(18)31311-4 (2018).
Menéndez-Arias, L. & Delgado, R. Update and latest advances in antiretroviral therapy. Trends Pharmacol Sci43, 16-29, doi:10.1016/j.tips.2021.10.004 (2022).
Joint United Nations Programme on HIV and AIDS. HIV estimates with uncertainty bounds 1990-Present.
Bull, M. et al. Genital Shedding of Human Immunodeficiency Virus Type-1 (HIV) When Antiretroviral Therapy Suppresses HIV Replication in the Plasma. J Infect Dis222, 777-786, doi:10.1093/infdis/jiaa169 (2020).
Homans, J. et al. Permissive and protective factors associated with presence, level, and longitudinal pattern of cervicovaginal HIV shedding. J Acquir Immune Defic Syndr60, 99-110, doi:10.1097/QAI.0b013e31824aeaaa (2012).
Uvin, S. C. & Caliendo, A. M. Cervicovaginal human immunodeficiency virus secretion and plasma viral load in human immunodeficiency virus-seropositive women. Obstet Gynecol90, 739-743, doi:10.1016/s0029-7844(97)00411-0 (1997).
Hart, C. E. et al. Correlation of human immunodeficiency virus type 1 RNA levels in blood and the female genital tract. J Infect Dis179, 871-882, doi:10.1086/314656 (1999).
Cu-Uvin, S. et al. Effect of highly active antiretroviral therapy on cervicovaginal HIV-1 RNA. Aids14, 415-421, doi:10.1097/00002030-200003100-00015 (2000).
Kovacs, A. et al. Determinants of HIV-1 shedding in the genital tract of women. Lancet358, 1593-1601, doi:10.1016/s0140-6736(01)06653-3 (2001).
King, C. C. et al. Prevalence, Magnitude, and Correlates of HIV-1 Genital Shedding in Women on Antiretroviral Therapy. J Infect Dis216, 1534-1540, doi:10.1093/infdis/jix550 (2017).
France, M., Alizadeh, M., Brown, S., Ma, B. & Ravel, J. Towards a deeper understanding of the vaginal microbiota. Nat Microbiol7, 367-378, doi:10.1038/s41564-022-01083-2 (2022).
Ravel, J. et al. Vaginal microbiome of reproductive-age women. Proc Natl Acad Sci U S A108 Suppl 1, 4680-4687, doi:10.1073/pnas.1002611107 (2011).
Gajer, P. et al. Temporal dynamics of the human vaginal microbiota. Sci Transl Med4, 132ra152, doi:10.1126/scitranslmed.3003605 (2012).
Gosmann, C. et al. Lactobacillus-Deficient Cervicovaginal Bacterial Communities Are Associated with Increased HIV Acquisition in Young South African Women. Immunity46, 29-37, doi:10.1016/j.immuni.2016.12.013 (2017).
Liu, F. et al. Comparative metagenomic analysis of the vaginal microbiome in healthy women. Synth Syst Biotechnol6, 77-84, doi:10.1016/j.synbio.2021.04.002 (2021).
Anahtar, M. N. et al. Cervicovaginal bacteria are a major modulator of host inflammatory responses in the female genital tract. Immunity42, 965-976, doi:10.1016/j.immuni.2015.04.019 (2015).
Gudnadottir, U. et al. The vaginal microbiome and the risk of preterm birth: a systematic review and network meta-analysis. Sci Rep12, 7926, doi:10.1038/s41598-022-12007-9 (2022).
Borgdorff, H. et al. Cervicovaginal microbiome dysbiosis is associated with proteome changes related to alterations of the cervicovaginal mucosal barrier. Mucosal Immunol9, 621-633, doi:10.1038/mi.2015.86 (2016).
Łaniewski, P. et al. Linking cervicovaginal immune signatures, HPV and microbiota composition in cervical carcinogenesis in non-Hispanic and Hispanic women. Sci Rep8, 7593, doi:10.1038/s41598-018-25879-7 (2018).
Norenhag, J. et al. The vaginal microbiota, human papillomavirus and cervical dysplasia: a systematic review and network meta-analysis. Bjog127, 171-180, doi:10.1111/1471-0528.15854 (2020).
Brotman, R. M. et al. Bacterial vaginosis assessed by gram stain and diminished colonization resistance to incident gonococcal, chlamydial, and trichomonal genital infection. J Infect Dis202, 1907-1915, doi:10.1086/657320 (2010).
Borgdorff, H. et al. Lactobacillus-dominated cervicovaginal microbiota associated with reduced HIV/STI prevalence and genital HIV viral load in African women. Isme j8, 1781-1793, doi:10.1038/ismej.2014.26 (2014).
Sha, B. E. et al. Female genital-tract HIV load correlates inversely with Lactobacillus species but positively with bacterial vaginosis and Mycoplasma hominis. J Infect Dis191, 25-32, doi:10.1086/426394 (2005).
Thurman, A. R. et al. Vaginal Microbiota and Mucosal Pharmacokinetics of Tenofovir in Healthy Women Using a 90-Day Tenofovir/Levonorgestrel Vaginal Ring. Front Cell Infect Microbiol12, 799501, doi:10.3389/fcimb.2022.799501 (2022).
Dabee, S. et al. Genital microbiota of women using a 90 day tenofovir or tenofovir and levonorgestrel intravaginal ring in a placebo controlled randomized safety trial in Kenya. Sci Rep12, 12040, doi:10.1038/s41598-022-13475-9 (2022).
Mazibuko-Motau, N. et al. Vaginal microbial shifts are unaffected by oral pre-exposure prophylaxis in South African women. Sci Rep12, 16187, doi:10.1038/s41598-022-20486-z (2022).
Klatt, N. R. et al. Vaginal bacteria modify HIV tenofovir microbicide efficacy in African women. Science356, 938-945, doi:10.1126/science.aai9383 (2017).
Donahue Carlson, R. et al. The Female Genital Tract Microbiome Is Associated With Vaginal Antiretroviral Drug Concentrations in Human Immunodeficiency Virus-Infected Women on Antiretroviral Therapy. J Infect Dis216, 990-999, doi:10.1093/infdis/jix420 (2017).
Taneva, E. et al. Vaginal microbiome modulates topical antiretroviral drug pharmacokinetics. JCI Insight3, doi:10.1172/jci.insight.99545 (2018).
Nicol, M. R. et al. Depot Medroxyprogesterone Acetate and the Vaginal Microbiome as Modifiers of Tenofovir Diphosphate and Lamivudine Triphosphate Concentrations in the Female Genital Tract of Ugandan Women: Implications for Tenofovir Disoproxil Fumarate/Lamivudine in Preexposure Prophylaxis. Clin Infect Dis70, 1717-1724, doi:10.1093/cid/ciz443 (2020).
Cheu, R. K. et al. Impact of vaginal microbiome communities on HIV antiretroviral-based pre-exposure prophylaxis (PrEP) drug metabolism. PLoS Pathog16, e1009024, doi:10.1371/journal.ppat.1009024 (2020).
Lantz, A. M. et al. Vaginal microbiome, antiretroviral concentrations, and HIV genital shedding in the setting of hormonal contraception initiation in Malawi. Aids37, 2185-2190, doi:10.1097/qad.0000000000003686 (2023).
Liu, C. M. et al. The Effect of Antiretroviral Therapy Initiation on the Vaginal Microbiome in HIV-Infected Women. Open Forum Infect Dis6, ofz328, doi:10.1093/ofid/ofz328 (2019).
Keita, A. et al. Evidence of HIV-1 Genital Shedding after One Year of Antiretroviral Therapy in Females Recently Diagnosed in Bamako, Mali. Microorganisms9, doi:10.3390/microorganisms9102164 (2021).
Marques, M. P. C. et al. Protease inhibitor effects on prevalence of bacterial vaginosis and human papillomavirus-related lesions. J Obstet Gynaecol Res46, 899-906, doi:10.1111/jog.14245 (2020).
Politch, J. A. et al. Highly active antiretroviral therapy does not completely suppress HIV in semen of sexually active HIV-infected men who have sex with men. Aids26, 1535-1543, doi:10.1097/QAD.0b013e328353b11b (2012).
Gianella, S. et al. Shedding of HIV and human herpesviruses in the semen of effectively treated HIV-1-infected men who have sex with men. Clin Infect Dis57, 441-447, doi:10.1093/cid/cit252 (2013).
Lawn, S. D. et al. Correlation between human immunodeficiency virus type 1 RNA levels in the female genital tract and immune activation associated with ulceration of the cervix. J Infect Dis181, 1950-1956, doi:10.1086/315514 (2000).
Wylie, K. M. et al. Metagenomic analysis of double-stranded DNA viruses in healthy adults. BMC Biol12, 71, doi:10.1186/s12915-014-0071-7 (2014).
Jakobsen, R. R. et al. Characterization of the Vaginal DNA Virome in Health and Dysbiosis. Viruses12, doi:10.3390/v12101143 (2020).
Kaelin, E. A. et al. Cervicovaginal DNA Virome Alterations Are Associated with Genital Inflammation and Microbiota Composition. mSystems7, e0006422, doi:10.1128/msystems.00064-22 (2022).
Madere, F. S. et al. Transkingdom Analysis of the Female Reproductive Tract Reveals Bacteriophages form Communities. Viruses14, doi:10.3390/v14020430 (2022).
Li, Y. et al. Altered vaginal eukaryotic virome is associated with different cervical disease status. Virol Sin38, 184-197, doi:10.1016/j.virs.2022.12.004 (2023).
Siqueira, J. D. et al. Composite Analysis of the Virome and Bacteriome of HIV/HPV Co-Infected Women Reveals Proxies for Immunodeficiency. Viruses11, doi:10.3390/v11050422 (2019).
Wylie, K. M. et al. The vaginal eukaryotic DNA virome and preterm birth. Am J Obstet Gynecol219, 189.e181-189.e112, doi:10.1016/j.ajog.2018.04.048 (2018).
Happel, A. U. et al. Cervicovaginal Human Papillomavirus Genomes, Microbiota Composition and Cytokine Concentrations in South African Adolescents. Viruses15, doi:10.3390/v15030758 (2023).
Happel, A. U. et al. Presence and Persistence of Putative Lytic and Temperate Bacteriophages in Vaginal Metagenomes from South African Adolescents. Viruses13, doi:10.3390/v13122341 (2021).
Sung, H. et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin71, 209-249, doi:10.3322/caac.21660 (2021).
Taylor, L. J., Keeler, E. L., Bushman, F. D. & Collman, R. G. The enigmatic roles of Anelloviridae and Redondoviridae in humans. Curr Opin Virol55, 101248, doi:10.1016/j.coviro.2022.101248 (2022).
Mallick, H. et al. Multivariable association discovery in population-scale meta-omics studies. PLOS Computational Biology17, e1009442, doi:10.1371/journal.pcbi.1009442 (2021).
Beghini, F. et al. Integrating taxonomic, functional, and strain-level profiling of diverse microbial communities with bioBakery 3. eLife10, e65088, doi:10.7554/eLife.65088 (2021).
Christensen, J. K. et al. Prevalence and prognostic significance of infection with TT virus in patients infected with human immunodeficiency virus. J Infect Dis181, 1796-1799, doi:10.1086/315440 (2000).
Shibayama, T. et al. Inverse relationship between the titre of TT virus DNA and the CD4 cell count in patients infected with HIV. Aids15, 563-570, doi:10.1097/00002030-200103300-00004 (2001).
Thom, K. & Petrik, J. Progression towards AIDS leads to increased Torque teno virus and Torque teno minivirus titers in tissues of HIV infected individuals. J Med Virol79, 1-7, doi:10.1002/jmv.20756 (2007).
Schmidt, L. et al. Torque Teno Virus plasma level as novel biomarker of retained immunocompetence in HIV-infected patients. Infection49, 501-509, doi:10.1007/s15010-020-01573-7 (2021).
Honorato, L. et al. The Torque Teno Virus Titer in Saliva Reflects the Level of Circulating CD4(+) T Lymphocytes and HIV in Individuals Undergoing Antiretroviral Maintenance Therapy. Front Med (Lausanne)8, 809312, doi:10.3389/fmed.2021.809312 (2021).
Li, L. et al. AIDS alters the commensal plasma virome. J Virol87, 10912-10915, doi:10.1128/jvi.01839-13 (2013).
Liu, K. et al. HIV-1 Infection Alters the Viral Composition of Plasma in Men Who Have Sex with Men. mSphere6, doi:10.1128/mSphere.00081-21 (2021).
van Rijn, A. L., Roos, R., Dekker, F. W., Rotmans, J. I. & Feltkamp, M. Torque teno virus load as marker of rejection and infection in solid organ transplantation - A systematic review and meta-analysis. Rev Med Virol33, e2393, doi:10.1002/rmv.2393 (2023).
Shkoporov, A. N. & Hill, C. Bacteriophages of the Human Gut: The "Known Unknown" of the Microbiome. Cell Host Microbe25, 195-209, doi:10.1016/j.chom.2019.01.017 (2019).
Blostein, F., Gelaye, B., Sanchez, S. E., Williams, M. A. & Foxman, B. Vaginal microbiome diversity and preterm birth: results of a nested case-control study in Peru. Ann Epidemiol41, 28-34, doi:10.1016/j.annepidem.2019.11.004 (2020).
Carrillo-Ng, H. et al. Variations in cervico-vaginal microbiota among HPV-positive and HPV-negative asymptomatic women in Peru. BMC Res Notes14, 4, doi:10.1186/s13104-020-05422-6 (2021).
Holliday, M. et al. Insights into the Genital Microbiota of Women Who Experienced Fetal Death in Utero. Microorganisms11, doi:10.3390/microorganisms11081877 (2023).
Borgdorff, H. et al. The association between ethnicity and vaginal microbiota composition in Amsterdam, the Netherlands. PLoS One12, e0181135, doi:10.1371/journal.pone.0181135 (2017).
Marconi, C. et al. Characterization of the Vaginal Microbiome in Women of Reproductive Age From 5 Regions in Brazil. Sex Transm Dis47, 562-569, doi:10.1097/olq.0000000000001204 (2020).
Dabee, S., Passmore, J. S., Heffron, R. & Jaspan, H. B. The Complex Link between the Female Genital Microbiota, Genital Infections, and Inflammation. Infect Immun89, doi:10.1128/iai.00487-20 (2021).
Soria, J. et al. Transmitted HIV resistance to first-line antiretroviral therapy in Lima, Peru. AIDS Res Hum Retroviruses28, 333-338, doi:10.1089/aid.2011.0131 (2012).
Reichelderfer, P. S. et al. Effect of menstrual cycle on HIV-1 levels in the peripheral blood and genital tract. WHS 001 Study Team. Aids14, 2101-2107, doi:10.1097/00002030-200009290-00005 (2000).
Bull, M. E. et al. Monotypic human immunodeficiency virus type 1 genotypes across the uterine cervix and in blood suggest proliferation of cells with provirus. J Virol83, 6020-6028, doi:10.1128/jvi.02664-08 (2009).
Patel, A. et al.Virus and prokaryote enumeration from planktonic aquatic environments by epifluorescence microscopy with SYBR Green I., <https://cdr.lib.unc.edu/downloads/tx31qt983?locale=en> (2007).
Liang, G. et al. The stepwise assembly of the neonatal virome is modulated by breastfeeding. Nature581, 470-474, doi:10.1038/s41586-020-2192-1 (2020).
Bushnell, B. BBMap., <https://sourceforge.net/projects/bbmap/> (2019).
Breitwieser, F. P., Baker, D. N. & Salzberg, S. L. KrakenUniq: confident and fast metagenomics classification using unique k-mer counts. Genome Biol19, 198, doi:10.1186/s13059-018-1568-0 (2018).
Davis, N. M., Proctor, D. M., Holmes, S. P., Relman, D. A. & Callahan, B. J. Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data. Microbiome6, 226, doi:10.1186/s40168-018-0605-2 (2018).
Oksanen, J. et al. (2022).
McMurdie, P. J. & Holmes, S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One8, e61217, doi:10.1371/journal.pone.0061217 (2013).
R Core Team. (2022).
Beghini, F. et al. Integrating taxonomic, functional, and strain-level profiling of diverse microbial communities with bioBakery 3. Elife10, doi:10.7554/eLife.65088 (2021).
Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. 201117, 3, doi:10.14806/ej.17.1.200 (2011).
Nurk, S., Meleshko, D., Korobeynikov, A. & Pevzner, P. A. metaSPAdes: a new versatile metagenomic assembler. Genome Res27, 824-834, doi:10.1101/gr.213959.116 (2017).
Li, W. & Godzik, A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics22, 1658-1659, doi:10.1093/bioinformatics/btl158 (2006).
Tisza, M. J., Belford, A. K., Domínguez-Huerta, G., Bolduc, B. & Buck, C. B. Cenote-Taker 2 democratizes virus discovery and sequence annotation. Virus Evol7, veaa100, doi:10.1093/ve/veaa100 (2021).
Guo, J. et al. VirSorter2: a multi-classifier, expert-guided approach to detect diverse DNA and RNA viruses. Microbiome9, 37, doi:10.1186/s40168-020-00990-y (2021).
Camacho, C. et al. BLAST+: architecture and applications. BMC Bioinformatics10, 421, doi:10.1186/1471-2105-10-421 (2009).
Nayfach, S. et al. CheckV assesses the quality and completeness of metagenome-assembled viral genomes. Nat Biotechnol39, 578-585, doi:10.1038/s41587-020-00774-7 (2021).
Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat Methods9, 357-359, doi:10.1038/nmeth.1923 (2012).
Katoh, K. & Standley, D. M. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol30, 772-780, doi:10.1093/molbev/mst010 (2013).
Capella-Gutiérrez, S., Silla-Martínez, J. M. & Gabaldón, T. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics25, 1972-1973, doi:10.1093/bioinformatics/btp348 (2009).
Minh, B. Q. et al. IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era. Mol Biol Evol37, 1530-1534, doi:10.1093/molbev/msaa015 (2020).
Rambaut, A. FigTree v1.4.4, <https://github.com/rambaut/figtree>
Paradis, E. & Schliep, K. ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics35, 526-528, doi:10.1093/bioinformatics/bty633 (2019).
Huson, D. H. et al. MEGAN Community Edition - Interactive Exploration and Analysis of Large-Scale Microbiome Sequencing Data. PLoS Comput Biol12, e1004957, doi:10.1371/journal.pcbi.1004957 (2016).
Benjamini, Y. & Hochberg, Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society. Series B (Methodological)57, 289-300 (1995).
Elff, M. (2022).
Pinheiro, J., Bates, D. & R Core Team. (2022).

There is NO Competing Interest.

CVLExtendedDataFiguresNatMicrobiol.docx
SupplementaryTables.xlsx
Supplementary Tables 1-7

Download PDF

Version 1

posted

You are reading this latest preprint version

Longitudinal cervicovaginal microbiome and virome alterations during ART and discordant shedding in women living with HIV

Status:

Version 1

Abstract

Figures

Introduction

RESULTS

DISCUSSION

METHODS

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1