This cross-sectional study included HIV-1-ART patients who were recruited from the outpatient clinic at Karolinska University Hospital, Stockholm, Sweden, and HIV-1-negative controls. The inclusion criteria for the HIV-1-ART participants were age >18 years (median: 45 years; range: 38–62 years), and HIV-1 positive for at least 6 months. The exclusion criteria were ongoing HIV-1-related complications or antibiotic treatment during the previous 3 months. The patients received ART for a median of 7.7 years (range: 4.4–20.8 years). The HIV-1-negative controls were healthy individuals (median: 32 years; range: 24–51 years) who did not receive antibiotic treatment during the last 3 months.

This study was approved by the Regional Ethics Committee, Stockholm (2009/1485-31, 2013/1944-31/4, 2014/920-3).

A sterile tube for fecal sampling without preservation media was used when participants were able to donate feces at the clinic as previously described (Vesterbacka et al., 2017). Samples were frozen and stored at −80°C within 24 h. Stool collection tube with DNA stabilizer (Stratec Biomedical) was used for participants who collected the feces at home. The samples were delivered to the outpatient clinic by the participants or instantly sent by post and stored at −80°C according to the manufacturer’s instructions.

Clinical data of the HIV-1-ART participants were collected from the InfCareHIV database, including gender, age, current, and nadir CD4 + T-cell counts, ART regimen at sampling date, total ART duration (years), and transmission mode (Table 1).

First, we did literature review and selected and evaluated two widely used protocols, i.e., QIAamp PowerFecal Pro DNA Isolation kit (QP) (Qiagen, Germany) and the standardized International Human Microbiota Standards (IHMS) Protocol Q recommended by the International Human Microbiome Consortium¹ (Costea et al., 2017), which performed better than other methods in terms of both DNA quality and microbial diversity in human fecal microbiota studies (Angebault et al., 2018; Lim et al., 2018; Fiedorova et al., 2019), while never having been compared before. A subset of fecal samples was randomly selected from cases and controls and then subjected to the two protocols according to instructions with minor modifications. The quantity and quality of gDNA extracted by the two protocols were measured using NanoDrop^TM One (Thermo Scientific, United States) and Agilent 2200 TapeStation (Agilent Technologies, United States) including gDNA concentration, purity, and integrity. The protocol yielding higher gDNA quantity and better quality was then used to extract gDNA of the remaining samples. The gDNA samples were stored at −80°C until library preparation and sequencing.

Sequencing libraries were prepared with the Nextera DNA Flex kit (Illumina, Inc.) following the manufacturer’s instructions. One paired-end (PE) library with insert size of approximately 320 bp was constructed for each sample. Libraries were normalized with Qubit assay; the pooled library was then sequenced on one lane on NovaSeq6000 platform (NovaSeq Control Software 1.6.0/RTA v3.4.4) with a 2 × 150 setup using “NovaSeqXp” workflow in “S1” mode flowcell. The Bcl to FastQ conversion was performed using bcl2fastq_v2.20.0.422 from the CASAVA software suite. The quality scale used is Sanger/phred33/Illumina 1.8 +.

A total of 146 GB raw sequencing data were generated for 24 fecal samples, approximately 6.1 GB per sample was obtained. The raw sequencing data were preprocessed using our in-house bioinformatics pipelines (Figure 1). Briefly, the adapter and low-quality reads (a quality score of less than Q30) were removed using Trim galore (v0.6.4).² After the quality trimming, Bowtie2 (v2.3.5.1) (Langmead and Salzberg, 2012) was used in combination with SAMtools (v1.19) (Li et al., 2009) and BEDtools (v2.29.2) (Quinlan and Hall, 2010) to identify and remove human DNA sequences. The human unmapped reads were then used for downstream analysis. After filtering, an average of 41 million 150-bp PE reads per sample was obtained, the adapter and low-quality reads accounted for 1.32% of total reads, and the average rate of human DNA contamination in all samples was 0.58% (Supplementary Table S1).

