Gut dysbiosis associates with cytokine production capacity in viral-suppressed people living with HIV

Background People living with human immunodeficiency virus (PLHIV) are exposed to chronic immune dysregulation, even when virus replication is suppressed by antiretroviral therapy (ART). Given the emerging role of the gut microbiome in immunity, we hypothesized that the gut microbiome may be related to the cytokine production capacity of PLHIV. Methods To test this hypothesis, we collected metagenomic data from 143 ART-treated PLHIV and assessed the ex vivo production capacity of eight different cytokines [interleukin-1β (IL-1β), IL-6, IL-1Ra, IL-10, IL-17, IL-22, tumor necrosis factor, and interferon-γ] in response to different stimuli. We also characterized CD4+ T-cell counts, HIV reservoir, and other clinical parameters. Results Compared with 190 age- and sex-matched controls and a second independent control cohort, PLHIV showed microbial dysbiosis that was correlated with viral reservoir levels (CD4+ T-cell–associated HIV-1 DNA), cytokine production capacity, and sexual behavior. Notably, we identified two genetically different P. copri strains that were enriched in either PLHIV or healthy controls. The control-related strain showed a stronger negative association with cytokine production capacity than the PLHIV-related strain, particularly for Pam3Cys-incuded IL-6 and IL-10 production. The control-related strain is also positively associated with CD4+ T-cell level. Conclusions Our findings suggest that modulating the gut microbiome may be a strategy to modulate immune response in PLHIV.


Introduction
Human immunodeficiency virus (HIV) infection, as shown in our previous study, induces chronic activation of the innate and adaptive immune systems, leading to a chronic inflammatory state (van der Heijden et al., 2021). Combination antiretroviral therapy (ART) significantly decreases immune activation and systemic inflammation but does not restore the homeostasis in the immune system to that seen in healthy populations (Hileman and Funderburg, 2017). The persistent inflammation, due, in part, to a dysbalanced cytokine system, contributes to a higher risk of non-AIDS-related morbidity in people living with HIV (PLHIV), including cardiovascular disease, neurocognitive disease, and certain HIVrelated cancers (Keating et al., 2012;Hileman and Funderburg, 2017). Previous studies have shown that HIV triggers the production of proinflammatory cytokines [tumor necrosis factor (TNF), interleukin-6 (IL-6), and IL-1] and anti-inflammatory cytokines (IL-10) (Osuji et al., 2018;Planès et al., 2018). Starting ART has been shown to decrease plasma concentrations of IL-10 and IL-6, but the concentrations of TNF and other proinflammatory cytokines remain elevated (Osuji et al., 2018). There are multiple potential causes for this dysregulated cytokine system, such as the effect of HIV, lymphoid tissue damage (Klatt et al., 2013), and, in particular, gut dysbiosis (Crakes and Jiang, 2019).
HIV infection induces significant changes in gut microbial composition and metabolic function (Vaźquez-Castellanos et al., 2018;Vujkovic-Cvijin and Somsouk, 2019). ART can partially restore the HIV-associated gut dysbiosis, but it cannot normalize the gut microbiome to a pattern resembling that of a healthy control (HC) population (Lozupone et al., 2013;Ray et al., 2021). Compared with HCs, long-term treated PLHIV still exhibit decreased alpha diversity, increased abundances of Enterobacteriaceae, and decreased abundances of Bacteroidetes and Alistipes (Crakes and Jiang, 2019) and butyrate-producing bacteria that help maintain healthy gut homeostasis (Pinto-Cardoso et al., 2017). However, no consistent pattern of gut dysbiosis has been defined in long-term treated PLHIV (Tuddenham et al., 2020). One major reason for this is that sexual behavior may have a strong influence on gut microbiota, and most previous studies were not able to control for this factor (Tuddenham et al., 2020). For example, an overrepresentation of Prevotella accompanied with a decrease in Bacteroides in PLHIV was recently found to be due to men who have sex with men (MSM) status rather than HIV infection (Vujkovic-Cvijin et al., 2013;Noguera-Julian et al., 2016;Armstrong et al., 2018;Vujkovic-Cvijin et al., 2020). Considering that a major part of the HIV-infected population in Europe and North America is composed of MSM, the biological importance of this feature needs to be explored.
In long-term treated PLHIV, a link between gut dysbiosis and host cytokine levels has been reported. For instance, in PLHIV with atherosclerosis, class Clostridia was positively correlated with plasma levels of IL-1b and interferon-g (IFN-g) (Ishizaka et al., 2021), whereas coproic acid, a gut bacteria-derived short-chain fatty acid (SCFA), was linked with decreased expression of IL-32 (El-Far et al., 2021). Mechanistically, the bacteria-cytokine association is potentially based on bacterial components and products (Dillon et al., 2012;Dillon et al., 2016;Dillon et al., 2017;Vaźquez-Castellanos et al., 2018). For example, heat-killed Escherichia coli induced a higher production of IL-17 and IFN-g in HIV-exposed mononuclear cells ex vivo (Dillon et al., 2012). In addition, lipopolysaccharide (LPS), a Gram-negative bacterial cell wall component, induced a depletion of CD4 + cells by increasing expression of HIV coreceptor C-C chemokine receptor type 5 on CD4 + cells (Dillon et al., 2016). In addition, butyrate, a product of saccharolytic fermentation of dietary fibers by gut microbiota, decreased gut T-cell activation in an ex vivo human intestinal cell culture model (Dillon et al., 2017). The downregulation of antiinflammatory bacterial pathways, such as SCFA biosynthesis or indole production, also contributes to gut inflammation in longterm treated PLHIV (Vaźquez-Castellanos et al., 2018).
The present study documents a detailed profile of gut microbial composition and function at both species and strain levels using metagenomic sequencing in PLHIV. We then characterize the association of gut dysbiosis in relation to HIV clinical phenotypes and peripheral blood mononuclear cell (PBMC) cytokine production capacity. Notably, the association between cytokine production capacity and gut microbiome has only been studied in healthy populations (Schirmer et al., 2016) and in a limited number of PLHIV (Lozupone et al., 2013). In addition, we control for sexual behavior-related factors in the association analysis. Finally, we identify two P. copri strains with different genetic repertoires that exhibit enrichment in PLHIV and HCs, respectively. The control-related P. copri strain showed stronger associations with the PBMC production capacity of IL-6 and IL-10, as well as the CD4 + T-cell level.

