Rocafort et al. Microbiome (2019) 7:73 https://doi.org/10.1186/s40168-019-0687-5 RESEARCH Open Access Evolution of the gut microbiome following acute HIV-1 infection Muntsa Rocafort1,2 , Marc Noguera-Julian1,2,3, Javier Rivera1,3, Lucía Pastor2,4,5,6, Yolanda Guillén1,2, Jost Langhorst7,8, Mariona Parera1,2, Inacio Mandomando6, Jorge Carrillo1,2, Víctor Urrea1,2, Cristina Rodríguez1,2, Maria Casadellà1,2, Maria Luz Calle3, Bonaventura Clotet1,2,3,9, Julià Blanco1,2,3,5, Denise Naniche4,6 and Roger Paredes1,2,3,9* Abstract Background: In rhesus macaques, simian immunodeficiency virus infection is followed by expansion of enteric viruses but has a limited impact on the gut bacteriome. To understand the longitudinal effects of HIV-1 infection on the human gut microbiota, we prospectively followed 49 Mozambican subjects diagnosed with recent HIV-1 infection (RHI) and 54 HIV-1-negative controls for 9–18 months and compared them with 98 chronically HIV-1- infected subjects treated with antiretrovirals (n = 27) or not (n = 71). Results: We show that RHI is followed by increased fecal adenovirus shedding, which persists during chronic HIV-1 infection and does not resolve with ART. Recent HIV-1 infection is also followed by transient non-HIV-specific changes in the gut bacterial richness and composition. Despite early resilience to change, an HIV-1-specific signature in the gut bacteriome—featuring depletion of Akkermansia, Anaerovibrio, Bifidobacterium, and Clostridium—previously associated with chronic inflammation, CD8+ T cell anergy, and metabolic disorders, can be eventually identified in chronically HIV-1-infected subjects. Conclusions: Recent HIV-1 infection is associated with increased fecal shedding of eukaryotic viruses, transient loss of bacterial taxonomic richness, and long-term reductions in microbial gene richness. An HIV-1-associated microbiome signature only becomes evident in chronically HIV-1-infected subjects. Keywords: Microbiome, HIV-1, acute HIV-1 infection, HIV-1 pathogenesis, AIDS Background like type 2 diabetes, cardiovascular diseases, and frailty There is cumulative evidence that the gut microbiota syndrome in people living with HIV (PLWH) [3, 10–12]. plays an important role in HIV-1 pathogenesis. The All such disorders have been associated with gut micro- gut-associated lymphoid (GALT) and epithelial tissues biota alterations in non-HIV-infected subjects [13–15]. are severely and rapidly damaged following HIV-1 infec- It is thus conceivable that the gut microbiome might tion [1]. The resulting local and systemic inflammation also be involved in the premature aging of PLWH. and loss of CD4+ T cells, chronic immune activation, However, the precise mechanisms governing the inter- and immune deregulation [2–9] are not fully restored play between the host immune system and HIV-1, as with antiretroviral therapy (ART). Chronic inflamma- well as the exact changes occurring in the gut micro- tion, immune activation, and endotoxemia are linked to biome following HIV-1 infection, remain to be defined. the precocious development of aging-related diseases Cross-sectional studies consistently show reduced bac- terial diversity and enrichment in Proteobacteria in * Correspondence: rparedes@irsicaixa.es PLWH, which are linked to lower CD4+ T cell counts Presented at the Conference on Retroviruses and Opportunistic Infections [16], higher inflammation, and increased immune activa- (CROI) 2018, Boston, MA, USA, March 4–7, 2018. tion [5, 17–19]. Previously reported associations between 1 IrsiCaixa AIDS Research Institute, Ctra de Canyet s/n, 08916 Badalona, Catalonia, Spain Prevotella predominance and HIV-1 infection have been 2 Universitat Autònoma de Barcelona, 08193 Bellaterra, Catalonia, Spain shown to be confounded by HIV-1 risk group [19–22]. Full list of author information is available at the end of the article © The Author(s). 2019 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Rocafort et al. Microbiome (2019) 7:73 Page 2 of 17 Using fecal shotgun metagenomics, we previously found findings are analyzed in the context of extensive clinical and that nadir CD4+ T cell count, the main predictor of clin- microbiological metadata to rule out potential confounders ical complications in PLWH, is also a major independent and are associated with a comprehensive panel of immune predictor of reduced microbial gene richness and gut mi- markers, including inflammation, activation, gut barrier in- crobial shifts in this population [23]. Such shifts include tegrity, and bacterial translocation. Finally, a whole metage- depletion of methanogens, sulfate-reducing bacteria, and nomic analysis is performed on a subset of patients to gain other oxygen-sensitive syntrophic microbes, coupled with further insights into gut microbial species and functional dy- enrichment in ROS/RNS-resistant microbes like Bacter- namics following HIV-1 infection. oides and Proteobacteria. Longitudinal studies in non-human primates, which Results allow to control for confounders affecting human studies Study subjects have consistently reported increases in taxa from Out of 4011 subjects screened for fever-like illness or phylum Proteobacteria, such as Actinobacillus spp. and undergoing voluntary HIV-1 testing in the Manhiça Dis- Aggregatibacter spp., as well as in potential pathogens trict Hospital, Mozambique, 85 (2.1%) fulfilled criteria for from the Mycoplasmacetae family and Staphylococ- recent HIV-1 infection (RHI). Forty-nine of them plus 55 cus spp. following simian immunodeficiency virus (SIV) HIV-negative (NEG) individuals consented to participate infection [24, 25]. However, no consistent changes are in this microbiome study (Additional file 1: Table S1). In observed in other major phyla like Bacteroidetes or Fir- addition, 98 chronically HIV-1-infected subjects, 27 of micutes, suggesting that gut microbiome changes after them receiving ART (CHI_ART), mainly consisting on SIV infection are pleiotropic [24]. Furthermore, in NVP, 3TC and AZT, and 71 ART-naive (CHI_noART), non-human primate models, pathogenic but not provided a single fecal and blood sample for non-pathogenic SIV infection is linked to an expansion cross-sectional comparisons with RHI and NEG (Add- of the enteric virome, including Adenovirus, Picorna- itional file 1: Figure S1). For such cross-sectional compari- virus, and Anellovirus, which is associated with gut epi- son, the longitudinal follow-up of RHI was divided in two thelial damage, tissue inflammation, and gut leakiness periods: before and after the first 6 months of follow-up [26]. (RHI ≤ 6 and RHI > 6), and only data from the first sam- In the only longitudinal study available to date in ple available from each period, as well as the first sample humans [27], the gut microbiome of 59 subjects diag- available from each NEG subject, were included. nosed with HIV-1 infection during Fiebig I–IV stages At the first study visit (baseline, Table 1) chronically was enriched in phylum Bacteroidetes and depleted in HIV-1-infected subjects were older, had slightly higher phyla Firmicutes and Proteobacteria compared to 26 BMI scores, and were less likely to have suffered malaria HIV-negative controls. Antiretroviral treatment (ART) in the previous month than the remaining participants initiation was associated with relative increases in Fuso- (Table 1). As expected, compared with HIV-negative in- bacteria, Proteobacteria, and Tenericutes and decreases dividuals, RHI subjects had lower CD4+ T cell, higher in Bacteroidetes and Firmicutes. Such findings are in CD8+ T cell, and lower platelet counts. There were no contrast with non-human primate data [24, 25] which significant differences between groups in the prevalence report opposite dynamics regarding Proteobacteria, as of hepatitis B virus, syphilis, and multiple other micro- well as with data from cross-sectional human studies [5, bial determinations. Longitudinally, the RHI group 18, 28, 29]. Moreover, they provide limited information showed statistically significant decreases in CD8+ T cell below the phylum taxonomic level, which is essential to counts and viral load, together with a significant increase understand the pathophysiology of gut dysbiosis in HIV in CD4+/CD8+ ratio, during the first 6 months after infection and to devise future interventions on the gut HIV-1 acquisition (Additional file 1: Figure S2). microbiome. The initiation of ART shortly after diagno- sis did not enable a more prolonged evaluation of the 16S rRNA gene sequencing gut microbiome dynamics following HIV-1 infection. Fi- Alpha diversity nally, the study did not provide insights as to the role of Genus richness and diversity were not different between intestinal virus expansion after HIV-1 infection. cross-sectional comparison groups except for the ACE Here, we present the results of a prospective and con- richness numeric parameter, which was significantly trolled cohort study investigating the evolution of the gut lower in HIV-negative subjects than in the remaining bacteriome and of several clinically relevant eukaryotic vi- groups (Fig. 1a). Longitudinally, richness and diversity ruses following HIV-1 infection in Mozambique. Longitu- increased during the first 6 and 4 months in both RHI dinal microbiome data are compared with cross-sectional and NEG subjects, respectively (Fig. 1b), with similar samples from treated and non-treated chronically HIV-1-in- slopes between the two groups. Only gender was signifi- fected subjects from the same setting. Gut microbiome cantly associated to differences in microbial richness and Rocafort et al. Microbiome (2019) 7:73 Page 3 of 17 Table 1 Baseline subject’s characteristics Recent HIV-1 infection Chronic HIV-1 infection HIV-Negative p value ART-naive On ART RHI vs neg Overall N = 202 49 71 27 55 Female gender 32 (65%) 52 (73%) 15 (55%) 43 (78%) 0.189 0.154 Age (years) 26 (20; 30) 35 (29; 45.5) 42 (36.5; 46.5) 25 (21; 37) 0.210 < 0.001 Weight (kg) 55 (49; 62) 63 (53.5; 68.5) 61 (60; 68.5) 58.5 (52; 67) 0.093 0.007 Height (cm) 160 (155; 169.7) 160 (154.2; 164) 164 (157.5; 167.5) 163.9 (160; 170.8) 0.192 0.025 BMI 21 (19.1; 23.5) 23.7 (21; 27.6) 23.3 (21.5; 25) 22.2 (18.5; 25) 0.637 < 0.001 Pregnancy 3 (6.1%) 0 0 7 (12.7%) 0.730 0.003 + 3 CD4 T cell/mm 572 (409; 677) 533 (430.5; 726) 460 (366; 594.5) 928 (741; 1142.8) < 0.001 < 0.001 CD8+ T cell/mm3 1229 (740; 1683) 988 (763.5; 1304) 834 (657.5; 1130) 591 (387.5; 689.7) < 0.001 < 0.001 CD4+/CD8+ 0.42 (0.23; 0.70) 0.56 (0.38; 0.79) 0.52 (0.33; 0.88) 1.65 (1.21; 2.24) < 0.001 < 0.001 HIV-1 RNA at timepoint 1 (copies/ 108,220 (31,898; 273,650) 33,400 (7288; 105,200) 75 (75; 272.5) – – < 0.001 mL) HIV-1 RNA at screening (copies/ 2,752,700 (179,220; 15,078,000) – – – – – ml) Time since HIV-1-diagnosis, days – 1344 (977; 1737) 1814 (1463; 2302) – – 0.018 Time on ART, days – – 1257 (678; 1640) – – – Fiebig stage at screening* III 24 (49%) – – – – IV 7 (14.2%) – – – – V 7 (14.2%) – – – – VI 11 (22.6%) – – – – White blood cells 5.8 (4.4; 6.7) 5.2 (4.2; 6.4) 3.7 (3.2; 5) 5.2 (4.4; 6.5) 0.413 < 0.001 Hemoglobin 12.1 (10.9; 13.1) 11.8 (10.4; 12.7) 11.4 (10.8; 12.7) 12.3 (11.6; 13.2) 0.409 0.114 Hematocrite 36.9 (33.5; 39.8) 35.6 (32.7; 38.5) 34.4 (33; 38) 37.7 (34.7; 40.1) 0.320 0.044 Platelets 184.5 (140.7; 230.7) 205 (174.5; 258) 192 (173.5; 216) 226 (185; 265.5) 0.004 0.028 Hepatitis B 4 (8.1%) 2 (2.8%) 3 (11.1%) 1 (1.8%) 0.173 0.125 Syphilis 2 (4.1%) 0 1 (3.7%) 4 (7.3%) 0.682 0.089 Diarrhea, previous week 7 (14.3%) 3 (4.2%) 1 (3.7%) 4 (7.3%) 0.343 0.214 Temperature (°C) 36.3 (36.2; 36.5) 36.4 (36.2; 36.4) 36.4 (36.2; 36.4) 36.2 (36.1; 36.5) 0.956 0.509 Fever, previous 24 h 7 (14.3%) 3 (4.2%) 1 (3.7%) 3 (5.4%) 0.186 0.186 Ritchie test 2 (4.1%) 2 (2.8%) 0 1 (1.81%) 0.600 0.864 Entamoeba hystolitica 2 (4.1%) 6 (8.4%) 1 (3.7%) 2 (3.6%) 1 0.733 Giardia spp 2 (4.1%) 2 (2.8%) 0 5 (9.1%) 0.438 0.236 Cryptosporidium spp 0 2 (2.8%) 0 0 1 0.635 Clostridium difficile 2 (4.1%) 0 1 (3.7%) 1 (1.8%) 0.602 0.458 Strongyloides 1 (2.0%) 0 0 1 (1.8%) 1 0.636 Malaria, prev. month 7 (14.3%) 1 (1.4%) 0 9 (16.4%) 1 < 0.001 Current malaria test 2 (4.1%) 1 (1.4%) 0 0 – 0.071 Malaria test severity 3 (1.5; 3) 4 (4; 4) NA NA – 0.071 Differences in continuous variables were evaluated using ANOVA test (except for the time since HIV-1 diagnosis, which is evaluated using a Student T test). Differences in categorical variables were tested using the Fisher's test. The statistical significance threshold was set to P = 0.05 *Fiebig stage determination at screening visit is described in [30]. Malaria test was only performed in subjects reporting febrile symptoms. Malaria severity is measured in a scale of 1 to 5, from light to severe Rocafort et al. Microbiome (2019) 7:73 Page 4 of 17 Fig. 1 Microbial genus richness and diversity indices obtained with 16S rRNA gene sequencing. a Median ± IQR values in different cross-sectional comparison groups. Kruskal–Wallis p values are shown at the bottom of each plot. Asterisks (p values *< 0.1 and **< 0.05) highlight statistically significant Tukey post hoc pairwise differences between groups, corrected for multiple comparisons (FDR < 0.05). b, c Linear mixed models of the dynamics of richness and diversity indices in recently HIV-infected (RHI, red) and HIV-negative (NEG, green) subjects over time. Horizontal axes show months after study enrollment. Each dot corresponds to a sample, and samples from the same individual through follow-up are line- connected. Single dots correspond to individuals with no longitudinal follow-up. Thick black lines correspond to the modeled slope of each parameter. Vertical dashed lines show the inflection point at month 6 used for modelling in b and c. Statistically, significant differences from 0 (flat slope) are shown with asterisks. p values *< 0.1 and **< 0.05 diversity considering the baseline dataset, with men variables showed a significant differential effect on gut showing significantly lower values than female. Longitu- microbiome composition groups in PERMANOVA ana- dinally, men showed stronger increases along follow-up lyses (not shown). However, there was a statistically sig- compared to females (data not shown). nificant reduction in genus beta-diversity from RHI < 6 to RHI > 6 (Additional file 1: Figure S5). Microbial genus composition In cross-sectional group comparisons (Fig. 2), CHI_- Using non-metric multi-dimensional scaling (NMDS) on noART were significantly depleted in Akkermansia, Anaero- microbial phyla and genera composition data matrices, vibrio, Bifidobacterium, and Clostridium, relative to NEG there were no significant clustering patterns between individuals. Antiretroviral treatment exposure (CHI_ART) cross-sectional comparison groups, regardless of the was not associated with changes in abundance of such gen- distance metric used (Additional file 1: Figure S4). Parti- era, relative to CHI_noART. Individuals with RHI < 6 tion around medoids (PAM) algorithm did not support months were enriched in Haemophilus and Veillonella microbiome clustering, and none of the metadata relative to HIV-negative subjects. Finally, there was an Rocafort et al. Microbiome (2019) 7:73 Page 5 of 17 Fig. 2 Differences in bacterial genera between groups. Box plots show median (± IQR) abundance of bacterial genera. Bacterial genera named “unclassified” are identified by their closest taxonomic level identification. Only bacterial genera with a significantly different abundance between groups (Kruskal–Wallis p value < 0.05) are shown. Statistically significant post hoc pairwise differences (Tukey post hoc pairwise tests corrected for multiple comparisons, FDR < 0.05) are shown with asterisks. Only the first microbiome measurement obtained RHI < 6, RHI > 6, and NEG was used for cross-sectional comparisons with CHI_ART and CHI_noART enrichment in Odoribacter and depletion in Dorea in RHI > abundance over time, but confirmed the aforemen- 6 months relative to RHI < 6 months. tioned observations regarding reduced Akkermansia, To evaluate the longitudinal evolution of microbial gen- Anaerovibrio, and Clostridium in CHI_noART, rela- era, 7 coabundant microbial genus clusters were identified tive to HIV-negative individuals. Moreover, SP2 was (SP1 to 7, Fig. 3a, b). Only clusters SP1, 4, and 7 showed negatively correlated with markers of microbial significant longitudinal changes during follow-up, and translocation (sCD14), inflammation (IP10), and gut only in subjects with RHI (Fig. 3b). Cluster SP1, which integrity (IgG ASCA) (Fig. 3c). was associated with increased plasma levels of cytokines involved in TNF-mediated innate responses (IL-10, IL-12, Whole shotgun metagenome sequencing IP-10, BAFF, CD27, Fas ligand, and TNF receptor 2), de- Alpha diversity creased during the first 6 months of follow-up (Fig. 3c). Shotgun sequencing data also showed increasing trends Conversely, clusters SP4 and SP7, which included Odori- in microbial richness and Shannon diversity over time bacter, Rikenellaceae unclassified, and Barnesiaellaceae (Additional file 1: Figure S3) in both RHI and NEG unclassified and Clostridiales unclassified, Butyricimonas, groups, suggesting that such changes in ecological pa- Faecalibacterium, and Succinivibrio, respectively, in- rameters were not HIV-1-specific. Differences were not creased during the same period (Fig. 3b). There was no statistically significant, however, due to limited power. evidence, however, that the longitudinal evolution of such Previous publications in European subjects have linked clusters was significantly different between RHI and NEG. low values of gut microbial gene richness to metabolic dys- Clusters SP2 and SP3 did not significantly change in regulation, obesity, immune activation, and inflammation Rocafort et al. Microbiome (2019) 7:73 Page 6 of 17 Fig. 3 (See legend on next page.) Rocafort et al. Microbiome (2019) 7:73 Page 7 of 17 (See figure on previous page.) Fig. 3 Dynamics of bacterial clusters following HIV-1 infection. a Within group co-abundant bacterial genus clusters (SP1 to 7) obtained using 16S rRNA gene sequencing. The color gradient is proportional to the mean of scaled individual relative abundance values (mean = 0, sd = 1) per bacterial genera and study group. Dots show statistically significant differences in genus abundance relative to HIV-negative subjects (NEG). b Linear mixed models of the longitudinal evolution of bacterial clusters in subjects with recent HIV-1 infection (RHI) and HIV-negative (NEG) individuals. Horizontal axes show months after study enrollment. Each dot corresponds to a sample, and samples from the same individual through follow-up are line-connected. Single dots correspond to individuals with no longitudinal follow-up. Thick black lines correspond to the modeled slope of each bacterial cluster. Statistically significant differences from 0 (flat slope) are shown with asterisks. p values *< 0.1 and **< 0.05. c Spearman’s correlation between bacterial clusters and immune markers measured in blood. The color gradient is proportional to the Spearman’s rho value. Only unadjusted statistically significant correlations (p value < 0.05) are shown. CHI_ART, CHI_noART, and first available samples from individuals in any of the RHI < 6, RHI > 6, and NEG groups were used to compute correlation values. Other immune markers measured in blood include IgA, IgM, IgG2, and IgG4 for serological makers; EndoCab IgG and IgA ASCA for gut permeability; FABP2 for bacterial translocation; IL7, IL13, GCSF, RANTES, MIP1 alpha, and beta for T cell function; IFN gamma, TNF alpha, and IL8 for Th1 pro-inflammatory responses; TGF beta for anti-inflammatory responses; CD40 ligand and IL21 for B cell function, Eotaxin, IL5, sCD163, and IL15 for innate cells; CXCL16 and IL1 beta for inflammation; B7H1, PDL2, and IL2R for immune activation; and EGF and VEGF for angiogenesis. Several markers were also measured in feces although none of them showed significant correlations with bacterial clusters: sIgA, ANCA, and ASCA for serological markers; EDNEPX, calprotectin, PMNE, lactoferrin, and S100A12 for neutrophil and eosinophil activation; and HBD2, zonulin, and alpha 1 antitrypsin for enterocyte damage and gut permeability. No correlations were found between bacterial clusters and levels of CD4+ and CD8+ T cell activation, exhaustion, and senescence in blood. [13, 31, 32]. In agreement with previous reports, the gut models only identified significant increases in Eubacter- microbial gene richness also followed a bimodal distribu- ium biforme and Collinsella aerofaciens in RHI and in- tion in Mozambican individuals. Samples with gene counts creases in Dorea unclassified and Ruminococcus flavencis above and below a gene count threshold of 444.219 genes in HIV-negative subjects. Consistent with previous re- were categorized as high (HGC) and low gene counts ports [23], subjects with RHI had higher counts of bac- (LGCs), respectively, (Fig. 4). Whereas 63% (22/35) of the terial antioxidant enzymes like catalase, bacterioferritin, RHI samples were included in the LGC category, only 37% glucose-6-phosphate dehydrogenase, or peroxiredoxin, (7/19) of HIV-negative samples were LGC (Fisher’s test although this shotgun analysis was underpowered to p value 0.003) (Fig. 4a). At the first timepoint, 10/13 identify more granular differences (Additiona1 file 1: (76.9%) of RHI vs 4/8 (50%) of HIV-negative subjects (Fish- Figure S6). er’s test p value 0.172) were already LGC (Fig. 4b). Eight of 13 (61.5%) subjects from the RHI group, but none of the HIV-negative individuals showed low gene count values in Fecal virus shedding all follow-up timepoints. Moreover, 5/8 HIV-negative sub- HIV-1 infection was associated with increased fecal jects (62.5%) showed increases in gene richness during shedding of Adenovirus measured by RT-PCR. This was follow-up. observed in RHI (53% of subjects) but also in chronically HIV-1-infected individuals, regardless of whether they Microbial species composition were ART-naïve (51%) or received ART (44%) (all Like in previous European populations, microbial species p values < 0.05, relative to 20% of HIV-1-negative sub- positively correlated with higher gene richness included jects) (Table 2). Conversely, increased shedding of Cyto- Subdoligranulum spp., methanogenic archaea, several bu- megalovirus and Enterovirus was only observed in tyrate producers from the Ruminococcus, Dorea and Eubac- CHI_noART subjects (6% and 21%, respectively), relative terium genera, and Butyrivibrio crossotus (Fig. 5a). to NEG (2% and 4%, respectively) (all p values < 0.05). Prevotella copri was negatively correlated with gene rich- Human herpes viruses 6A, 6B, and 8 were not detected ness in this study, which is consistent with the compos- in any individual, despite using proper extraction and itional nature of the microbiome in this Prevotella-rich positive PCR controls. African population [33, 34]. An unclassified species from Using 16S rRNA gene sequencing to assess bacterial genus Dorea was the unique taxa negatively correlated to taxa, fecal Adenovirus shedding was linked to increased gene richness in the RHI group. relative abundance of an unclassified Mogibacteriaceae In a longitudinal evaluation, several microbial species genus and lower relative abundance of an unclassified associated with gene richness were depleted in RHI Erysipelotrichaceae genus. Cytomegalovirus shedding compared to HIV-negative individuals (Fig. 5b). In gen- was only linked to lower abundance of family Clostridia- eral, such depletion occurred already at the first time- ceae. Enterovirus shedding was linked to lower abun- point and in some cases (i.e., methanogenic archaea, dances of Parabacteroides, Akkermansia, and an Ruminococcus champanellensis, R. flavencis) did not re- unclassified genus from the bacterial family Rikenella- cover during the available follow-up. Linear mixed ceae as well as to higher abundance of an unclassified Rocafort et al. Microbiome (2019) 7:73 Page 8 of 17 a b Fig. 4 Microbial gene richness in recently HIV-1-infected and HIV-negative subjects using shotgun metagenomics. a The leftmost density plot shows a bimodal distribution of all samples according to their observed gene richness value, which enables their classification into low (LGC) and high gene count (HGC). The rightmost density plot shows that HGCs are enriched in NEG, whereas RHI predominate in LGCs. b Longitudinal evolution of microbial gene richness in RHI (in red) and NEG (in green). Each box represents an individual with its longitudinal follow-up samples (timepoints 1, 4, and 9). The dark gray-colored area represents the LGC zone, whereas light gray-colored area represents the HGC zone genus from the family Elusimicrobiaceae (Additional file during chronic HIV-1 infection and does not resolve 1: Figure S7). with ART. In genus-level analyses, HIV-1 infection is Adenovirus shedding was significantly associated with followed by transient changes in the gut microbial richness higher FABP2 and IgG ASCA levels in blood and to and composition, which are not HIV-specific. An HIV-1-as- lower levels of IL10 and IgG2 (Additional file 1: sociated microbiome signature—featuring depletion of Figure S8). Enterovirus detection in feces was positively Akkermansia, Anaerovibrio, Bifidobacterium, and Clostrid- associated to higher levels of IL13, IL8, FABP2, CXCL16, ium—only becomes evident in chronically HIV-1-infected CD40L, and EGF in blood and zonulin levels in feces subjects. Finer species-level shotgun metagenomics analyses and to lower levels of MIG and sCD1643 in blood. Cyto- suggest an early and sustained depletion of methanogenic ar- megalovirus shedding was associated to higher IgA and chaea and several fiber-consuming butyrate-producers fol- MCP1 and lower levels of IL15 in blood. lowing HIV-1 infection. By using RT-PCR, we purposely performed a targeted Discussion analysis of a discrete set of eukaryotic viruses with This study shows that acute HIV-1 infection is followed known pathogenic potential, rather than a comprehen- by increased fecal adenovirus shedding, which persists sive shotgun virome analysis. Most of the gut virome is Rocafort et al. Microbiome (2019) 7:73 Page 9 of 17 Fig. 5 Evolution of richness-associated microbial species in recently HIV-1-infected vs. HIV-negative subjects using shotgun metagenomics. a Microbial species significantly (p < 0.05) associated with microbial gene richness in both study groups, Spearman’s correlation rho values (red is for positive, blue is for negative correlation). b Relative abundance of richness-associated bacterial species over time (months 1, 4, and 9). p values within each box compare month 1 to 4 area under the curves of subjects with recent HIV-1 infection (RHI) versus HIV-negative individuals (NEG). Asterisks show slope values significantly different from 0 in linear mixed models for each bacteria and study group composed by bacteriophages of uncertain pathogenicity, RT-PCR, we observed that fecal Adenovirus shedding and previous analyses using shotgun sequencing found not only was present in more than half of RHI subjects, no relationship between the gut bacteriophage compos- but also in subjects with chronic HIV-1 infection, being ition and HIV-1 [25] or SIV infection [35]. rare in HIV-negative controls. In contrast, fecal CMV Our findings, however, are consistent with previous and Enterovirus shedding was mostly observed in un- studies in non-human primates, where SIV infection was treated chronically HIV-1-infected individuals, suggest- associated with expansion of the enteric virome with ing that CMV shedding might require more prolonged minimal changes in the bacteriome [26, 35]. Using immune dysregulation to occur. Several intestinal Rocafort et al. Microbiome (2019) 7:73 Page 10 of 17 Table 2 Prevalence of fecal virus shedding Recent HIV-1 Chronic HIV-1 infection HIV-negative infection ART-naive On ART Adenovirus 26/49 (53.2%)** 36/71 (50.7%)** 12/27(44.4%)** 11/55 (20.0%) Cytomegalovirus 3/49 (6.1%) 4/71 (5.6%)* 1/27 (3.7%) 1/55 (1.8%) Enterovirus 1/43 (2.4%) 4/19 (21.1%)* 1/25 (4.0%) 2/45 (4.4%) Human herpes virus 6A, 6B, and 8 0/49 (0%) 0/71 (0%) 0/71 (0%) 0/55 (0%) Numbers (percent) of subjects with detectable virus in feces by qualitative commercial real-time PCR. To avoid ascertainment bias due to different follow-up between groups, only the first fecal sample available for testing is used for comparison *p < 0.1, **p < 0.05; Fisher’s pairwise comparisons relative to HIV-1 negative eukaryotic viral families including Adenovirus, CMV, The gut microbial composition also experienced transient and Enteroviruses promote immune activation in the longitudinal changes, with reductions in bacterial clusters context of HIV-1/SIV infection [26]. Indeed, Adenovirus, associated with TNF-mediated innate immune responses Enterovirus, and CMV fecal shedding were linked to sev- and increases in certain butyrate-producing bacteria like eral immune markers of mucus alteration (FABP2, IgG Odoribacter. Increased abundance of potentially pathogenic ASCA, IgA, and MCP1), epithelial disruption (zonulin Proteobacteria has been reported in chronic HIV-1 infected and EGF), and T cell proliferation, activation, and infil- individuals [5, 18, 28, 29] and non-human primates with tration (CXCL16, IL10, and CD40L), hence reinforcing acute SIV infection [24], and members of the Enterobacteri- their role in the inflammatory response to HIV-1 infec- aceae family have been associated to immune activation tion. The persistence of fecal Adenovirus shedding in and inflammation in vivo and ex vivo [5, 28]. However, ART-treated subjects is worrying for clinicians, because none of the gut microbiome changes observed in our study it might be a relentless source of immune activation dif- was clear enough to be considered HIV-specific. ficult to curtail with available treatments. Of note, This study was performed in a Prevotella-rich African HHV-6 and HHV-8 were not related at all with gut population. The lack of Prevotella to Bacteroides shifts microbiome changes following HIV-1 infection. No dif- or changes in the Firmicutes:Bacteroidetes ratio in longi- ferences in eukaryotic viruses fecal shedding were tudinal follow-up provided further evidence against the observed between males and females. previously reported association between these parame- An intrinsic limitation of our study design is that we ters and HIV-1 infection [20, 21, 36, 37]. We did not did not prospectively screen HIV-negative subjects at collect information on same-sex relationships, a poten- risk for HIV-1 infection. Thus, we missed the very early tial confounder in microbiome studies [22], because this events occurring between HIV-1 infection and diagnosis. was a sensitive issue in local culture and all men re- Bearing this limitation in mind, the “increase—plateau” ported heterosexual sex in previous questionnaires. microbial richness dynamics observed in RHI and NEG There is extensive data demonstrating Prevotella pre- groups could be explained by a return to baseline values dominance in African populations [33, 34], so the effect after a transient richness depletion due to an external in- of same-sex intercourse on the gut microbiome would sult. Microbiome richness kinetics in HIV-infected sub- have been hard to discern in this study. Although signifi- jects mirrored that of CD4+ T cell counts and were cant differences in diet have been described between Af- inversed to HIV-1-RNA kinetics. However, the observation rican and European populations and its role on the gut of similar microbiome dynamics in HIV-negative individuals flora [34], collecting information on diet was not in- shows that such changes were not HIV-1-specific. Of note, cluded in the original project design. Hence, we were most HIV-negative controls in our study were seeking med- unable to test the effect of diet in relation to the gut ical care due to fever and/or mononucleosis-like symptoms microbiota in this study. An interesting collateral, previ- when they were recruited into the study, and this possibly ously unexplored aspect of our study is that the micro- affected microbial richness as well. More sensitive shotgun bial species associated with gene richness in this African metagenomics analyses suggested that early microbial gene population were almost identical to those described in richness reductions resolved earlier in HIV-negative individ- European populations [13, 23]. uals than in HIV-1-infected subjects, but larger studies are The only statistically significant association between the needed to confirm this observation. Previous studies have gut microbiome composition and HIV-1 infection was a linked lower microbial gene richness to lower nadir CD4+ T reduction in Akkermansia, Anaerovibrio, Bifidobacterium, cell counts [32], but such variable was not collected in our and Clostridium in chronically HIV-1-infected subjects study. Nonetheless, we did not find any significant correl- relative to HIV-negative individuals. Interestingly, Bifido- ation between gene richness and most recent CD4+ T cell bacterium and Clostridium have anti-inflammatory ef- counts, neither for the RHI and NEG groups. fects, and have been previously associated with HIV-1 Rocafort et al. Microbiome (2019) 7:73 Page 11 of 17 infection [38–40], so reductions in these genera might with known pathogenic potential. As mentioned before, this contribute to chronic inflammation in HIV-1-infected strategy has advantages and limitations. subjects. Akkermansia plays a key role in improved human fat and sugar metabolism [41] and promotes CD8+ T Conclusions cell-mediated immune responses to anti-cancer immune In conclusion, our study shows that HIV-1 infection is modulators like PD-1 and CTLA4 agonists [42]. Defects followed by increased fecal Adenovirus shedding and by in Akkermansia might thus be related to HIV-1-associated transient, non-HIV-specific changes in the gut bacteriome. metabolic disorders as well as to CD8+ T cell anergy. In Despite early resilience to change, an HIV-1-specific sig- our study, Akkermansia and Clostridiaceae abundance nature in the gut microbiome including depletion of were negatively correlated with markers of microbial Akkermansia, Anaerovibrio, Bifidobacterium, and Clos- translocation (sCD14), inflammation (IP10), and gut integ- tridium—previously associated with chronic inflamma- rity (fecal calprotectin). Interestingly, ART initiation was tion, CD8+ T cell anergy, and metabolic disorders—can be not associated with improved abundance of such mi- identified in chronically HIV-1-infected subjects. Longitu- crobes. However, proper longitudinal studies are still dinal studies addressing the ability of ART to prevent and/ needed to assess the impact of ART on the microbiome. or recover such changes are needed. To our knowledge, this is the most comprehensive prospective study investigating the longitudinal impact of HIV-1 infection in the gut microbiome in a hyper- Methods endemic area such as the Manhiça district in Experimental models and subjects’ details Mozambique, which was previously reported to have Human subjects a 40% HIV-1 prevalence in adult population [43]. It This was a sub-study of a prospective observational cohort also incorporates cross-sectional comparisons with study (The GAMA cohort Study) [30, 45] originally de- chronic HIV-infected subjects, which have been highly signed to identify systemic and gastrointestinal inflamma- informative of the longer-term effects of HIV-1 infec- tion biomarkers of recent HIV-1 infection. The study tion in the gut microbiome. However, this study has a population was enrolled between April 2013 and May number of limitations. 2014 at the Manhiça District Hospital (MDH), district of First, HIV-1-infected subjects were followed more Manhiça, Southern Mozambique. The study was approved often than HIV-negative individuals and there was a sig- by local institutional review boards at the Barcelona Clinic nificant amount of loss to follow-up, particularly among Hospital (2011/6264) and by the Ministry of Health of HIV-negative subjects. This could have led to ascertain- Mozambique (461/CNBS/12). Written informed consent ment biases and affected our ability to detect further was obtained from patients before participation. associations. Following healthy individuals in clinical tri- This microbiome sub-study had two components als for prolonged periods of time in resource-limited set- (Additional file 1: Figure S1): tings is very difficult due to a variety of structural constraints [44]. Linear mixed models can partially ac- (a) First, we conducted a prospective observational cohort count for information loss, but it is evident that lost to study to characterize the longitudinal changes follow-up adds uncertainty to our evaluations and might occurring in the fecal microbiome of individuals with have affected our observations. A second limitation is recent HIV-1 infection (RHI) and HIV-negative con- that, due to the study design, we missed the hyperacute trols (NEG) during at least 9 months after their initial phase of HIV infection. Such phases can only be studied HIV-1 assessment. Adults older than 18 years, who either through very large and costly HIV-negative pro- were residents of the established district surveillance spective screening programs, which are out of our reach, system study area and who presented to the outpatient or through controlled experiments in non-human pri- clinic of the MDH for nonspecific febrile symptoms or mates, which might not always reflect events in humans. voluntary HIV counseling and testing were invited to A third limitation is that most HIV-negative subjects were participate in the study. not completely healthy when they were recruited into the (b) Second, data from this longitudinal cohort was study. However, this helped us recognize that the transient compared to a single evaluation from chronically HIV- changes observed in the gut microbiome of HIV-1-infected 1-infected subjects followed at the outpatient HIV-1 subjects were not, in fact, HIV-1-specific. Fourth, shotgun clinic at the MDH, who were receiving antiretroviral metagenomic data was only available for a subset of pro- therapy (CHI_ART) or not (CHI_noART). For this spectively followed HIV-1-infected and HIV-negative sub- cross-sectional comparison, we divided the longitu- jects and was not available for chronic HIV-1-infected dinal follow-up of RHI in two periods: before and after subjects. Finally, we did not perform a formal comprehen- the first 6 months of follow-up (RHI≤6 and RHI > 6), sive virome analysis but, instead, focused on specific viruses and only included data from the first sample available Rocafort et al. Microbiome (2019) 7:73 Page 12 of 17 from each period, as well as the first sample available HIV-1 RNA levels; hemogram; body temperature; from each HIV-1 negative subject (NEG). whether fever was reported or not within the previous 24 h of sample collection; whether the patient reported Pregnant women were not allowed to enter the study. or not diarrhea the week before sample collection and its severity; and information on pregnancy, antiretroviral Diagnosis of recent HIV-1 Infection therapy, and antibiotic use. Blood was collected by finger prick for HIV rapid antibody The dataset also included information on a malaria testing with Determine HIV 1/2 (Abbott Laboratories, optic microscopic assessment performed the day of sam- Chicago, IL, USA). Positive results were confirmed with a ple collection, whether the participant reported malaria more specific Uni-Gold rapid test (Trinity Biotech Co., within the previous month of sample collection, the Wicklow, Ireland). Individuals with positive HIV-1 ser- presence of occult blood in feces, and the result of ology in both rapid tests were not eligible for enrollment microbiological tests in feces, including Giardia spp., into the study and were referred for clinical management. Cryptosporidium spp., Entamoeba spp., Clostridium dif- Recent HIV infection (RHI) was defined as a negative or ficile toxin, as well as serologies for hepatitis B and C, indeterminate rapid test serology (first test negative or syphilis and Strongyloides stercoralis. first test positive and second test negative) and positive HIV-1 viremia by reverse transcriptase–polymerase chain reaction on frozen plasma (Abbott Real-Time HIV-1 16S rRNA gene sequencing: DNA extraction, amplicon Assay, limit of detection 150 copies/mL). HIV-1 RNA test- generation and sequencing ing was performed by applying a multilevel pooling Fecal DNA was extracted using the PowerSoil DNA Ex- scheme of 10 samples/pool as described before [46]. A traction Kit (MO BIO Laboratories, Carlsbad, CA, USA). subgroup of HIV-negative (NEG) time-matched controls The V3–V4 variable regions from the 16S rRNA gene was selected by computer randomization, which were also (amplicon size expected ~ 460 bp) were PCR-amplified prospectively followed and served as controls for RHI. using the MiSeq rRNA Amplicon Sequencing protocol and primers 16S_F 5′-(TCG TCG GCA GCG TCA GAT Study follow-up GTG TAT AAG AGA CAG CCT ACG GGN GGC Subjects with RHI were seen at months 1, 2, 3, 4, 6, 9, 12, WGC AG)-3′ and 16S_R 5′-(GTC TCG TGG GCT 15, and 18 after study enrollment, whereas HIV-negative CGG AGA TGT GTA TAA GAG ACA GGA CTA CHV controls were seen at months 1, 4, and 9. This was consid- GGG TAT CTA ATC C)-3′. Amplification was per- ered the only follow-up realistically feasible in otherwise formed in 25 reactions containing 2.5 μl of non-diluted healthy subjects in the Mozambican setting. Medical con- DNA template, 0.8 μl of each primer at 10 μM, 8.4 μl of sultation and HIV counseling were provided at each med- DNA and RNA-free water, and 12.5 μl of KAPA HiFi ical visit. Study participants received antibiotic treatment HotStart Ready Mix (containing KAPA HiFi HotStart as required according to their clinical status. They also DNA Polymerase, buffer, MgCl2, and dNTPs, KAPA began ART when needed, according to the Mozambican Biosystems Inc., Wilmington, MA, USA). Thermal cyc- ART guidelines applicable at the time of the study (i.e., ling conditions consisted of an initial denaturation step when CD4+ T cell counts were lower than 250 cells/mm3 (5 min at 95 °C), followed by 30 cycles of denaturation or clinical complications requiring ART initiation oc- (20 s at 98 °C), annealing (15 s at 69 °C), and extension curred). If subjects had to start antibiotics or ART, they (15 s at 72 °C). These were followed by a final extension provided a last fecal sample and the study follow-up was step of 1 min at 72 °C. Once the desired amplicon was immediately terminated. Similarly, the study follow-up confirmed in 1% agarose gel electrophoresis, the ampli- was interrupted if women became pregnant. Only one fied DNA was stored at − 30 °C until library preparation. clinical and analytical assessment was performed in Amplified DNA templates were cleaned from non-DNA chronically HIV-1-infected adults. In this pilot exploratory molecules and Illumina sequencing adapters, and dual study, no formal sample size calculation was performed. indices were attached using the Nextera XT Index Kit (Illumina, Inc.) followed by the corresponding PCR amp- Method details lification program as described in the MiSeq 16S rRNA Demographic and clinical data gene Amplicon Sequencing protocol. After a second Demographic and clinical data were collected in a spe- round of cleanup, amplicons were quantified using the cific questionnaire which included the study group to Quant-iT™ PicoGreen® dsDNA Assay Kit (Invitrogen, which the patient belonged (RHI, NEG, CHI_ART, and Carlsbad, MA, USA) and Nanodrop 1000 (Thermo Sci- CHI_noART); the timepoint of follow-up in months; the entific, DE, USA) and diluted in equimolar concentra- participant’s gender, age, weight, and height; his/her tions (T4 nM) for further pooling. Sequencing was CD4+ and CD8+ T cell counts; the CD4+/CD8+ ratio; performed on an Illumina MiSeqTM platform (Illumina, Rocafort et al. Microbiome (2019) 7:73 Page 13 of 17 Inc.) according to the manufacturer’s specifications to gen- above Q20. Paired filtered reads were used for taxo- erate paired-end reads of 300 base length in each direction. nomic characterization of microbial communities using MetaPhlan2 software with the default parameters [51]. Whole genome sequencing (WGS): library preparation and Entire sets of filtered reads were aligned using Bowtie2 to sequencing the integrated reference catalog of the human gut micro- In addition to 16S rRNA gene sequencing, a subset of 54 biome (IGC) (http://meta.genomics.cn/meta/home) [52]. samples was evaluated using shotgun fecal metagenomic Resulting alignments were then filtered with Samtools sequencing. Subjects chosen for shotgun sequencing [53] so that only reads with an alignment quality above analysis had fecal samples available for testing at least at Q20 were kept. Additionally, chimeric and secondary months 1 and 4. If available, month 9 from these same alignments were also removed. Eventually, a subset of six subjects was also sequenced. Extracted, non-diluted million aligned reads was created per sample alignment to DNA was fragmented with the Nextera-XT Illumina kit ensure that results of gene KEGG (http://www.genome.jp/ (Illumina, Inc.) following the manufacturer’s instruc- kegg/) and metabolic pathway richness were comparable tions. One library of approximately 300 bp clone insert between samples. sizes was constructed per sample. Samples were se- quenced in an Illumina HiSeq sequencer (Illumina, Inc.). Detection of fecal virus shedding by qualitative RT-PCR We selectively evaluated the presence of viruses with hu- 16S rRNA gene sequence analysis man pathogenic potential—i.e., Adenovirus (ADV), Cyto- Raw Illumina MiSeq sequences were filtered using a megalovirus (CMV), Human Herpesvirus (HHV) 6A and minimum quality threshold of Q20 in at least 50% of the 6B, and 8 and Enterovirus (ETV)—in cryopreserved stool bases and a minimum sequence average quality of Q20. DNA and RNA aliquots using commercial RT-PCR kits The 1.9.1 version of Qiime (Quantitative Insights Into and following the manufacturer’s instructions, including Microbial Ecology) software pipeline [47] was used for adjustment of quantification using a sensitivity control taxonomic classification of 16S rRNA gene sequences at 1 copy/μl. For ADV, the Adenovirus R-gene® kit (Bio- contained in each sample, using GreenGenes 13.8 as the Mérieux, Marcy-l’Étoile, France) was used. For CMV, reference database [48]. The first step was to join for- HHV 6A, 6B, and 8, specific probe/primer mixes for ward and reverse reads per each sample using the join_ RT-PCR assays (Virusys Corporation, TaneyTown, MD, paired_ends.py script. The maximum percent differences USA) were used in combination with the TaqMan® Gene accepted in the overlapping region were set to 15. Reads Expression Master Mix (Thermo Scientific, DE, USA) at were then filtered to contain a maximum of 2 ambigu- 1/20 dilution. For ETV, fecal RNA was first extracted ous bases (N) using the split_libraries_fastq.py script. An from RNAlater®-cryopreserved fecal samples using the open-reference OTU picking approach was used to clus- Stool total RNA Purification Kit (Norgen Biotek, Cor- ter reads at 97% sequence similarity and construct oper- poration, Thorold, Canada). Non-diluted fresh RNA ational taxonomic units (OTUs). Briefly, in this aliquots were directly used to perform qualitative pick_open_reference_otus.py approach, reads are clus- RT-PCR testing of enteroviruses using the Enterovirus tered against a reference sequence collection and reads R-gene® Kit (BioMérieux, Marcy-l’Étoile, France). that do not hit any reference sequence are subsequently clustered using a de novo clustering approach. We then Cytokine quantification in plasma used the ChimeraSlayer [49] method in the parallel_i- A total number of 51 cytokines were quantified in dentify_chimeric_seqs.py script to identify chimeric se- plasma samples including general markers of inflamma- quences using a PyNast reference from the GreenGenes tion, cell death and growth factors, angiogenesis, mono- 13.8 database. Finally, make_otu_table.py, filter_align- cyte function and mobility, T and B cell function, ment.py and make_phylogeny.py scripts were used to neutrophil and eosinophil activation, enterocyte damage, create final biom and phylogenetic tree files. intestinal permeability, and microbial translocation. De- tection of IL-1R antagonist, MCP-1, GCSF, IFN gamma, Metagenomic sequence analysis IL-12, IL-13, IL-7, VEGF, MIG, RANTES, Eotaxin, MIP1 Raw HiSeq sequences were filtered using Trimmomatic beta, IP-10, IL-2R, IFN alpha, IL-15, GMCSF, TNF alpha, [50]. Nextera adapters were removed, and reads were IL-1β, IL-2, IL-4, IL-5, IL-6, IL-10, MIP1 alpha, IL-17, trimmed enabling a minimum quality of Q30, a global IL-8, and EGF were determined by using Human Cyto- minimum length of 100 bp and using a sliding window kine Magnetic 30-Plex panel (Invitrogen). Detection of set at a minimum quality of Q20 for each 30-bp-long CD40 ligand and IL21 was performed by using the consecutive segments. Human contamination was re- Bio-Plex Pro Human Th17 cytokine assay (Bio-Rad). De- moved by mapping filtered sequences against the human tection of BAFF, CD27, and TNFR2 was performed by genome and removing reads with an alignment quality using a Human Magnetic Luminex Screening Assay Rocafort et al. Microbiome (2019) 7:73 Page 14 of 17 (R&D). Detection of sCD14, LBP, FABP 2, CRP, sCD163, collapsed to different taxonomic levels using Phyloseq CXCL16, sCD23, B7H1, PD ligand 2, TRAIL, Fas ligand, [56] R package and the corresponding relative abun- TGF beta 1 (R&D); IgG ASCA and IgA ASCA (Orgen- dance values were calculated. tec); EndoCab IgG and EndoCab IgA (Hycult biotech); Differences in overall microbiota composition were and zonulin (Immundiagnostik) was assessed by evaluated using Non-Metric Multidimensional Scaling enzyme-linked immunosorbent assays (ELISA) proce- (NMDS) ordination analysis on Bray-Curtis, Euclidean, dures. All quantifications were performed according to weighted UniFrac and unweighted UniFrac distances as the manufacturer’s instructions. Overflow and under available in Phyloseq [56] and Vegan [54] R packages. limit of detection values were validated for every analyte Permutational analysis of variance using Bray-Curtis dis- as the double and the half of the detection limit, tances as implemented in adonis function from the respectively. Vegan R package was used to test the influence of differ- Additional quantification of general antibody types ent metadata variables on microbiota composition simi- and subtypes was also performed by ELISA using the larity measures. The PAM algorithm implemented in following antibodies: polyclonal affinity pure goat Cluster [57] R package was used to test for the existence anti-human IgA antibody α-Chain Specific (Jackson of clusters of samples based on Bray-Curtis distance Immunoresearch), mouse anti-human IgM antibody matrices computed among samples. (clone G20-127) (BD Pharmingen), mouse anti-human Clustering of bacterial genera into co-abundance IgG1 antibody (clone HP6069) (Life technology), mouse groups was performed using hclust algorithm in R stats anti-human IgG2 antibody (clone G18-21) (BD Pharmin- package with ward.D2 methodology on an Euclidean dis- gen™), mouse anti-human IgG3 Hinge (clone HP6050 ) tance matrix computed from a relative abundance table (SouthernBiotech), mouse anti-human IgG4 antibody of bacterial genera across groups. For each genus and (clone G17-4) (BD Pharmingen™), and polyclonal affini- group, we calculated the mean value of all previously pure F(ab')2 fragment goat anti-human IgG, Fcγ frag- scaled individual abundance percentages (force mean = ment specific. 0 and sd = 1). A phylogenetic distance threshold of 0.8 allowed us to identify 7 different bacterial clusters on Cytokine quantification in feces the resulting dendrogram. Determinations in stool samples were performed by The LEfSe algorithm [58] was used to describe which ELISA commercial assays according to the manufac- bacterial genera were significantly enriched or depleted turer’s instructions for calprotectin, PMN-elastase, zonu- in association with positive and negative detection of lin, EDN/EPX, HBD2, secretory IgA1, α1–antitrypsin, Adenovirus, Enterovirus and Cytomegalovirus in feces. S100A12, and claudin (Immundiagnostik, Bensheim, Comparisons were done using both the soft and the Germany); lactoferrin, ASCA, and pANCA (TechLab, strict statistical criterion regarding eukaryotic viral de- Blacksburg, USA). tection and HIV-1 infection status. To increase sample size in tested groups for the latter variable, samples from Quantification and statistical analyses RHI, CHI_ART and CHI_noART groups were catego- Ecological analysis of gut microbiota rized as “Positive” meanwhile the samples form NEG Bacterial taxonomic richness (observed richness and the remained as “Negative”. numeric richness estimators Chao1 and ACE) as well as diversity/evenness measurements (Shannon and Simp- son indices) were calculated using Vegan [54] and Biodi- Analysis on gut microbiota’s gene richness and metabolic versityR [55] R packages, correspondingly. We modified functions the Simpson index’s formula so that it was computed as Filtered metagenomic sequences were mapped against 1/(1-Simpson concentration index). OTUs present the Integrated Gene Catalog (IGC) [52] using the bwa uniquely in one sample were filtered out to reduce noise. software [59]. Unique alignments with a minimum For each of the samples, a subset of 5000 counts was quality of Q20 were selected for subsequent analyses. randomly selected using the rrarefy function from To be able to compare the microbial gene content Vegan, as representative of the entire sample sequence across samples we used a downsampling size of six set. Samples with less than 5000 counts were not consid- million aligned sequences. ered for this ecological analysis. Gene richness was measured as the total number of different genes present in the sample regardless of their Analysis of gut microbiota composition abundance and length. A minimum of one filtered All reads available were used to assess bacterial compos- mapped sequence was set to consider the presence of a ition. OTUs present in less than 10 counts in 10% of gene. The copy number of each gene was estimated by samples were filtered out. Then, OTU counts were dividing the total reads mapping to a gene divided by Rocafort et al. Microbiome (2019) 7:73 Page 15 of 17 the gene’s length. A gene’s relative abundance was mea- phase were fitted using a random intercept, whereas the sured as its copy number divided by the sum of the total last point from the first LMM was used as the intercept gene copies in the sample. As in previous studies [13], for the second LMM. Separate statistics were then per- gene richness followed a bimodal distribution in a prob- formed for each LMM slope. In both cases, p values in- ability density function. The local minimum between the dicated that the slope was significantly different from 0. two modes was considered the threshold value to classify For the cross-sectional comparisons, intra-patient individuals in two groups according to their microbial comparisons at different time points were tested using gene richness content: high gene counts (HGC) and low paired Wilcoxon Rank sum test, whereas comparisons gene counts (LGC). between two groups were tested using non-paired Wil- coxon rank sum test. Comparisons between multiple Microbial functional profile groups were statistically tested using Kruskal–Wallis and Genes found in shotgun metagenomes downsampled at post hoc pairwise tests corrected for multiple compari- 6 M sequences were associated to one or multiple KEGG sons using Benjamini–Hochberg or false discovery rate categories according to the IGC reference database [52]. as needed. The Spearman’s correlation test was used to To measure KEGG abundances, all genes associated to test relation between numerical variables. the same KEGG category were collapsed and their copy numbers were summed. If a gene was associated to more Additional file than one KEGG category, all categories were considered separately. The relative abundance of KEGG categories Additional file 1: Table S1. Study participant’s flow and samples available for testing. Figure S1. Diagram of study design with the was measured as in the relative gene abundance ap- prospective longitudinal and the cross-sectional components. Figure S2. proach detailed above. To analyze the metabolic path- Longitudinal evolution of CD4+ and CD8+ T cell counts and HIV-1 RNA ways represented in all samples, we fed the HUMAnN levels. Figure S3. Changes in observed richness and Shannon’s diversity using shotgun sequencing. Figure S4. Ordination plot of cross-sectional pipeline [60] with the KEGG copy number table. The sample dataset using gut microbiota composition. Figure S5. Beta- HUMAnN output included the relative abundances of diversity of the fecal microbiomes during the first (RHI < 6) and second 6 metabolic modules and pathways for each subject. months (RHI > 6) following HIV-1 infection. Figure S6. Evolution of reactive oxygen species (ROS)-associated enzymes in recently HIV-1- KEGG functions that were not associated to prokaryotic infected vs. HIV-1-negative subjects. Figure S7. Differences in bacterial or archaeal metabolism according to KEGG catalog in- genera relative abundance between adenovirus-, cytomegalovirus- and formation [61] were manually filtered out, to ensure that enterovirus-positive and -negative groups using LEfSe. Figure S8. Cytokines and immune measurements in blood and feces differentially only microbial functions were analyzed. abundant between adenovirus, enterovirus, and cytomegalovirus-positive and -negative groups. (DOCX 2327 kb) Fecal virus shedding analysis Results derived from RT-PCR testing were used as pres- Acknowledgements ence/absence of viral DNA/RNA in feces. Comparison We thank the members of the IGTP Genomics and Bioinformatics Core Facilities (Maria Pilar Armengol, Lauro Sumoy, and Iñaki Martínez) for their of number of individuals with virus in feces among contribution to this publication. ISGlobal, IrsiCaixa, and IGTP are members of groups at study entry was performed using Fisher’s exact the CERCA Programme, Generalitat de Catalunya. statistical test. To summarize results from tested sam- Funding ples from the same individual at different time points to This study was supported by a grant from the Fondo de Investigaciones a single value per individual, we used two different cri- Sanitarias (FIS), Instituto de Salud Carlos III and the European Regional teria: 1) if the first sample available from that given indi- Development Fund (ERDF) (PI13/02514), the Bill and Melinda Gates Foundation, and through philanthropic donations from the Fundació Glòria vidual was considered positive, that subject was positive, Soler, the People in Red – Barcelona 2016 edition and the Fundació 2) if that given individual had at least 1 sample testing Catalunya La Pedrera. IrsiCaixa is supported by the Red de SIDA RD16/0025/ positive along its follow-up, that subject was 0041, cofinanced by the ISCIII and the ERDF, “Investing in your future.” M.R. is funded through a FI-DGR grant (FI-B00184) from Agència de Gestió d’Ajuts considered positive. Universitaris i de Recerca (AGAUR) at the Secretaria d’Universitats i Recerca del Departament d’Economia i Coneixement de la Generalitat de Catalunya. Statistical analysis M.L.C. is funded through the grant MTM2015-64465-C2-1-R, Spanish Ministry of Economy and Competitiveness, Spain. Longitudinal changes in the RHI and NEG groups were analyzed using linear mixed models. After visual inspec- Availability of data and materials tion of raw data, we allowed for a biphasic modeling of The metadata and raw 16S rRNA gene and shotgun metagenome sequences such longitudinal changes, allowing for an inflection have been deposited and are publicly available at the Bioproject PRJNA450025. point at month 6 for RHI—as this is the accepted dur- Authors’ contributions ation of recent HIV-1 infection [62], and at 4 months for RP, JB, DN, and BC conceived and designed the study. LP and IM recruited the NEG subjects, because the parental GAMA study chose study participants and performed their clinical evaluations. Fecal DNA was extracted, amplified, and sequenced by MR, MCa, and MP under the to follow subjects at 1, 4, and 9 months and no 6-month supervision of MN and RP. MR, JR, YG, VU, and CR performed the bioinformatic data was available. Linear mixed models from the first and statistical analyses, with the supervision of MN, MLC, JB, and DN. RP and LP Rocafort et al. Microbiome (2019) 7:73 Page 16 of 17 performed the inflammation and immune activation analyses under the 14. Koeth RA, Wang Z, Levison BS, Buffa JA, Org E, Sheehy BT, et al. Intestinal supervision of JCa and DN. JB and JL performed the inflammation analyses in microbiota metabolism of l-carnitine, a nutrient in red meat, promotes feces. MR performed the RT-PCR assessments of fecal viruses. MR, MN, and RP atherosclerosis. Nat Med. 2013;19:576–85. wrote the paper, which was reviewed, edited, and approved by all authors. All 15. Tang WHW, Wang Z, Levison BS, Koeth RA, Britt EB, Fu X, et al. Intestinal authors read and approved the final manuscript. microbial metabolism of phosphatidylcholine and cardiovascular risk. N Engl J Med. 2013;368:1575–84. Competing interests 16. Nowak P, Troseid M, Avershina E, Barqasho B, Neogi U, Holm K, et al. Gut The authors declare that they have no competing interests. The sponsors of microbiota diversity predicts immune status in HIV-1 infection. Aids. 2015;29: the study had no role in study design, data collection, data analysis, data 2409–18. interpretation, or writing of the report. 17. Dillon SM, Manuzak JA, Leone AK, Lee EJ, Rogers LM, MD MC, et al. HIV-1 infection of human intestinal lamina propria CD4+ T cells in vitro is enhanced by exposure to commensal Escherichia coli. J Immunol. 2012;189: Publisher’s Note 885–96. Springer Nature remains neutral with regard to jurisdictional claims in 18. Vujkovic-Cvijin I, Dunham RM, Iwai S, Maher MC, Albright RG, Broadhurst MJ, published maps and institutional affiliations. et al. Dysbiosis of the gut microbiota is associated with HIV disease progression and tryptophan catabolism. Sci Transl Med. 2013;5:193ra91. Author details 19. Wu GD, Chen J, Hoffmann C, Bittinger K, Chen Y, Keilbaugh SA, et al. 1 Linking Long-term Dietary Patterns with Gut Microbial Enterotypes. Science. IrsiCaixa AIDS Research Institute, Ctra de Canyet s/n, 08916 Badalona, Catalonia, Spain. 2Universitat Autònoma de Barcelona, 08193 Bellaterra, 2011;334:105–8. Catalonia, Spain. 3Universitat de Vic-Universitat Central de Catalunya, C. 20. David LA, Materna AC, Friedman J, Campos-Baptista MI, Blackburn MC, Perrotta Sagrada Família 7, 08500 Vic, Catalonia, Spain. 4ISGlobal, Barcelona Centre for A, et al. Host lifestyle affects human microbiota on daily timescales. Genome International Health Research (CRESIB), Hospital Clínic–Universitat de Biol. 2014;15:R89. Barcelona, 08036 Barcelona, Catalonia, Spain. 5Institut Germans Trias i Pujol 21. Kelley CF, Kraft CS, de Man TJ, Duphare C, Lee H-W, Yang J, et al. The rectal (IGTP), Hospital Germans Trias i Pujol, Universitat Autonoma de Barcelona, mucosa and condomless receptive anal intercourse in HIV-negative MSM: 08916 Badalona, Catalonia, Spain. 6Centro de Investigação em Saúde da implications for HIV transmission and prevention. Mucosal Immunol. 2016 Manhiça (CISM), 1929 Maputo, Mozambique. 7Integrative Gastroenterology, [cited 2017 Mar 30]; Available from: https://www.nature.com/doifinder/10. Kliniken Essen-Mitte, University of Duisburg-Essen, Essen, Germany. 8Chair for 1038/mi.2016.97 Integrative Medicine and translationale Gastroenterology, Klinikum Bamberg, 22. Noguera-Julian M, Rocafort M, Guillén Y, Rivera J, Casadellà M, Nowak P, et University of Duisburg-Essen, Essen, Bavaria, Germany. 9Infectious Diseases al. Gut microbiota linked to sexual preference and HIV infection. Service, Hospital Universitari Germans Trias i Pujol, Ctra de Canyet s/n, 08916 EBioMedicine. 2016;5:135–46. Badalona, Catalonia, Spain. 23. Guillén Y, Noguera-Julian M, Rivera J, Casadella M, Rocafort M, Parera M, et al. Nadir CD4+ T-cell count strongly predicts gut disbyosis in HIV infection. Received: 17 December 2018 Accepted: 22 April 2019 2017; 24. Glavan TW, Gaulke CA, Santos Rocha C, Sankaran-Walters S, Hirao LA, Raffatellu M, et al. Gut immune dysfunction through impaired innate References pattern recognition receptor expression and gut microbiota dysbiosis in 1. Hunt PW, Sinclair E, Rodriguez B, Shive C, Clagett B, Funderburg N, et al. Gut chronic SIV infection. Mucosal Immunol. 2016,9:677–688. epithelial barrier dysfunction and innate immune activation predict 25. Monaco CL, Gootenberg DB, Zhao G, Handley SA, Ghebremichael MS, Lim mortality in treated HIV infection. J Infect Dis. 2014;210:1228–38. ES, et al. Altered virome and bacterial microbiome in human 2. Appay V, Sauce D. Immune activation and inflammation in HIV-1 infection: immunodeficiency virus-associated acquired immunodeficiency syndrome. causes and consequences. J Pathol. 2008;214:231–41. Cell Host Microbe. 2016;19:311–22. 3. Brenchley JM, Price DA, Schacker TW, Asher TE, Silvestri G, Rao S, et al. 26. Handley SA, Thackray LB, Zhao G, Presti R, Miller AD, Droit L, et al. Microbial translocation is a cause of systemic immune activation in chronic Pathogenic simian immunodeficiency virus infection is associated with HIV infection. Nat Med. 2006;12:1365–71. expansion of the enteric virome. Cell. 2012;151:253–66. 4. Brenchley JM, Schacker TW, Ruff LE, Price DA, Taylor JH, Beilman GJ, et al. 27. Sortino O, Phanuphak N, Chomchey N, Belkaid Y, Davis J, Quinones M, et al. CD4+ T cell depletion during all stages of HIV disease occurs predominantly Acute HIV-infection and human gut microbiome before and after anti- in the gastrointestinal tract. J Exp Med. 2004;200:749–59. retroviral therapy goals background and study rationale results. 2017 5. Dinh DM, Volpe GE, Duffalo C, Bhalchandra S, Tai AK, Kane AV, et al. Available from: https://www.croiconference.org/sites/default/files/posters- Intestinal microbiota, microbial translocation, and systemic inflammation in 2017/214_Sortino.pdf. chronic HIV infection. J Infect Dis. 2015;211:19–27. 28. Dillon SM, Lee EJ, Kotter C V, Austin GL, Dong Z, Hecht DK, et al. An altered 6. Guadalupe M, Reay E, Sankaran S, Prindiville T, Flamm J, McNeil A, et al. intestinal mucosal microbiome in HIV-1 infection is associated with mucosal Severe CD4+ T-cell depletion in gut lymphoid tissue during primary human and systemic immune activation and endotoxemia. Mucosal Immunol. 2014; immunodeficiency virus type 1 infection and substantial delay in restoration 7:983–994. following highly active antiretroviral therapy. J Virol. 2003;77:11708–17. 29. Yang L, Poles MA, Fisch GS, Ma Y, Nossa C, Phelan JA, et al. HIV-induced 7. Klatt NR, Funderburg NT, Brenchley JM. Microbial translocation, immune immunosuppression is associated with colonization of the proximal gut by activation and HIV disease. Trends Microbiol. 2013;21(1):6–13. environmental bacteria. AIDS. 2016;30:19–29. 8. Marchetti G, Tincati C, Silvestri G. Microbial translocation in the 30. Pastor L, Parker E, Carrillo J, Urrea V, Fuente-Soro L, Respeito D, et al. A pathogenesis of HIV infection and AIDS. Clin Microbiol Rev. 2013;26:2–18. cytokine pattern that differentiates preseroconversion from 9. Mehandru S, Poles MA, Tenner-Racz K, Horowitz A, Hurley A, Hogan C, et al. postseroconversion phases of primary HIV infection. J Acquir Immune Defic Primary HIV-1 infection is associated with preferential depletion of CD4+ T Syndr. 2017;74:459–66. lymphocytes from effector sites in the gastrointestinal tract. J Exp Med. 31. Cotillard A, Kennedy SP, Kong LC, Prifti E, Pons N, Le Chatelier E, et al. 2004;200:761–70. Dietary intervention impact on gut microbial gene richness. Nature. 2013; 10. Douek DC. Disrupting T-cell homeostasis: how HIV-1 infection causes 500(7464):585–8. disease. AIDS Rev. 2003;5:172–7. 32. Guillén Y, Noguera-Julian M, Rivera J, Casadella M, Zevin AS, Rocafort M, et 11. Paiardini M, Müller-Trutwin M. HIV-associated chronic immune activation. al. Low nadir CD4+ T-cell counts predict gut dysbiosis in HIV-1 infection. Immunol Rev. 2013;254:78–101. Mucosal Immunol. 2018;12(1):232–46. 12. Sandler NG, Douek DC. Microbial translocation in HIV infection: causes, 33. Ou J, Carbonero F, Zoetendal EG, DeLany JP, Wang M, Newton K, et al. Diet, consequences and treatment opportunities. Nat Rev Microbiol. 2012;10:655–66. microbiota, and microbial metabolites in colon cancer risk in rural Africans 13. Le Chatelier E, Nielsen T, Qin J, Prifti E, Hildebrand F, Falony G, et al. and African Americans. Am J Clin Nutr. 2013;98:111–20. Richness of human gut microbiome correlates with metabolic markers. 34. De Filippo C, Cavalieri D, Di Paola M, Ramazzotti M, Poullet JB, Massart S, et Nature. 2013;500:541–6. al. Impact of diet in shaping gut microbiota revealed by a comparative Rocafort et al. Microbiome (2019) 7:73 Page 17 of 17 study in children from Europe and rural Africa. Proc Natl Acad Sci U S A. 59. Li H, Durbin R. Fast and accurate short read alignment with Burrows- 2010;107:14691–6. Wheeler transform. Bioinformatics. 2009;25:1754–60. 35. Handley SA, Desai C, Zhao G, Droit L, Monaco CL, Schroeder AC, et al. SIV 60. Abubucker S, Segata N, Goll J, Schubert AM, Izard J, Cantarel BL, et al. infection-mediated changes in gastrointestinal bacterial microbiome and Metabolic reconstruction for metagenomic data and its application to the virome are associated with immunodeficiency and prevented by human microbiome. Eisen JA, editor. PLoS Comput Biol. 2012;8:e1002358. vaccination. Cell Host Microbe. 2016;19:323–35. 61. Kanehisa M, Goto S. KEGG: Kyoto Encyclopedia of Genes and Genomes. 36. Noguera-Julian M, Rocafort M, Guillén Y, Rivera J, Casadellà M, Nowak P, et Nucleic Acids Res. 2000;28:27–30. al. Gut microbiota linked to sexual preference and HIV infection. 62. Cohen MS, Gay CL, Busch MP, Hecht FM. The detection of acute HIV EBioMedicine. 2016;5:135–146. infection. J Infect Dis. 2010;202:S270–7. 37. Wu GD, Chen J, Hoffmann C, Bittinger K, Chen Y-Y, Keilbaugh SA, et al. Linking long-term dietary patterns with gut microbial enterotypes. Science (80- ). 2011;334:105–108. 38. Gori A, Tincati C, Rizzardini G, Torti C, Quirino T, Haarman M, et al. Early impairment of gut function and gut flora supporting a role for alteration of gastrointestinal mucosa in human immunodeficiency virus pathogenesis. J Clin Microbiol. 2008;46:757–8. 39. Pérez-Santiago J, Gianella S, Massanella M, Spina CA, Karris MY, Var SR, et al. Gut Lactobacillales are associated with higher CD4 and less microbial translocation during HIV infection. AIDS. 2013;27:1921–31 40. Mutlu EA, Keshavarzian A, Losurdo J, Swanson G, Siewe B, Forsyth C, et al. A compositional look at the human gastrointestinal microbiome and immune activation parameters in HIV infected subjects. PLoS Pathog. 2014;10:e1003829. 41. Schneeberger M, Everard A, Gómez-Valadés AG, Matamoros S, Ramírez S, Delzenne NM, et al. Akkermansia muciniphila inversely correlates with the onset of inflammation, altered adipose tissue metabolism and metabolic disorders during obesity in mice. Sci Rep. 2015;5:16643. 42. Routy B, Le Chatelier E, Derosa L, Duong CPM, Alou MT, Daillère R, et al. Gut microbiome influences efficacy of PD-1-based immunotherapy against epithelial tumors. Science (80- ). 2018;359:91–7. 43. González R, Munguambe K, Aponte JJ, Bavo C, Nhalungo D, Macete E, et al. High HIV prevalence in a southern semi-rural area of Mozambique: a community-based survey. HIV Med. 2012;13:581–8. 44. Schwitters A, Lederer P, Zilversmit L, Gudo PS, Ramiro I, Cumba L, et al. Barriers to health care in rural Mozambique: a rapid ethnographic assessment of planned mobile health clinics for ART. Glob Heal Sci Pract. 2015;3:109–16. 45. Pastor L, Urrea VV, Carrillo J, Parker E, Fuente-Soro L, Jairoce C, et al. Dynamics of CD4 and CD8 T-cell subsets and inflammatory biomarkers during early and chronic HIV infection in Mozambican adults. Front Immunol. 2017;8:1925. 46. Serna-Bolea C, Muñoz J, Almeida JM, Nhacolo A, Letang E, Nhampossa T, et al. High prevalence of symptomatic acute HIV infection in an outpatient ward in southern Mozambique: identification and follow-up. AIDS. 2010;24:603–8. 47. Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. QIIME allows analysis of high-throughput community sequencing data. Nat Methods. 2010;7:335–336. 48. DeSantis TZ, Hugenholtz P, Larsen N, Rojas M, Brodie EL, Keller K, et al. Greengenes, a chimera-checked 16S rRNA gene database and workbench compatible with ARB. Appl Environ Microbiol. 2006;72:5069–72. 49. Haas BJ, Gevers D, Earl AM, Feldgarden M, Ward DV, Giannoukos G, et al. Chimeric 16S rRNA sequence formation and detection in Sanger and 454- pyrosequenced PCR amplicons. Genome Res. 2011;21:494–504. 50. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20. 51. Truong DT, Franzosa EA, Tickle TL, Scholz M, Weingart G, Pasolli E, et al. MetaPhlAn2 for enhanced metagenomic taxonomic profiling. Nat Methods. 2015;12:902–3. 52. Li J, Jia H, Cai X, Zhong H, Feng Q, Sunagawa S, et al. An integrated catalog of reference genes in the human gut microbiome. Nat Biotech. 2014; advance on:834–41. 53. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9. 54. Oksanen J, Blanchet FG, Kindt R, Legendre P, Minchin PR, O’Hara RB, et al. Package ‘vegan.’ R Packag ver 20–8; 2013. p. 254. 55. Roeland Kindt. Package “BiodiversityR.” 2016. 56. McMurdie PJ, Holmes S, Gregory Jordan SC. Package “phyloseq.”; 2016. 57. M Maechler, P Rousseeuw, A Struyf , M Hubert, K Hornik, Studer, M, Roudier P. Package “cluster.” 58. Segata N, Izard J, Waldron L, Gevers D, Miropolsky L, Garrett WS, et al. Metagenomic biomarker discovery and explanation. Genome Biol. 2011; 12(6):R60.