The taxonomic assignment and abundance estimation were performed with MetaPhlAn 2.0 (Truong et al., 2015) using default parameters. MetaPhlAn 2.0 relies on ∼1 million unique clade-specific marker genes identified from ∼17,000 reference genomes (∼13,500 bacterial and archaeal, ∼3,500 viral, and ∼110 eukaryotic genomes) (Truong et al., 2015). This computational tool provides the relative abundances of each microbial clade with species-level resolution. The filtered reads were mapped to all reference genome sequences to determine the presence and abundance of different taxonomic groups present in human fecal samples. Given that this study focused on bacterial communities, downstream analysis was performed on bacterial taxa.

Filtered reads were mapped against the integrated gene catalog (IGC) (Li et al., 2014) using the bwa software (Li and Durbin, 2009). Gene richness was measured as the total number of different genes present in the sample regardless of their abundance and length. A minimum of one filtered mapped read was set to consider the presence of a gene, as described previously (Guillen et al., 2019).

Metagenomics functional analysis was performed using the HMP Unified Metabolic Analysis Network 2 (HUMAnN2), which identifies the species profile from metagenomic shotgun sequencing data and aligns reads to their pan-genomes, performs translated search on unclassified reads, and quantifies gene families and pathways (Franzosa et al., 2018). By default, gene families were annotated using a comprehensive protein database UniRef90 (Suzek et al., 2015) and metabolic pathways using MetaCyc database (Caspi et al., 2016). The UniRef90 gene family abundance from HUMAnN2 was then regrouped to Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000) orthology (KO). To gain a comprehensive profile of functional pathways, the filtered reads were further mapped to the KEGG database, which contains more compounds while having less reactions and pathways than does the MetaCyc database (Altman et al., 2013).

To characterize the presence of bacterial virulence factors genes and AMR genes, we mapped filtered reads from each sample against the VFDB (Liu et al., 2019) and the CARD (Alcock et al., 2020), respectively, using Bowtie2 (v2.3.5.1). The mapped read segments were estimated using SAMtools idxstats. The copy number of each gene was estimated by dividing the total reads mapping to a gene divided by the gene’s length. The relative abundance of genes was calculated using R function make_relative in R package funrar (v1.4.1).

Rarefaction curves were generated using R function rarecurve (R package vegan v2.5.6) to assess the saturation of samples at different sequencing depth for recovery of bacterial species, bacterial virulence factors genes, and AMR genes.

Species richness and Shannon diversity were calculated using R function estimate_richness. Beta diversity was measured by Bray–Curtis, weighted UniFrac, and unweighted UniFrac distances using R package Phyloseq (v1.30.0) (McMurdie and Holmes, 2013). Samples were clustered according to their species composition using non-metric multidimensional scaling (NMDS) approach based on Bray–Curtis, weighted UniFrac, and unweighted UniFrac distances in Phyloseq (v1.30.0).

Due to the small sample size in this study, different methods were used to determine and verify the results obtained from the other method. Differences in abundance of bacterial species, pathways, bacterial virulence genes, and AMR genes between two groups were determined using Wilcoxon rank sum test with a significance level of <0.05. The linear discriminant analysis (LDA) effect size (LEfSe) algorithm (Segata et al., 2011) was used to identify specific bacterial taxa and metabolic pathways as taxonomic and functional biomarkers and to verify the Wilcoxon rank sum test. Kruskal–Wallis test was used to process the dataset with LEfSe alpha values set at 0.05. The threshold used to consider a discriminative feature for the logarithmic LDA score was set at >2. The biomarker discovery was performed at bacterial order, family, genus, and species levels and all functional levels. All the correlation analyses were performed using the Spearman’s rank correlation coefficient (library “psych”; function “corr.test”). Spearman’s correlation (rho) with p < 0.05 was considered statistically significant. Statistical correction for multiple testing was not applied because of the low sample size and the exploratory focus of this study.

The study subjects were categorized into case (HIV-1-ART individuals, n = 13) and control (HIV-1-negative controls, n = 11) groups. Cases were classified into groups based on median value (740 cells/mm³) of current CD4 + T-cell counts: ≥740 cells/mm³ (n = 7) and <740 cells/mm³) (n = 6).