Study cohorts
The HIV cohort used in this study was described in our previous study (van der Heijden et al., 2021), in which we recruited 211 PLHIV from the HIV clinic of the Radboud University Medical Center between December 2015 and February 2017. Caucasian individuals who were 18 years of age or older were on combination antiretroviral therapy (cART) for more than 6 months with an HIV-RNA load greater than or equal to 200 copies/ mL and showed no signs of opportunistic infections or active hepatitis B/C were included. Venous blood was collected in sterile 10-mL ethylenediaminetetraacetic acid (EDTA) and 8-mL serum BD Vacutainer tubes (Becton Dickinson) and processed within 1-4 h. Isolation of PBMCs and monocytes was performed on freshly collected blood by density centrifugation over Ficoll-Paque (VWR) or by the Pan Monocyte Isolation Kit (Miltenyi Biotec), respectively. Study participants self-collected stool at their homes no more than 24 h prior to study visits, using an empty sterile 50-mL tube with a screwcap with integrated spoon inside. Samples were stored in a refrigerator until being brought in for their visits and then were aliquoted and frozen at −80°C. For this study, 143 metagenomic sequencing samples were available. Participants were excluded if they reported any antibiotics usage in the 3 months prior to fecal sample collection. One fecal sample was analyzed per individual. We included 190 age-and sex-matched HCs from the Dutch Microbiome Project (DMP) cohort as the control group . To replicate our findings in an independent cohort, we also included 173 sex-matched HCs from the 500 Functional Genomics Project (500FG) cohort (Schirmer et al., 2016).