Fecal samples were randomly selected from cases (n = 3) and controls (n = 7) to evaluate two gDNA extraction protocols. Using the same fecal material mass, IHMS Protocol Q produced significantly higher gDNA yield (p < 0.0001, median DNA concentration 17.26 times higher) than the QP kit (Figure 2 and Supplementary Table S2). The 260/280 absorbance ratio was used to assess the gDNA purity, within an acceptable range of 1.8–2.0 as good-quality DNA. gDNA from the IHMS Protocol Q had 260/280 ratio ranging from 1.89 to 1.96, whereas four out of 10 samples extracted from the QP kit had 260/280 ratio below 1.8. The gDNA integrity was assessed with DNA Integrity Number (DIN) ranging from 1 to 10, where 1 indicates highly degraded gDNA and 10 represents highly intact gDNA (Nguyet et al., 2014). No significant difference in DIN values was found between the two protocols (p = 0.82); however, one sample with lower 260/280 ratio of 1.55 from the QP kit had also appreciably low gDIN value of 2.8. Our study indicated that the IHMS Protocol Q yielded significantly higher gDNA quantity and more stable gDNA quality than the QP kit. The remaining 14 samples were therefore extracted by the IHMS Protocol Q. All samples (n = 24) extracted by the IHMS Protocol Q were further subjected to shotgun metagenome sequencing and downstream analysis.

Surprisingly, the average proportion of filtered reads assigned to known bacterial species accounted for 39.9% of total reads, less than 0.1% of reads were assigned to archaea and virus, whereas 60% of reads could not be mapped to any known reference genome of different taxonomic groups by MetaPhlAn 2.0 (Supplementary Figure 1a and Supplementary Table S3). Given that the focus of this study was on bacterial communities, the downstream analysis was performed on bacteria only. Intriguingly, the bacterial species-level richness rapidly reached a plateau for all 24 samples at less than 5 million reads (Supplementary Figure 1b). MetaPhlAn 2.0 identified 288 bacterial species in the 24 fecal samples (Supplementary Table S4), belonging to 109 genera, 53 families, 23 orders, 14 classes, and seven phyla (Figure 3A). Clostridiales and Bacteroidales were the most abundant orders (Supplementary Figure 2a). Ruminococcaceae and Bifidobacteriaceae were the most abundant families; Faecalibacterium and Bifidobacterium were the most abundant genera (Figures 3B,C). The most abundant bacterial species observed in all fecal samples was Faecalibacterium prausnitzii (average relative abundance 15.97% in case group vs. 20.16% in control group), followed by Bifidobacterium adolescentis (11.17% case vs. 7.11% control), Collinsella aerofaciens (5.96% case vs. 6.51% control), and Ruminococcus bromii (4.16% case vs. 6.10% control). Sixteen species showed average relative abundance above 1% (Figure 3D). However, no significant difference of these abundant species was found between HIV-1-ART cases and HIV-1-negative controls.

Despite the small sample size, three families, five genera, and 13 species biomarkers were identified between the case and control groups (Figure 3E and Supplementary Figures 2b,c) using LEfSe biomarker discovery tool. Among the 13 species biomarkers, seven were identified as biomarkers for the case group and six for the control group (Figure 3E). Wilcoxon rank sum test was consistent with LDA with two exceptions (Corynebacterium pseudogenitalium and Streptococcus anginosus) (Supplementary Table S5A).

A decreased tendency of microbial richness and Shannon diversity at species level was observed in HIV-1-ART cases compared to HIV-1-negative controls, though differences were not statistically significant (p > 0.05) (Figure 4A). Beta diversity assessed with Bray–Curtis, unweighted UniFrac, and weighted UniFrac dissimilarities was significantly higher in the case group compared to that in the control group (p = 0.0025, p = 0.0031, p = 0.047, respectively) (Figure 4B). NMDS with Bray–Curtis distance (Figure 4C) and UniFrac distances (data not shown) showed no significant separation between cases and controls, although samples from cases were more diversely distributed as compared to controls. Based on IGC analysis, no statistically significant difference was found in microbial gene richness between cases and controls.

To assess the potential functionality of gut microbiota in HIV-1 infection, microbial genes were annotated and analyzed using different databases and functional systems: KEGG and MetaCyc pathways, KO. We found that 94.6% of filtered reads were unmapped or unintegrated against the KEGG database. The most abundant KEGG pathways identified in this study were biosynthesis of ansamycins (ko01051), valine, leucine and isoleucine biosynthesis (ko00290), ribosome (ko03010), D-Glutamine, and D-glutamate metabolism (ko00471) (Supplementary Figure 3a). Similarly, 95.4% of filtered reads were unmapped or unintegrated against the MetaCyc pathway database; the top MetaCyc pathways were adenosine ribonucleotide de novo biosynthesis (PWY-7219) (MetaCyc pathways class: nucleoside and nucleotide biosynthesis), UMP biosynthesis (PWY-5686) (nucleoside and nucleotide biosynthesis), UDP-N-acetylmuramoyl-pentapeptide biosynthesis II (lysine-containing) (PWY-6386) (cell structure biosynthesis), L-isoleucine biosynthesis I (from threonine) (ILEUSYN-PWY) (amino acid biosynthesis), pyruvate fermentation to isobutanol (PWY-7111) (generation of precursor metabolites and energy), and L-valine biosynthesis (VALSYN-PWY) (amino acid biosynthesis) (Supplementary Figure 3b). No statistically significant difference of these main pathways was found between the case and control groups. We found that one MetaCyc pathway, Bifidobacterium shunt (P124-PWY), also known as fructose-6-phosphate pathway, was statistically significantly enriched in the case group compared to the control group (Supplementary Table S5B). Additionally, the microbial communities in the case group were significantly enriched in genes encoding for enzymes such as xylulose-5-phosphate/fructose-6-phosphate phosphoketolase (K01621), which is key enzyme in fructose-6-phosphate pathway, as well as pyridoxamine 5’-phosphate oxidase (K00275), and dihydropteroate synthase (K00796) (Supplementary Table S5C).

The relative abundance of bacterial species, functional genes, and pathways in relation to the CD4 + T-cell counts were analyzed by LDA, Wilcoxon rank sum test, as well as Spearman’s correlation analysis. LDA identified five species biomarkers (shown in Supplementary Figure 4 including Lactococcus lactis) for the lower CD4 + T-cell counts (LC) group and one (Alistipes putredinis) for the higher CD4 + T-cell counts (HC) group. Wilcoxon rank sum test was consistent with LDA with two exceptions (Supplementary Table S6A). Consistently, A. putredinis was positively correlated to both nadir CD4 + T-cell count (rho = 0.603, p = 0.0290) and current CD4 + T-cell count (rho = 0.779, p = 0.0017), while Lactococcus lactis was inversely correlated to nadir CD4 + T-cell count (rho = −0.733, p = 0.0044) and current CD4 + T-cell count (rho = −0.717, p = 0.0058) using Spearman’s correlation analysis (Figure 5A). The total ART duration was positively correlated to Ruminococcus lactaris (rho = 0.580, p = 0.0375) while inversely correlated to Streptococcus sanguinis (rho = −0.638, p = 0.0191), Streptococcus gordonii (rho = −0.554, p = 0.0497), Megamonas unclassified (rho = −0.620, p = 0.0237), and Aggregatibacter unclassified (rho = −0.559, p = 0.0470) (Figure 5A).