Metagenomic data generation and profiling
The same protocol for fecal DNA isolation and metagenomic sequencing was used for both HIV samples and heathy control samples. Fecal DNA isolation was performed using the QIAamp Fast DNA Stool Mini Kit FSK; Qiagen, cat. no. 51604). The stool samples are lysed using the proteinase K and then heated in the temperature of 95°C. Fecal DNA was sent to Novogene to conduct library preparation and perform whole-genome shotgun sequencing on the Illumina HiSeq platform. We filtered out the low-quality reads by trimming the bases with PHRED quality below 30 and discarded the reads with the length below 70 base pairs after trimming. The reads aligning to human genome were removed by mapping the data to the human reference genome (version HCBI37) using KneadData (v0.7.4). After filtering, the average read depth was 26.8 million for 143 HIV samples, 23.1 million for 190 DMP samples, and 25.1 million for 173 500FG samples. Microbial taxonomic and functional profiles were determined using Metaphlan3 (v3.0.7) (Beghini et al., 2021) and HUMAnN3 (v3.0.0.alpha.3) (Beghini et al., 2021). MetaPhlAn3 is a bioinformatic method to determine microbial composition by mapping reads against a database of nearly 1 million unique clade-specific marker genes identified from around 17,000 reference genomes (13,500 bacterial and archaeal, 3,500 viral, and 110 eukaryotic) (https://huttenhower.sph.harvard.edu/ metaphlan3/). The reads identified by MetaPhlAn3 are mapped to species-specific pangenomes with UniRef90 annotations, and the MetaPhlAn3-unclassified reads are translated and aligned to a protein database. Bacteria/pathways present in < 20% of the samples from one cohort were discarded.

Strain profiles and analysis
We used Pangenome-based Phylogenomic Analysis3 (PanPhlAn3) (Beghini et al., 2021) to identify the gene composition at the strain level. A total of 4,998 gene families from the P. copri pangenome were detected across 282 samples from the three cohorts. A Jaccard distance matrix was built according to the presence/absence pattern of gene families. The strain cluster tree was constructed using the R basic function hclust with the hierarchical clustering method "complete". Three strain clusters were defined at a tree height of 0.3. Visualizations were generated using the dendextend R package. Differentially abundant gene families were obtained using a logistic regression model, controlling for sex, age, and read counts. The subsequent Pfam enrichment analysis was conducted using the clusterProfiler R package (v. 3.18.1) (P. copri Pfam annotation from PhanPhlan3), where the p-value for enrichment can be calculated by hypergeometric distribution. In the association analysis between cytokine production and HIV-related parameters and the two P. copri strains, we used linear regression, controlling for age, sex, read counts, and sexual behavior, as well as using false discovery rate (FDR) < 0.1 as the significant threshold. We identified five P. copri peptides with immune function in Immune Epitope Database and Analysis Resource (IEDB) (Vita et al., 2019) (Supplementary Table 25) and checked their abundance across the two strains in PLHIV and HCs using ShortBRED (Kaminski et al., 2015) using linear regression and controlling for age, sex, and read counts.

Microbial compositional and differential abundance analysis
The relative abundance data obtained from MetaPhlan3 were used to calculate bacterial diversity using the vegan R package (v. 2.5-7). Alpha diversity was calculated using the diversity function. The Bray-Curtis distance n-by-n matrix was built using the vegdist function, and, then, the PERMANOVA statistical test was applied on the matrix. We used the principal coordinates analysis method to visualize the dissimilarities of beta diversity between different cohorts. To obtain the differentially abundant bacterial species and pathways between PLHIV and HCs, we first transformed the relative abundance data using centered log-ratio (CLR) transformation, as described before . Bacteria/pathways present in < 20% samples in at least one cohort were then discarded, and the remaining data were inverse-rank-transformed to follow a normal distribution. A linear regression model was then fitted, controlling for body mass index (BMI), smoking status, and read counts. Benjamini-Hochberg correction was used to correct for multiple hypothesis testing (using FDR < 0.05 as the significance threshold). Network of bacterial pathways was conducted using igraph R package (v. 1.2.6), with layout method of "layout_with_fr," and the bacterial pathways are annotated into super pathways based on Metacyc Database. Dysbiosis index (DI) and function imbalance (FI) scores were calculated as the log2 ratio between geometric means of relative abundances of species/ pathways that were enriched in PLHIV (linear regression FDR < 0.05 and HIV cohort beta > 0 in PLHIV vs. HCs) and depleted in PLHIV (linear regression FDR < 0.05 and HIV cohort beta < 0 in PLHIV vs. HCs). DI and FI scores for HCs from the 500FG cohort were calculated using the same species/pathways from the comparison between the HIV and DMP cohorts.

Measurement and analysis of ex vivo PBMC cytokine production, intestinal damage, and monocyte activation
The detailed methods of ex vivo PBMC stimulation, cytokines measurements, and plasma markers measurements have been described before (van der Heijden et al., 2021). Venous blood was collected and processed within 1-4 h in sterile 10-mL EDTA and 8-mL serum BD Vacutainer containers (Becton Dickinson). In short, density centrifugation was performed on freshly collected venous blood to obtain the isolation of PBMCs. The freshly isolated cells were then incubated with different bacterial, fungal, and viral stimuli at 37°C and 5% CO 2 for either 24 h or 7 days. IL-1b, IL-6, IL-1Ra, IL-10, and TNF were determined in the supernatants of the 24-h PBMC or monocyte stimulation experiments using enzyme-linked immunosorbent assay (ELISAs). IL-17, IL-22, and IFN-g were measured after the 7-day stimulation of PBMCs. Cytokine production data for the 500FG cohort were obtained using the same method, and the measurements that overlapped with those in the HIV cohort are summarized in Supplementary Table 9. For the comparison of cytokine production between PLHIV and HCs from 500FG, we used different samples compared with the previous studies (van der Heijden et al., 2021;Van de Wijer et al., 2021), so we conducted a reanalysis. A linear regression model was used, with correction for sex and age, using FDR < 0.05 as the significant threshold. Differentially abundant cytokine production and eight kinds of anti-inflammatory cytokine production (IL-10 and IL-1Ra) were included in the subsequent analysis. Furthermore, in the HIV cohort, the plasma levels of intestinal fatty acid-binding protein (iFABP), a marker of microbial translocation, and monocyte activation markers (sCD14 and sCD163) were measured using ELISA (Duoset or Quantikine, R&D Systems).

Microbial associations to HIVrelated variables, cytokine production capacity, intestinal damage, and monocyte activation
For association analysis between gut microbiome and HIV-related variables, cytokine production capacity, intestinal damage, and monocyte activation, we included bacterial alpha diversity (Shannon index), beta diversity, Prevotella-to-Bacteroides (P/B) ratio, DI and FI scores, as well as 76 species and 163 pathways that were significantly different between PLHIV and HCs from DMP cohort. Before Spearman correlation analysis, we first inverse-rank-transformed the data to follow a standard normal distribution and then adjusted the bacteria-phenotype associations for confounding factors (age and sex) and a technical variable (read counts). Significant associations were defined at FDR < 0.1 level. In additional analyses, we further adjusted for sex behavior factors [sexual orientation (SO) and the number of sexual partners in the previous year (Num-P)].

Differences in microbial composition and function
The present study included 143 PLHIV and 190 healthy individuals with matched age and sex from the DMP cohort (hereafter referred to as matched HCs). Participants' baseline characteristics are shown in Table 1. Included PLHIV were on long-term ART (median of 6.35 years) and were virally suppressed (plasma HIV-RNA < 200 copies/mL). However, when compared with the matched HCs, the gut microbial composition of the PLHIV still showed a significant decrease in alpha diversity (species-level Shannon index, Wilcoxon rank sum test, p = 2.4 × 10 −5 ; Figure 1A) and a distinct composition that was reflected by a significant difference in beta diversity (PERMANOVA, p < 1.0 × 10 −3 , R 2 = 0.07; Figure 1B), even when adjusting for BMI, read counts, and smoking status. These observations are consistent with findings of previous studies (Nowak et al., 2015;Zilberman-Schapira et al., 2016;Gootenberg et al., 2017). With the aid of metagenomics data, we also observed that PLHIV showed a lower diversity of bacterial metabolic pathways (pathway-level Shannon index, Wilcoxon rank sum test, p = 4.0 × 10 −4 ; Supplementary Figure 1) and different functional profiles (PERMANOVA, p < 1.0 × 10 −3 , R 2 = 0.04; Figure 1C). To explore the bacterial taxa and pathways that showed significantly different abundance in PLHIV, we confined the differential abundance analysis to the 123 common species and 331 metabolic pathways present in ≥ 20% of samples in at least one cohort. A linear regression model with correction for read counts, BMI, and smoking status revealed 76 (62.6%) species and 163 pathways (49.2%) with differential abundances between PLHIV and HCs (FDR < 0.05; Supplementary Tables 1, 2). Notably, sexual behavior is a known factor that influences the gut microbiome (Tuddenham et al., 2020), and the PLHIV was enriched for c (MSM) (n = 96, 67%). We further included a control group of 53 men from the DMP cohort who reported to have male partners (DMP-MSM) during sample collection time (Table 1). Of the 76 differentially abundant species and 163 pathways, 51 species (89.5%) and 102 pathways (62.6%) were replicated at FDR < 0.05 level (Supplementary Tables 3, 4, Supplementary Figures 2, 3), suggesting that most of our reported differences between PLHIV and healthy individuals were not cofounded by sex behavior.

Differentially abundant species
Fifty-seven of the 76 differentially abundant species showed enrichment in PLHIV. Increased abundances of Prevotella and Prevotellaceae in PLHIV were widely observed by previous studies PLHIV show a distinct gut microbiome composition and function compared with HCs. (A) Comparison of microbial alpha diversity between the PLHIV cohort and HCs from the DMP cohort. Y-axis refers to the Shannon index at the species level. (B, C) Beta diversity based on Bray-Curtis distance of species and pathway abundance is shown in a principal coordinates analysis (PcoA) plot with centroids for PLHIV and HCs. The coordinates of the centroids are set as the mean value of the principal components for each cohort. (D) Heatmap depicting the relative abundance of the top 30 species that differed significantly between PLHIV and HCs. Data are CLR-transformed and then inverse-rank-transformed to follow a normal distribution. Differentially abundant species were selected using a linear regression model with correction for BMI, read counts, and smoking status. (E) Network of bacterial pathways that were significantly different between PLHIV and HCs. Rectangular nodes represent super-pathways, with the colors showing enrichment (pink) or depletion (blue) in PLHIV. Circular nodes show sub-pathways belonging to the super-pathways. The color of circles shows the log2 value of fold change between the relative abundance of pathway in PLHIV and HCs, where a gradient is applied depending on fold change. Circle size indicates p-value. Lines connect each pathway to its respective super-pathway. Only super-pathways including two or more pathways and sub-pathways with FDR < 0.05 are shown. Differentially abundant pathways were selected using a linear regression model with correction for BMI, read counts, and smoking status. (F-H) Density curves of the P/B ratio, DI score, and FI score for the three cohorts depicting the different distribution of these bacterial signatures in these cohorts. Significance was tested using Dunn's test.

Differentially abundant pathways
At the metabolic pathway level, we observed differential abundances in several amino acid biosynthesis pathways, including enriched L-tryptophan biosynthesis (PWY-6629) and depleted L-ornithine and L-citrulline biosynthesis pathways (ARGININE-SYN4-PWY and CITRULBIO-PWY) in PLHIV (Supplementary Table 2). Importantly, both tryptophan and citrulline play critical roles in inflammation (Murray, 2003;Papadia et al., 2010), whereas ornithine can later be turned into nitric oxide (NO), which is important for vascular function (Baumgartner et al., 2000;Dirajlal-Fargo et al., 2017). In addition, the reductive tricarboxylic acid (TCA) cycle I (P23-PWY) was enriched in PLHIV. The reductive TCA cycle is a carbon dioxide-fixation pathway significant for the production of organic molecules for the biosynthesis of sugars, lipids, amino acids, and pyrimidines (Smith and Morowitz, 2004). Moreover, after clustering bacterial pathways to the same biochemical processes, we found that PLHIV showed lower abundances of bacterial pathways involved in fatty acid, lipid, and carbohydrate biosynthesis but higher abundances of bacterial pathways involved in amino acids biosynthesis and TCA cycle ( Figure 1E).

Dysbiosis index
Together, our data show dysbiosis in both gut microbial composition and metabolic function in PLHIV. Previous studies have suggested the P/B ratio as a landmark parameter for PLHIV (Dillon et al., 2014;Ling et al., 2016), and our finding of a significantly higher P/B ratio in PLHIV compared with HCs confirms these observations (Dunn's test, p = 9.4 × 10 −13 ; Figure 1F, Supplementary Table 5, Supplementary Figure 5A).
We also constructed a DI based on the differentially abundant species by calculating the log2 ratio of the geometric mean of PLHIV-enriched species (57 species) to PLHIV-depleted species (19 species). This DI score was significantly higher in PLHIV than in HCs (p = 9.9 × 10 −33 ; Figure 1G, Supplementary Table 5). We then sought to validate the DI score in an independent cohort with similar metagenomic data and the same DNA isolation method. Whereas no such data were available for an independent cohort of PLHIV, data whereas available for a separate healthy cohort: 500FG (Table 1). We found that the DI score of 500FG was not different from the DMP controls (p = 0.29; Supplementary Figure 5B) but was significantly lower than that of our cohort of PLHIV (p = 9.6 × 10 −37 ). We also constructed a FI score using the log2 ratio of the geometric mean of PLHIV-enriched bacterial pathways (87 pathways) to PLHIV-depleted bacterial pathways (76 pathways) ( Figure 1H). This FI score was also significantly higher in PLHIV than in DMP HCs (p = 1.1 × 10 −32 ; Supplementary Table 5, Supplementary Figure 5C), and this was supported by the observation that the 500FG FI was also significantly lower than that of PLHIV (p = 6.7 × 10 −18 ).

HIV-associated gut dysbiosis associates with clinical phenotypes
We conducted a systematic association analysis between HIVrelated variables and microbial alpha diversity, beta diversity, P/B ratio, DI score, and FI score, correcting for sex, age, and read counts (Supplementary Figures 6, 7, Supplementary Table 6). HIV clinical parameters were also taken into account, including the time between HIV diagnosis and inclusion in the study or cART initiation, CD4 + T-cell counts (nadir and latest, recovery after cART), plasma viral loads, and HIV-1 reservoir measurements in circulating CD4 + T cells, including the CD4 + T-cell-associated HIV-1 DNA (CA-HIV-DNA) and CD4 + T-cell-associated HIV-1 RNA (CA-HIV-RNA) levels. Sexual behavior was also considered part of the HIV clinical phenotype, including Num-P and SO, e.g., MSM and receptive anal intercourse (RAI).
We did not observe any significant association between the Shannon index of bacterial species and any HIV clinical phenotypes (Supplementary Table 7). However, six parameters were associated with bacterial beta diversity, four with P/B ratio, three with DI score, and four with FI score at FDR < 0.1 level (Figure 2A, Supplementary  Table 7). Notably, although most of the reported microbial differences between PLHIV and HCs were not driven by sexual behavior (Supplementary Tables 3, 4), within PLHIV, sex behavior was still among the strongest associated factors. For example, SO was the top factor positively associated with the P/B ratio (Spearman correlation, p = 5.7 × 10 −12 ) and DI score (p = 6.4 × 10 −12 ), followed by Num-P and RAI (Supplementary Table 7). Moreover, MSM and RAI, as well as larger Num-P, were associated with increased P/B ratio, DI score, and FI score. We also observed associations for other HIV clinical phenotypes, such as an association between beta diversity and the HIV reservoir parameters CA-HIV-DNA and CA-HIV-RNA levels (p = 6.0 × 10 −3 and 8.0 × 10 −3 , respectively) and an association between P/B ratio and CD4 recovery relative rate (p = 0.02). However, these associations were largely dependent on sexual behavior, as none remained significant after correcting for SO and Num-P (Supplementary Table 7, Supplementary Figure 8). We further assessed whether HIV clinical parameters were associated with individual species and metabolic pathways. After correcting for age, sex, read counts, and sexual behavior, no significant associations were detected for species or pathways, although we did observe suggestive associations for Firmicutes bacterium CAG 95 and Prevotella sp. CAG 5226 with HIV reservoir measurements at a normal significance level ( Figures 2B-D, Supplementary  Table 8). At the metabolic pathway level, HIV duration and cART duration tended to be the strongest factors in addition to sexual behavior factors linked with metabolic pathways (Supplementary Table 9). For example, the pathways of stearate, octanoyl, and oleate biosynthesis showed a positive correlation with HIV duration, whereas the top result for cART duration was a negative association with the polyamine biosynthesis pathway (Supplementary Table 9). We did not observe any significant associations between the gut microbiome and cART regimens after correcting for sexual behavior (Supplementary Tables 10, 11).

Prevotella copri and Bacteroides vulgatus associate with cytokine production capacity
Despite the long-term cART, the immune responses of PLHIV are known to be different than those of HCs ( Van de Wijer et al., 2021). For instance, our previous study reported a significant increase in the production of proinflammatory cytokines in PLHIV (van der Heijden et al., 2021). In the present study, we newly compared the cytokine production data of 143 PLHIV (van der Heijden et al., 2021), for whom we collected microbiome data for this study, to 173 HCs from the 500FG cohort for whom the cytokine production capacity was measured previously (Schirmer et al., 2016) (Supplementary Table 12). We identified 24 cytokine abundances that were significantly different between PLHIV and healthy population ( Figure 3A, Supplementary Table 13). PLHIV showed a significant increase of IL-1b, IL-6, and TNF production upon stimulation with Pam3Cys (TLR2 ligand), LPS (TLR4 ligand), and C. albicans hyphae but decreased IFN-g production upon stimulation with S. aureus and C. albicans hyphae. Moreover, cytokine production capacities were also related to some HIVrelated phenotypes. MSM status was associated with increased Association between HIV-associated gut dysbiosis and HIV-related variables. (A) Heatmap depicting the associations between HIV-associated bacterial signature (Shannon index, beta diversity, P/B ratio, DI score, and FI score) and HIV-related variables, using the Spearman correlation test, with correction for age, sex, and read counts. Box color indicates Spearman correlation rho if the association p-value < 0.05. The associations at FDR < 0.1 level are highlighted with *. R 2 is calculated using PERMANOVA based on Bray-Curtis distance of species and then multiplied by 10 to rescale. (B-D) Associations between individual species and HIV reservoir level, with correction for age, sex, read counts, and sexual behavior. proinflammatory cytokine responses (e.g., Pam3Cys-induced TNF and S. aureus-induced IL-22 production), whereas higher Num-P was linked with decreased anti-inflammatory cytokine responses (e.g., Pam3Cys-and LPS-induced IL-10 production; Supplementary  Figure 9, Supplementary Table 14).
Interestingly, we observed significant associations of cytokine production capacity with microbial alpha and beta diversity and P/B ratio, as well as four associations with DI and FI scores and individual species (Supplementary Tables 15, 16). No associations were observed with metabolic pathway abundance (Supplementary Table 17). In particular, IL-10 production upon stimulation with Pam3Cys or LPS was negatively associated with P/B ratio and FI score but positively associated with bacterial Shannon index (linear regression, p < 0.05; Figure 3B, Supplementary Table 15), after correction for age, sex, and sexual behavior. At the individual species level, the top association was between Pam3Cys-induced IL-10 production and the relative abundance of P. copri (rho = −0.37, p = 9.1 × 10 −6 ; Supplementary Table 16), after correcting for age, sex, and read counts. After further adjustment for SO and Num-P, the association between P. copri and IL-10 production remained significant (rho = −0.34, p = 1.9 × 10 −4 ; Figure 3C, Supplementary Table 16). In addition, we also detected a positive significant association between B. vulgatus and Pam3Cysinduced IL-1b production (rho = 0.33, p = 7.4 × 10 −5 ), and this remained significant after controlling for age, sex, read counts, and sexual behavior (rho = 0.32, p = 3.7 × 10 −4 ; Figure 3D). Notably, this association was not significant in 500FG (rho = −7.8 × 10 −3 , p = 0.93; Figure 3E), showing a significant heterogeneity effect (Cochran's Qtest, p = 0.003; Supplementary Table 16).

Prevotella copri strains in PLHIV are genetically different
In microbial species, strain-level genomic makeup is critical in determining their functional properties within human bodies (Scholz et al., 2016). We therefore wondered whether PLHIV harbor different strains of P. copri and B. vulgatus, thereby affecting cytokine production capacity. To examine this, we performed an analysis on the basis of the presence or absence of their gene repertoire using PanPhlan3 (Beghini et al., 2021). As B. Association between gut microbiome and inflammatory cytokine production. (A) Heatmap showing ex vivo cytokine production enriched (red) and depleted (blue) in PLHIV as compared with HCs from 500FG. A linear regression model (age and sex included as covariates) was used to calculate the P-values. Box color indicates the Spearman correlation rho for those associations that were significant at p-value < 0.05 level. The associations at FDR < 0.1 level were highlighted with *. White box indicates associations with p-value > 0.05, and gray box indicates no measurement.
(B) Heatmap showing Spearman correlation rho between cytokine production and HIV-associated bacterial signature (Shannon index, beta diversity, P/B ratio, and DI and FI scores), with correction for age, sex, read counts, and sexual behavior. White box indicates P > 0.05. The row labels indicating the cytokine produced and the name of stimulation. (C, D) Association between relative abundance of species and cytokine production in PLHIV: (C) Prevotella copri with Pam3Cys-induced IL-10 production and (D) Bacteroides vulgatus with Pam3Cys-induced IL-1b production. (E) Association between relative abundance of Bacteroides vulgatus and Pam3Cys-induced IL-1b production in HCs from 500FG. The relative abundance of the species was CLR-and inverse-rank-transformed. Association was corrected for age, sex, read counts, and sexual behavior using linear regression. Pam3Cys, synthetic TLR2 ligand; LPS.1ng, PBMC cytokine production in response to LPS (1 ng/mL) stimulation; LPS.100ng, PBMC cytokine production in response to LPS (100 ng/mL) stimulation; RPMI, serum-free medium; Poly IC, TLR3 ligand.
vulgatus did not show any distinct clusters according to genetic content, we did not study it further. We generated P. copri genetic repertoire profiles for 282 samples from the PLHIV cohort (n = 102), DMP cohort (n = 111), and 500FG cohort (n = 69; Supplementary Table 18). Hierarchical clustering analysis based on the Jaccard distance of the presence or absence of gene families revealed three distinct clusters at a distance cutoff of 0.3 ( Figure 4A). Cluster 1 was relatively small (n = 13) and was therefore excluded in the following analysis. In contrast, clusters 2 and 3 were larger (n = 150 and 119, respectively; Supplementary Table 19). Interestingly, cluster 2 was significantly enriched in healthy individuals, and cluster 3 was enriched in PLHIV (Fisher exact test, p = 3.7 × 10 −22 ; Supplementary Table 19). We hereafter refer to these two distinct P. copri strains as the "control-related strain" and "PLHIVrelated strain," respectively. We further assessed the difference in the gene content of the two strains. A total of 2,821 of the 4,998 gene families were differentially abundant (logistic regression, FDR < 0.05; Figure 4A, Supplementary Table 20), and seven protein families (Pfam) (Finn et al., 2010) were enriched in the PLHIVrelated strain (hypergeometric test, FDR < 0.05; Figure 4B, Supplementary  Table 22). Interestingly, the frequency of different P. copri strains was also associated with sex behavior, and PLHIV-related strain showed an enrichment in PLHIV-MSM than that in PLHIV-MSW (Fisher exact test, p = 0.01). However, after controlling for the behavior of

Prevotella copri strain profile is associated with cytokine production capacity
We further compared the associations of different P. copri strains with cytokine production capacity. The control-related strain showed stronger associations with cytokine production than PLHIV-related strain, especially for Pam3Cys-induced IL-6 and IL-10 production in PLHIV ( Figures 5A, B, Supplementary Table 23). However, in the 500FG cohort, the control-related strain showed an opposite association with IL-6 production (Supplementary Figure 11A, Supplementary Table 24). In contrast, PLHIV-related strain also showed a negative association with IL-10 production in response to Candida albicans in PLHIV; however, the association did not show heterogeneity against control-related strain (Supplementary Figure 11B, Supplementary Table 23).
We also checked the association between P. copri strains and HIV-related parameters, as well as with monocyte activation markers (sCD14 and sCD163) and microbial translocation marker (IFABP). The control-related strain showed a positive association with CD4 + T-cell counts (CD4 counts) and CD4 + Tcell absolute recovery level after cART (CD4-recovery-abs) ( Figures 5C, D, Supplementary Table 25). These strain-CD4 associations did not show heterogeneity between the control-and PLHIV-related strains (Supplementary Table 25). However, we did not observe any significant associations between P. copri strains and monocyte activation and microbial translocation markers (Supplementary Table 26).
To explore the potential mechanism behind the distinct immune functions of the two P. copri strains, we searched the IEDB and found five epitope peptides derived from different P. copri proteins (Supplementary Table 27). These peptides can be presented by antigen-presenting cells and induce IFN-g production by T cells and antibody secretion by B cells. We then compared the abundances of these five peptides between the two P. copri strains from PLHIV and 500FG cohort (Methods). Interestingly, only one peptide from P. copri, glutamate 5-kinase protein, was found to be significantly enriched in the samples with the control-related strain as compared with the samples with the PLHIV-related strain (linear regression, p = 6.1 × 10 −5 ; Figure 5E, Supplementary Table 28), suggesting that it is potentially contributing to the immune function of the control-related strain.

Discussion
We performed metagenomics-based microbiome associations of the gut microbiome with ex vivo cytokine production capacity in PLHIV. First, we observed a remarkable microbial dysbiosis in PLHIV. The microbial diversity was significantly lower, and the abundance levels of 76 species and 163 metabolic pathways were significantly different from HCs. Sex behavior was identified as an important factor that could influence the gut microbiome, but sex behavior could not explain all the PLHIV-associated microbial differences. Second, PLHIV-associated microbial dysbiosis was associated with clinical parameters and cytokine production capacity in PLHIV, including HIV reservoir, CD4 recovery relative rate, and production capacity of IL-10 and IL-1b. Interestingly, P. copri species in PLHIV and HCs showed distinct genetic differences, suggesting that two different strains were present in PLHIV and HCs, respectively. We also found that the control-related P. copri strain showed a stronger negative association with IL-6 and IL-10 production and a positive association with CD4 counts.
In the comparison of gut microbial composition between PLHIV and HCs, we observed a lower alpha diversity and higher P/B ratio in PLHIV ( Figures 1A, F), as well as a depletion of Alistipes species ( Figure 1D) (Crakes and Jiang, 2019;Vujkovic-Cvijin and Somsouk, 2019). In contrast to previous studies, we found some SCFA-producing species and beneficial species to be enriched in PLHIV in our study, including Acidaminococcus fermentans (Chang et al., 2010) and Faecalibacterium prausnitzii (Parada Venegas et al., 2019), as well as the equol-forming species Slackia isoflavoniconvertens (Schröder et al., 2013) and the antitumorigenic species Holdemanella biformis (Zagato et al., 2020). This inconsistency may be due to the fact that the PLHIV in our cohort had a longer duration of ART intake (median of 6.4 years) compared with previous studies, which is supported by earlier evidence that ART can partially restore the gut microbial composition (Lozupone et al., 2013). Functionally, PLHIV in our study showed an increased microbial capacity of L-tryptophan biosynthesis and a decreased de novo biosynthesis of ornithine from 2-oxoglutarate, functions that have been related to inflammation and vascular function (Murray, 2003;Dirajlal-Fargo et al., 2017), respectively. This observation supports the idea that those functional changes in the gut microbiome in PLHIV contribute to the persistent inflammation seen in these individuals. Importantly, we also showed that sex orientation can influence the gut microbiome but the microbial dysbiosis between PLHIV and HCs is not driven by sex orientation only.
Infected long-lived memory CD4 + T cells make up the majority of cell types constituting the HIV reservoir (Churchill et al., 2016). To the best of our knowledge, no other studies have ever reported associations between the microbiome and the HIV reservoir. We analyzed HIV-1 DNA and RNA levels in isolated circulating CD4 + T cells, which reflect the size of the HIV reservoir (Rutsaert et al., 2019). We observed that the CD4 + T-cell counts were significantly associated with bacterial beta diversity (Figure 2A). We also found Prevotella copri strains with different genetic content. (A) Heatmap showing the gene family profiles of P. copri strains in samples from the three cohorts (282 samples total: 102 PLHIV, 111 DMP, and 69 500FG). Each column represents a sample. Each row represents the presence or absence of a gene family. The clustered tree above the heatmap shows the three clusters of P. copri strains. Most samples from PLHIV were binned together into the PLHIV-related strain (right), but 16 PLHIV samples (middle) showed different profiles and were binned into the control-related strain. that Firmicutes bacterium CAG 95 negatively correlated with CA-HIV-DNA levels and that the relation between R. lactatiformans and CA-HIV-RNA level was positive, whereas Prevotella species showed negative associations with CA-HIV-DNA and CA-HIV-RNA ( Figures 2B-D). On one hand, how these bacterial species modulate the HIV reservoir remains speculative. R. lactatiformans has been linked with increased colonic IFN-g + T cells and immune activation (Frankel et al., 2019), whereas a decreased abundance of Firmicutes bacterium CAG 95 was found in subjects with hepatic steatosis (Zeybel et al., 2022). A previous study identified several immunogenic human leukocyte antigen (HLA) -DR isotypepresented P. copri peptides with the ability to induce T cells to produce the anti-viral cytokine IFN-g (Kak et al., 2018;Pianta et al., 2021). In our study, the P. copri HC-related strain was associated with decreased levels of IL-6 production, which can facilitate HIV-1 replication (Connolly et al., 2005). On the other hand, HIV infection leads to a severe gastrointestinal (GI) CD4 + T cells loss, which may affect immunological tolerance to specific gut microbes and thereby result in bacterial changes (Crakes and Jiang, 2019). This is supported by previous evidence that PLHIV with CD4 counts < 350 cells/mm 3 for 2 years of cART showed an enrichment of unclassified Subdoligranulum species and Coprococcus comes, compared with those with CD4 counts ≥ 350 cells/mm 3 (Lu et al., 2018). Importantly, the relative abundances of unclassified Subdoligranulum species and C. comes were positively correlated with CD8 + HLA-DR + T-cell count and CD8 + HLA-DR + /CD8 + percentage in PLHIV (Lu et al., 2018). However, we did not observe differences in bacterial composition between PLHIV with CD4 counts ≥ 500 cells/mm 3 and PLHIV with CD4 counts < 500 cells/mm 3 . Moreover, the HIV reservoir was also correlated with CD4 counts, HIV infection, and treatment duration (Supplementary Figure 7). Therefore, the causal direction between the gut microbiome and HIV reservoir remains inconclusive.
The gut microbiome may influence the host's immune response via their immunoregulatory metabolites or peptides. Moreover, gut microbes can also pass the intestinal barrier and translocate into the systematic circulation, thereby triggering an immune response (Brenchley et al., 2004;Mattapallil et al., 2005;Mudd and Brenchley, 2016). Indeed, microbial translocation has previously been observed in PLHIV (Nganou-Makamdop et al., 2021). It was Prevotella copri strains with different immune functions. (A, B) Distinct associations between control-related strain and PLHIV-related strain with IL-6 and IL-10 production capacity in PLHIV, using linear regression. Association is corrected for age, sex, read counts, and sexual behavior. Blue dots represent the PLHIV with the control-related strain. Dark pink indicates PLHIV with the PLHIV-related strain. Bacterial relative abundance is CLR-and inverse-rank-transformed. (C, D) Distinct associations between relative abundance of control-related strain and PLHIV-related strain with CD4 counts (C) and CD4-recovery-abs level (D). (E) Violin plot of the P. copri epitope peptide (from glutamate 5 kinase) level between the two kinds of strains in PLHIV and HCs from the 500FG cohort.
observed that translocated Serratia genera in blood drive innate and Th17 cytokine responses and associate with reduced gut barrier integrity (Nganou-Makamdop et al., 2021). We compared the genera and species enriched in the PLHIV in our study with previously reported translocated microbes in PLHIV, such as Serratia, Acidovorax, Sphingobium, and Burkholderia genera. However, abundances of these blood genera were either not detected in the fecal samples or not associated with HIV-related phenotypes and immune response in our HIV cohort. Specifically, we did not detect Serratia genera in our cohort, potentially because Serratia genera tend to colonize the respiratory and urinary tracts rather than the GI tract (Mahlen, 2011). We further detected two distinct P. copri strains between PLHIV and HCs. Interestingly, the control-related strain showed stronger negative associations with cytokine production capacity than the PLHIV-related strain in PLHIV, particularly for IL-6 and IL-10 production ( Figures 5A, B). A large proportion of PLHIV seem to have lost this control-related P. copri strain and show higher levels of IL-6 and IL-10 production. IL-6 and IL-10 both play a critical role in HIV pathogenesis and are used as markers of disease progression (Freeman et al., 2016). For example, a higher level of plasma IL-6 has been linked with increased morbidity and mortality and the failure to reconstitute CD4 counts (Rose-John et al., 2017). Higher PBMC production of IL-10 is linked to low CD4 counts (Vujkovic-Cvijin et al., 2013), and plasma IL-10 levels positively correlate with viral load (Chehimi et al., 1996). The association between P. copri strain and cytokine production may be due to specific P. copri immunogenic epitopes. Some P. copri immunogenic epitope peptides have been shown to induce IFN-g production of T cells (Pianta et al., 2021), and the control-related, immune-relevant P. copri strain identified in the current study indeed showed an enrichment of immunogenic epitopes ( Figure 5E). Furthermore, the altered P. copri strain may also show different microbial translation and different local interactions with the intestinal immune system. P. copri can even translocate to the systematic circulation and cause infection (Posteraro et al., 2019). Immunoglobulin G (IgG) specific to Prevotella species including Prevotella intermedia and Prevotella gingivalis has been found in circulation in patients with rheumatoid arthritis (Mikuls et al., 2012). However, more studies are required to understand the functionality differences between these two strains we identified.
We acknowledge several limitations of the current study. Although it is, to date, a large metagenomics-based study in PLHIV, the sample size is still relatively small. The analysis power can be small after correcting the number of tests performed. Another limitation of this study is the potential batch effects among the different cohorts due to non-biological factors such as technical differences. We also did not have an independent cohort of PLHIV to validate our findings, but we could replicate the DI and FI scores in the 500FG cohort (built in the same medical center and used the same sample collection method as the PLHIV cohort). In addition, we focus on the PLHIV-associated bacterial changes observed in PLHIV who are mostly MSM. Hence, our conclusions may not be generalizable to all PLHIV. However, because of the majority of PLHIV in Europe is MSM, it is still valuable to explore this PLHIV-associated bacterial changes that combined the effects of HIV infection and sexual behavior together. Another potential confounding effect would be related to disease comorbidity. PLHIV show a higher risk for various diseases, particularly cardiometabolic diseases. In our PLHIV cohort, 18.9% of patients have hypertension, 24.5% have mental problems, and 7% have diabetes. However, HCs were only matched on age, sex, and BMI, so we cannot exclude that the observed gut dysbiosis was due to disease comorbidities. We also acknowledge that IL-10 production data were not available for the HCs, so we could not compare changes in IL-10 production capacity between PLHIV and HCs. Furthermore, no blood microbiome data are available for our cohort; thus, we cannot directly assess microbial translocation in PLHIV and to what extent gut microbes can be translocated into the systematic circulation and thereby affect the host's immune response. Finally, ethnicity is also an important factor associated with gut microbial composition in PLHIV (Noguera-Julian et al., 2016). All PLHIV in our cohort have European ancestry. Hence, in the future study, we suggested collecting samples from across different ethnicities and controlling for its confounding effect. For the technical side, DNA isolation methods can also bias microbial profiling, making microbial data not comparable across different studies. The current study employed the FSK protocol that included a heating step in combination with the enzymatic lysis. Although previous studies have described that this cell lysis method can favor bacterial cell lysis by denaturizing the membrane proteins, this protocol seemed to yield a lower DNA concentration and a microbial community with a lower diversity, compared with DNA isolation method using a bead-beating step as the mechanical lysis (Fernańdez-Pato et al., 2022).
In conclusion, we observed differential microbial composition and function on the species and strain levels in long-term treated PLHIV. This HIV-associated bacterial signature was linked with HIV reservoir parameters and with PBMC production capacity of IL-1b and IL-10. A large fraction of the PLHIV have lost the control-related P. copri strain that was associated with IL-6 and IL-10 production capacity and with CD4 counts and CD4-recoveryabs. The loss of this control-related strain may contribute to a higher level of IL-6 and IL-10 production in PLHIV and to later immune activation and dysfunction. Our observation has provided deeper insight into the critical role of the gut microbiome during HIV infection. The immune-related P. copri strain may be used as part of treatment for chronic inflammation, particularly cytokine imbalance; however, more studies are warranted.

Data availability statement
Raw metagenomic sequencing data of HIV cohort are publicly available from SRA under project number PRJNA820547. Raw metagenomic sequencing data of Dutch Microbiome Project (DMP) and 500 Functional Genomics (500FG) cohort are publicly available from the European Genome-Phenome Archive via accession number EGAS00001005027 and from NCBI Short Read Archive (SRA) via accession number PRJNA942468, respectively. The code used for statistical analysis is available via GitHub: https://github.com/White-Shinobi/HIV-andgut-microbiome.

Ethics statement
The 200HIV study was approved by the Medical Ethical ReviewCommittee region Arnhem-Nijmegen (CMO2012-550). For theDutch Microbiome Project (DMP) cohort, the Lifelines study wasapproved by the medical ethical committee from the UniversityMedical Center Groningen (METc number: 2017/152). The 500Functional Genomics (500FG) study was approved by the EthicalCommitt ee of Radboud University Nijmegen (NL42561.091.12,2012/550). All informed consents were collected for allparticipants. Experiments were conducted in accordance with theprinciples of the Declaration of Helsinki.

Author contributions
JF and AV conceptualized and managed the study. WH, LV, QM, LJ, RW, AZ, and MN contributed to data generation. YZ, SA-S, DW, and RG analyzed the data. YZ, JF, and AV drafted the manuscript. SA-S, NV, DW, VM, WH, RG, RW, AZ, LV, QM, LJ, and MN reviewed and edited the manuscript. All authors contributed to the article and approved the submitted version. Organization for Scientific Research. RKW is supported by the Seerave Foundation and the Dutch Digestive Foundation (16-14). YZ is supported by a joint fellowship from the University Medical C e n t r e G r o n i n g e n a n d C h i n a S c h o l a r s h i p C o u n c i l (CSC202006170040). Moreover, AZ, LJ, MN, and JF are also supported by Dutch Heart Foundation IN-CONTROL (CVON2018-27).