LDA identified one MetaCyc pathway biomarker, i.e., phosphopantothenate biosynthesis I (PANTO-PWY) among patients in the lower current CD4 + T-cell counts (LC) group compared to those with higher current CD4 + T-cell counts (HC) group (data not shown). Using Wilcoxon rank sum test, 22 additional MetaCyc pathways and one KEGG pathway chlorocyclohexane and chlorobenzene degradation (ko00361) were found to be significantly enriched in the LC group (Supplementary Tables S6B,C). At the gene family level, 57 KOs with characterized function were found to be significantly different between LC and HC groups (Supplementary Table S6D). The Spearman’s correlation analysis showed that three MetaCyc pathways were inversely correlated to both nadir and current CD4 + T-cell counts (Figure 5B). Of these, phosphopantothenate biosynthesis I (PANTO-PWY) showed strongly inverse correlations with both nadir CD4 + T-cell count (rho = −0.784, p = 0.0015) and current CD4 + T-cell count (rho = −0.780, p = 0.0017). Consistently, phosphopantothenate biosynthesis I (PANTO-PWY) and glycerol degradation to butanol (PWY-7003) pathways were significantly enriched in the LC group (p = 0.0012 and p = 0.0082, respectively) by Wilcoxon rank sum test (Supplementary Table S6B). These data suggested that MetaCyc pathways related to phosphopantothenate biosynthesis I and glycerol degradation to butanol were associated with low CD4 + T-cell count. We also observed correlations between three KEGG pathways and CD4 + T-cell counts (Figure 5C). For instance, linoleic acid metabolism pathway (ko00591) was inversely correlated to both nadir CD4 + T-cell count (rho = −0.630, p = 0.0211) and current CD4 + T-cell count (rho = −0.577, p = 0.0391), which was consistent with an earlier study (Vesterbacka et al., 2017).

We identified multiple bacterial virulence factors genes and AMR genes in the full dataset (Supplementary Tables S7, S8). It is noteworthy that more than 200 virulence factor genes were detected exclusively in HIV-1-ART cases, though their relative abundance was low. These genes encode critical bacterial virulence factors such as adhesion, toxin, hemolysin, type III secretion system, etc. (Supplementary Table S7). We found that 173 AMR genes were present in cases only; in addition, cases were significantly enriched in genes associated with tetracycline antibiotic resistance and antibiotic efflux pumps (p < 0.05) (Supplementary Table S8). Rarefaction analysis showed that sequencing depth critically affected the recovery of bacterial virulence factors genes and AMR genes, which was dramatically different from taxonomic profiling. A sequencing depth of 57 million reads per sample was insufficient to capture full spectrums of bacterial virulence factors genes and AMR genes in human fecal samples in this study (Supplementary Figures 5a,b).

We observed correlations between virulence factors genes, AMR genes, and clinical variables (Figures 5D,E). Both current and nadir CD4 + T-cell counts were correlated inversely to virulence genes encoding 6-phosphogluconate dehydrogenase (gnd), E. coli common pilus structural subunit EcpA (ecpA), ferrienterobactin outer membrane transporter (fepA), and oxygen-independent coproporphyrinogen III oxidase (chuW) (Figure 5D). While nadir CD4 + T-cell count was inversely correlated to several important virulence genes such as type VI secretion system ATPase TssH (tssh) (rho = −0.625, p = 0.0224), fimbria adhesin EcpD (ecpD) (rho = −0.596, p = 0.0317), type II secretion system protein L (gspL) (rho = −0.702, p = 0.0074), polysialic acid transport protein KpsM (kpsM) (rho = −0.588, p = 0.0346), and enterobactin synthase component E (entE) (rho = −0.556, p = 0.0483). Additionally, current CD4 + T-cell count was inversely correlated to genes encoding outer membrane usher protein EcpC (ecpC) (rho = −0.652, p = 0.0157), type II secretion system protein C (gspC) (rho = −0.765, p = 0.0023), permease of iron compound ABC transport system (shuU) (rho = −0.592, p = 0.0330), and type 1 fimbriae regulatory protein FimB (fimB) (rho = −0.584, p = 0.0362). Our study also indicated that HIV-1 infection is correlated to changes in AMR genes. Both current and nadir CD4 + T-cell counts were inversely correlated to genes associated with beta-lactam antibiotic resistance (TEM-108 and TEM-219) and peptide antibiotic resistance (icr-Mo and ugd) (Figure 5E). While nadir CD4 + T-cell count was inversely correlated to tetracycline antibiotic resistance gene tetX (rho = −0.694, p = 0.0085), current CD4 + T-cell count was inversely correlated to mdtA (rho = −0.624, p = 0.0226), which encodes for resistance-nodulation-cell division (RND) antibiotic efflux pump, and glycopeptide resistance gene vanWG (rho = −0.584, p = 0.0362). Total ART duration was strongly positively correlated to aminoglycoside antibiotic resistance gene APH(3’)-Ia (rho = 0.725, p = 0.0050).