Dolutegravir restores gut microbiota in late-stage HIV-1 unlike darunavir: an open-label, randomized clinical trial

Domain: chemistry
Published: 2026-03-04
Editorial outcome: minor_revision

Abstract
Late presentation of HIV-1 infection is linked to gut dysbiosis, impaired immune reconstitution, excess inflammation, immune activation, and increased morbidity and mortality. It is unclear if antiretroviral therapy initiation can reverse HIV-associated gut dysbiosis at all, or if specific antiretroviral regimens are more effective in restoring the gut microbiota than others. This has important implications for the long-term health of individuals with HIV. In this multicenter, open-label, randomized clinical trial (NCT02337322), 88 antiretroviral-naïve individuals with advanced HIV-1 infection (median CD4+ T cells of 34 cells/mm3) were randomized (1:1) to initiate lamivudine/abacavir plus either dolutegravir or ritonavir-boosted darunavir, and were followed for 2 years. Both groups had similar HIV-1 suppression rates and recovery of CD4+ T cells. However, treatment with dolutegravir led to increased gut microbial richness and diversity and enrichment of specific microbial taxa and metabolic pathways. These changes were associated with reduced inflammation and lower immune activation, outcomes that did not occur with darunavir/ritonavir. After two years, participants on dolutegravir-based therapy had gut microbiota profiles more closely resembling those of people without HIV, compared to individuals taking darunavir/ritonavir. In summary, dolutegravir-based therapy restores the gut microbiota more effectively than darunavir/ritonavir in patients who present late with HIV.

Introduction
Nearly half of all new diagnoses of human immunodeficiency virus type 1 (HIV-1) infection occur in individuals with advanced disease and severe immunological deterioration, i.e., CD4+ T cell counts <350 cells/mm3 1,2. Initiating antiretroviral therapy as early as possible is crucial for improving prognosis 3,4. However, often this only results in partial immune reconstitution 5,6. Late presenters remain at high risk of opportunistic infections, malignancies, and death 5,7,8, even years after starting antiretroviral therapy 3,4. It is thus essential to implement measures for the early diagnosis and treatment of HIV-1 infections, identify optimal antiretroviral treatment options for late presenters, and uncover modifiable factors that may help mitigate the comorbidity gap associated with late presentation. Persistent gut dysbiosis is commonly found in individuals with low nadir CD4+ T cell counts 9, and may partially account for such a gap.

The gut-associated lymphoid tissue (GALT) plays a critical role in human health, as well as in HIV pathogenesis 10. It houses the majority of CD4+ T cells in the body, serves as the largest site for HIV-1 replication, and is essential for maintaining intestinal barrier integrity 11–13. A diverse and rich gut microbiome is crucial for human health 14–17, as commensal microbes help preserve mucosal integrity 18, modulate immune homeostasis 16, and promote enterocyte health 19,20. Disruption of enterocytes and GALT during acute HIV-1 infection leads to intestinal and systemic inflammation, immune dysregulation, and the translocation of bacterial products into the bloodstream, which exacerbates chronic systemic inflammation 21. Chronic gut inflammation and immune dysregulation result in decreased microbial richness, an increase in generalist microbes better adapted to oxidative stress, and a loss of highly specialized oxygen-sensitive species 9,22,23.

The initiation of antiretroviral treatment generally enhances immune reconstitution. However, its effect on restoring gut dysbiosis remains uncertain. Additionally, it is not clear if certain antiretroviral regimens with superior intestinal penetration, such as HIV integrase strand-transfer inhibitors, are more effective in restoring gut microbiota compared to those with lower intestinal penetration, like protease inhibitors 12,24,25. This has significant implications for the long-term well-being of individuals with HIV, as a rich and diverse gut microbiota is one of the key indicators of better overall health 26. Here, we examined these questions in a 2-year, two-arm randomized clinical trial comparing first-line antiretroviral therapy with dolutegravir versus ritonavir-boosted darunavir/ritonavir in subjects presenting with advanced-stage HIV. The primary and most secondary outcomes of this trial have been previously published 27; here, we report a pre-specified secondary outcome focusing on the longitudinal restoration of gut microbiota composition and function.

Results
Similar immune reconstitution, but greater reductions in inflammation and immune activation with dolutegravir—than with darunavir/ritonavir—including ART

Eighty-eight ART-naïve participants with baseline CD4+ T cell counts below 100 cells/mm3 (median, IQR: 34, 12–67) were randomized 1:1 to receive either darunavir/ritonavir or dolutegravir, each in combination with lamivudine and abacavir (Fig. 1). Both study groups were well balanced regarding demographic factors, body mass index, soluble immune activation markers, cytokines, and microbial gene richness at study entry (Table 1). Fecal specimens were collected at baseline (week 0) and at weeks 24, 48, and 96. Following DNA extraction and whole-genome shotgun metagenomic sequencing, samples failing DNA quantity requirements or predefined sequencing quality control thresholds were excluded. Metagenomic analyses were therefore restricted to participants with either a high-quality baseline sample or at least two high-quality follow-up samples available. All metagenomic analyses were conducted in the modified intention-to-treat population (mITT2) (Figs. 1 and 2A).

Both treatment arms achieved similar immune reconstitution and viral suppression rates over 96 weeks. CD4+ T-cell counts increased from baseline (Δ = +4.66-fold, 95% CI = 3.44–5.88, q value < 1 × 10−15 in darunavir/ritonavir; Δ = 3.33-fold, 95% CI = 2.23–4.43, q value = 2.5 × 10−11 in dolutegravir). Alongside these improvements, there were also significant CD8+ T-cell increases. Markers of immune activation and inflammation, including TNF-α, IL-6, and sCD14, decreased in both groups (all q value < 0.001). Participants on dolutegravir had a greater reduction in sCD14 at week 96 (Δ = −1.92-fold, 95% CI = −2.12 to −1.73) than those on darunavir/ritonavir (Δ = −1.63-fold, 95% CI = −1.83 to −1.43), with a between-group difference of −0.36-fold (95% CI = −0.62 to −0.10; q value = 0.015). CRP showed a similar pattern (dolutegravir Δ = −1.67-fold, 95% CI = −2.46 to −0.88; darunavir/ritonavir Δ = −0.43-fold, 95% CI = −1.26 to 0.40), but the between-group difference was not statistically significant after false discovery rate (FDR) adjustment (q value = 0.144). Body mass index increased slightly and similarly in both groups at week 96 (Δ = +0.16-fold; 95% CI = 0.89 to 2.261, q value < 0.001) (Supplementary Figs. 1 and 2).

Dolutegravir initiation leads to sustained improvements in microbial gene diversity and overall alpha diversity

Microbial Shannon diversity, observed microbial richness, Gini dominance, and gene richness were assessed next (Fig. 2B). Subjects initiating dolutegravir showed an increase in gene richness by week 48 relative to baseline (Δ = 0.27-fold, 95% CI = 0.06–0.49; q value = 0.004), which was sustained through week 96 (Δ = 0.29-fold, 95% CI = 0.07–0.51; q value = 0.004). Shannon diversity showed a smaller increase, which reached significance at week 96 (Δ = 0.13-fold, 95% CI = 0.01–0.24; q value = 0.025). Gini dominance decreased over time in dolutegravir at week 48 (Δ = −0.007-fold, 95% CI = −0.013 to −0.001; q value = 0.005) and week 96 (Δ = −0.009-fold, 95% CI = −0.015 to −0.003; q value < 0.001). In contrast, no temporal changes reached significance in the darunavir/ritonavir arm (all q value > 0.300). Between-group comparisons (dolutegravir vs darunavir/ritonavir) at weeks 48–96 were directionally positive for gene and observed richness but did not meet the study-wide adjusted threshold (gene richness q value ≈ 0.089; observed richness q value ≈ 0.105).

We then quantified within-arm consistency of longitudinal trajectories. For each participant, we estimated a subject-level slope of change from baseline (Δ per visit) and compared the dispersion of slopes across arms. Across alpha-diversity metrics, subjects on dolutegravir showed lower within-arm dispersion than those receiving darunavir/ritonavir (e.g., gene richness SD dolutegravir = 0.006 vs SD darunavir/ritonavir = 0.008; Shannon diversity SD dolutegravir = 0.003 vs SD darunavir/ritonavir = 0.004). Variance-ratio estimates (dolutegravir/darunavir/ritonavir) were consistently <1—0.495 (95% CI = 0.141–3.040) for gene richness, 0.536 (95% CI = 0.220–2.210) for Shannon diversity, 0.905 (95% CI = 0.252–2.970) for observed richness, and 0.719 (95% CI = 0.303–1.990)—indicating more homogeneous changes within the dolutegravir arm. Robust tests of equality of variances (Brown–Forsythe and Fligner–Killeen) were not significant after FDR correction (q value ≥ 0.317), consistent with limited power, but the directionally consistent pattern supports that mean changes are more readily detected under dolutegravir owing to reduced within-arm heterogeneity (Supplementary Fig. S3A).

Higher microbial richness associates with greater CD4+ counts and lower inflammation in dolutegravir recipients

To explore host–microbiome interactions over time, we then computed repeated-measure correlations between alpha diversity metrics and markers of systemic inflammation, immune activation, and bacterial translocation (Fig. 2C, D). Microbial gene richness correlated positively with CD4+ T cell counts (r = 0.44, q value = 0.004) and body mass index (r = 0.35, q value = 0.047), and inversely with plasma CRP levels (r = −0.43, q value = 0.047) in the dolutegravir arm. Observed richness also showed a positive association with body mass index (r = 0.37, q value = 0.023). These associations suggest that higher microbial richness is linked to improved immune reconstitution and lower systemic inflammation in dolutegravir recipients. No significant correlations were detected in the darunavir/ritonavir arm.

Dolutegravir induces progressive gut microbiota homogenization over 96 weeks

Between-sample beta diversity was first explored by non-metric multidimensional scaling (NMDS) of Bray–Curtis distances (Fig. 3A). Samples were clustered by treatment arm, with 95% confidence ellipses shown for each group and timepoint. To quantify the contribution of clinical covariates to the overall community structure, we performed a PERMANOVA analysis (Fig. 3B). Treatment arm influenced community structure (R² = 0.013, p value = 0.001), with a magnitude of effect similar to that of having an AIDS-defining event (R² = 0.013, p value = 0.001) or cytomegalovirus (CMV) seronegativity (R² = 0.006, p value = 0.001). Timepoints also explained a small but significant portion of the variance (R² = 0.014, p value = 0.003), while age, gender, HIV-1 load, and other covariates were not significant contributors (all p value > 0.100). This indicates that, although each of these factors influences the gut microbiota composition, none of them dominates on its own. Assessment of the taxonomic structure of the top 25 most abundant genera further showed that, as previously described 23,28, the microbiome of men who had sex with men (MSM) was enriched in Prevotella relative to Bacteroides, independent of ART regimen and other factors (Supplementary Fig. 3A).

Moreover, samples from dolutegravir-treated subjects progressively migrated toward a common region in the ordination space over the 96 weeks of follow-up—reflected by increasingly tight confidence ellipses—whereas samples from subjects receiving darunavir/ritonavir remained broadly dispersed. To determine whether this visual convergence represented a real reduction in community heterogeneity, we calculated for each sample the Euclidean distance to the centroid of its treatment- and time-specific group in the NMDS space (Fig. 3C). In the dolutegravir arm, centroid distances decreased markedly from baseline at week 24 (Δ = −0.35-fold, 95% CI = −0.59 to −0.11; q value = 0.001), week 48 (Δ = −0.40-fold, 95% CI = −0.62 to −0.17; q value < 0.001) and week 96 (Δ = −0.47-fold, 95% CI = −0.69 to −0.24; q value < 0.001), indicating progressive microbiota homogenization over time. In contrast, no significant temporal changes were detected in the darunavir/ritonavir arm (all q value > 0.800). Between-group comparisons confirmed that dolutegravir recipients exhibited lower centroid distances than darunavir/ritonavir at weeks 48 (Δ = −0.31-fold, 95% CI = −0.56 to −0.06; q value = 0.030) and 96 (Δ = −0.35-fold, 95% CI = −0.61 to −0.10; q value = 0.027), demonstrating a greater, sustained, and statistically robust reduction in microbial community dispersion under dolutegravir relative to darunavir/ritonavir.

To assess whether the longitudinal convergence observed under dolutegravir represented a shift toward a healthier gut microbiota configuration, we integrated our dataset with the metaHIV cohort 9,28, which includes fecal microbiomes from HIV-negative individuals and several clinical phenotypes of people with HIV (Fig. 3D, E). We computed each sample’s mean Bray–Curtis dissimilarity to the HIV-negative reference group. Both treatment arms remained significantly distinct from HIV-negative controls throughout follow-up (all q values < 0.001). However, dolutegravir-treated participants exhibited a consistent reduction in dissimilarity over time (Δ ≈ −0.05-fold between weeks 0–96), resulting in lower distances than those of darunavir/ritonavir at week 96. Mixed-effects modeling confirmed that, while all HIV-positive profiles—including darunavir/ritonavir and metaHIV subgroups such as concordant, discordant, and elite controllers—differed significantly from HIV-negative individuals (Δ = 0.08–0.14; q value < 0.001), dolutegravir samples at week 96 displayed the smallest deviation (Δ = 0.08 ± 0.02; q value = 8.6 × 10−5). These results indicate that dolutegravir-associated microbiota homogenization after 2 years of treatment partially restores the gut microbial community composition towards that of HIV-uninfected individuals, although full normalization is not achieved and might require longer treatment exposure.

To explore the functional implications of dolutegravir-specific homogenization, we assessed repeated-measures correlations between individual centroid distances and systemic immune or inflammatory markers (Fig. 3F, G and Supplementary Fig. 3B, C). In dolutegravir-treated participants, shorter centroid distances were strongly associated with higher CD4+ T cell counts (r = −0.48; q value = 1.8 × 10−5) and BMI (r = −0.53; q value = 2.5 × 10−6), and inversely correlated with markers of immune activation and inflammation including CD4+CD38+HLA-DR+ T cells (r = 0.50; q value = 0.018), TNF-α (r = 0.40; q value = 0.005), sCD14 (r = 0.38; q value = 0.008) and CRP (r = 0.35; q value = 0.023). No significant associations were observed in the darunavir/ritonavir arm (all q value > 0.180). This suggests that the extent of gut microbiota convergence under dolutegravir is associated with enhanced immune reconstitution and lower systemic inflammation, which aligns with a partial restoration of host–microbiome homeostasis.

Dolutegravir-enriched microbial taxa correlate with immunologic and inflammatory profiles

Next, we sought to understand which specific microbial taxa were involved in the observed variations in gut microbiome ecology and their associations with systemic inflammation and immune activation. A differential abundance analysis at multiple taxonomic levels using ANCOMBC2, and implementing a treatment * time_point interaction model that adjusted for baseline covariates and included patient identifier as a random effect to account for repeated measurements. Although no taxa reached study-wide significance after FDR correction under this stricter longitudinal framework, several genera displayed consistent effect sizes and narrow confidence intervals, suggesting biologically meaningful shifts between treatment arms.

Several dolutegravir-enriched species and genus showed 95% CIs excluding zero at specific visits (Fig. 4A and Supplementary Fig. 5A), including Methanobrevibacter smithii (week 48: log2FC = 2.35; 95% CI = 0.55–4.14), Parasutterella excrementihominis (week 48: log2FC = 3.63; 95% CI = 0.94–6.32), Olsenella scatoligenes (week 48: log2FC = 1.68; 95% CI = 0.23–3.13) and Mogibacterium diversum (week 48: log2FC = 2.49; 95% CI = 0.58–4.40) (Fig. 4A). Conversely, several Bacteroides-associated species were consistently depleted, notably Bacteroides thetaiotaomicron (week 24: log2FC = −3.51; 95% CI = −5.67 to −1.36; week 48: log2FC = −1.99; 95% CI = −3.91 to −0.08; week 96: log2FC = −2.15; 95% CI = −4.00 to −0.30) and Parabacteroides distasonis (week 24: −1.65; 95% CI = −3.03 to −0.28; week 48: −1.57; 95% CI = −2.92 to −0.23).

Although multiple-testing correction in this repeated-measures setting is conservative, the direction of effects is consistent across time, pointing to a reconfiguration characterized by the enrichment of fermentative lineages—including lactate- and succinate-producing genera such as Olsenella, Prevotella, and Parasutterella—alongside the methanogenic archaeon Methanobrevibacter, which can facilitate hydrogen disposal and sustain anaerobic fermentation; in parallel, we observe a relative depletion of Bacteroides-dominated profiles and reductions in Lachnoclostridium/Erysipelatoclostridium. This community shift is compatible with increased cross-feeding toward SCFA production and a more even, less dominance-driven structure, which aligns with the longitudinal increases in richness/evenness and the inverse associations with systemic inflammation.

To explore potential host–microbiome interactions, we next assessed the association between longitudinal changes in microbial species and systemic immunologic and inflammatory markers using repeated-measures correlation (Fig. 4B and Supplementary Fig. 5B). Despite modest overall effect sizes, several taxa exhibited consistent and biologically coherent relationships within treatment arms. Within the dolutegravir arm, Methanobrevibacter smithii abundance was positively correlated with BMI (r = 0.32; q value = 0.017), while Parasutterella excrementihominis showed an inverse association with plasma sCD14 levels (r = −0.42; q value = 0.012), consistent with its previously reported anti-inflammatory profile 29–31. Conversely, Lactobacillus sakei (r = −0.31; q value = 0.022) and Clostridium aldenense (r = −0.32; q value = 0.017) both exhibited negative correlations with BMI, and Olsenella scatoligenes correlated positively with BMI (r = 0.30; q value = 0.022) but inversely with CRP (r = −0.39; q value = 0.046), suggesting context-dependent metabolic and inflammatory associations. Alistipes indistinctus negatively correlated with activated CD8+CD38+HLA-DR+ T cells (r = −0.59; q value = 0.014), indicating a potential link between this taxon and reduced T cell activation. In contrast, under darunavir/ritonavir, only Parabacteroides distasonis displayed significant correlations, positively associated with BMI (r = 0.34; q value = 0.039). Together, these findings indicate that the immunometabolic landscape accompanying dolutegravir initiation is intertwined with species-specific microbial trajectories, particularly involving hydrogenotrophic archaea and commensals with fermentative or anti-inflammatory capacities.

Dolutegravir remodels gut microbial co-occurrence networks and increases local motif cohesion by week 96

At week 96, co-occurrence networks revealed clear topological differences by ART regimen (Fig. 5A). In participants treated with dolutegravir, the network divided into 32 connected components, whereas under darunavir/ritonavir treatment, it divided into 55 components. The largest connected component (LCC) included 53% of taxa in the dolutegravir group compared to 15% in the darunavir/ritonavir group, indicating a more consolidated core in the former. Global network metrics further highlight this increased cohesion: dolutegravir networks demonstrated a non-zero clustering coefficient (0.127 compared to 0.000), slightly greater density (0.0191 compared to 0.0076), and marginally greater natural connectivity (0.0151 compared to 0.0145). Eigenvector-based hub analysis at week 96 identified Bacteroides, Gordonibacter, Lachnoclostridium, and Parabacteroides as central nodes in the dolutegravir network, whereas darunavir/ritonavir hubs included Catenibacterium, Eggerthella, Erysipelatoclostridium, and Holdemanella, with fewer cross-phyla links. Together, these measures indicate that long-term dolutegravir treatment is associated with a more integrated and cohesive network structure.

To better understand how local structure in the network is organized, we performed a graphlet correlation analysis, which examines how often small patterns of connections—called graphlets—appear together (Fig. 5B). This approach allowed us to detect two distinct types of structural motifs in the dolutegravir network. One group of motifs was made up of linear, chain-like patterns, while the other formed more circular, loop-like structures. These motifs tended to co-occur consistently in the dolutegravir network, as shown by strong correlations between their frequencies (Spearman |ρ| > 0.8, p < 0.01). In contrast, the darunavir/ritonavir network also exhibited several strong motif–motif correlations, but these were fewer in number, and the overall correlation structure was more dispersed, indicating a less coherent local organization. Quantitative comparison of graphlet correlation distances confirmed that the week 96 dolutegravir network diverged significantly from both baseline and darunavir/ritonavir architectures, supporting a more pronounced structural reorganization under dolutegravir. This demonstrates that 96 weeks of dolutegravir not only increased global connectivity but also promoted the emergence of more stable and internally coherent higher-order interaction motifs. Network and graphlet analysis for weeks 0, 24, and 48 showed intermediate connectivity and motif organization that progressively converged toward the structures observed at week 96 (Supplementary Fig. 5A, B).

Dolutegravir-specific remodeling of metabolic and biosynthetic pathways is linked to host immune and inflammatory markers

Applying ANCOMBC2 to MetaCyc pathways with the same treatment * time_point model, including a random intercept for participant, no pathway met the pre-specified multiple-testing criteria after adding random effects (q value < 0.05; prevalence > 10%). However, several functions showed visit-specific effect sizes whose 95% CIs excluded zero (Fig. 6A). At Week 24, dolutegravir was associated with higher relative abundances of central-carbon and biosynthetic routes, including the incomplete reductive TCA cycle (P42-PWY; log2FC = 0.54; 95% CI = 0.03–1.06), phospholipid biosynthesis I (PHOSLIPSYN-PWY; log2FC = 0.52; 0.05–1.00), tetrapyrrole biosynthesis II (PWY-5189; log2FC = 0.97; 95% CI = 0.36–1.57), TCA cycle II (PWY-5690; log2FC = 0.50; 95% CI = 0.09–0.91), and ubiquinol-8 biosynthesis (UBISYN-PWY; log2FC = 1.36; 95% CI = 0.24–2.48). In the reverse comparison (darunavir/ritonavir > dolutegravir), we noted increased activity in the urea cycle (PWY-4984; log2FC = −0.89; 95% CI = −1.58 to −0.19), chondroitin sulfate degradation I (PWY-6572; log2FC = −1.48; 95% CI = −2.53 to −0.43), acetyl-CoA fermentation to butanoate II (PWY-5676; log2FC = −0.48; 95% CI = −0.94 to −0.01), and 4-hydroxybenzoate biosynthesis V (PWY-6435; log2FC = −1.44; 95% CI = −2.59 to −0.30).

By Week 48, dolutegravir showed higher L-glutamine biosynthesis III (PWY-6549; log2FC = 0.36; 95% CI = 0.05–0.66) and methanogenesis from H2/CO2 (METHANOGENESIS-PWY; log2FC = 1.30; 95% CI = 0.22–2.38), whereas darunavir/ritonavir showed higher polyisoprenoid biosynthesis (POLYISOPRENSYN-PWY; log2FC = −0.62; −1.06 to −0.18), urea cycle (PWY-4984; log2FC = −0.83; −1.44 to −0.22) y fatty-acid α-oxidation III (PWY66-388; log2FC = −1.82; 95% CI = −3.17 to −0.46). At Week 96, additional differences included lower creatinine degradation II (PWY-4722; log2FC = −0.95; 95% CI = −1.87 to −0.03), sulfoglycolysis (PWY-7446; log2FC = −1.47; 95% CI = −2.92 to −0.03) y L-glutamate → propanoate VIII (PWY-5088; log2FC = −0.73; 95% CI = −1.36 to −0.11) under dolutegravir (Fig. 6A).

At Week 96, dolutegravir-associated differences persisted but shifted towards pathways involved in nitrogen, sulfur, and amino acid metabolism. Specifically, dolutegravir showed lower relative abundances of creatinine degradation II (PWY-4722; log2FC = −0.95, 95% CI = −1.87 to −0.03), sulfoglycolysis (PWY-7446; log2FC = −1.47, 95% CI = −2.92 to −0.03) and L-glutamate → propanoate VIII (PWY-5088; log2FC = −0.73, 95% CI = −1.36 to −0.11), suggesting a late-phase relative depletion of nitrogen and sulfur turnover routes under dolutegravir.

Applying repeated-measures correlations, we identified treatment-specific couplings between microbial functions and host status (Fig. 6C). Under dolutegravir, CD4+ T cell counts tracked positively with central-carbon and biosynthetic pathways—TCA cycle II (ρ = 0.49, q value = 2.1 × 10−5), methanogenesis from H2/CO2 (ρ = 0.37, q value = 5.4 × 10−3), and L-glutamine biosynthesis III (ρ = 0.30, q value = 4.8 × 10−2). BMI also correlated with polyamine biosynthesis II (ρ = 0.43, q value = 1.6 × 10−3), L-glutamate → propanoate degradation (ρ = 0.37, q value = 5.9 × 10−3) and methanogenesis (ρ = 0.39, q value = 3.9 × 10−3). In contrast, under darunavir/ritonavir, BMI associated positively with the urea cycle (ρ = 0.38, q value = 2.2 × 10−2) and negatively with ketogenesis (ρ = −0.37, q value = 2.2 × 10−2) and creatinine degradation II (ρ = −0.34, q value = 4.2 × 10−2).

Collectively, these functional data indicate that dolutegravir induces an anabolic remodeling of gut microbial function, enriching central-carbon and biosynthetic pathways that align with higher CD4+ T cell counts, while darunavir/ritonavir shows narrower, largely catabolic shifts. In line with these findings, treatment-specific repeated-measures correlations indicate that dolutegravir-enriched functions are associated with improved immune status. In contrast, under darunavir/ritonavir therapy, body mass index (BMI) is correlated with urea cycle activity and inversely linked to ketogenesis, which may reflect a more catabolic and potentially pro-inflammatory state.

Methods
Study design
This substudy was nested within the ADVANZ-4 trial (EudraCT: 2014-002281-70, ClinicalTrials.gov: NCT02337322), a prospective, open-label, two-arm randomized clinical trial enrolling ART-naïve adults presenting with CD4+ T cell counts below 100 cells/mm3 at HIV diagnosis 27. Participants were randomized 1:1 to receive either dolutegravir (50 mg QD) or ritonavir-boosted darunavir (800/100 mg QD), each combined with a fixed-dose formulation of lamivudine (300 mg QD) and abacavir (150 mg QD). All participants received cotrimoxazole prophylaxis until achieving CD4+ counts >200 cells/mm3 for 3 consecutive months. Follow-up visits were conducted at baseline and weeks 24, 48, and 96. Fecal samples for microbiome analysis were collected at these time points. We employed a modified intention-to-treat 2 (mITT2) approach, restricting the microbiome analysis to participants providing fecal samples at baseline or at least two follow-up time points.

This study adheres to the STORMS reporting guidelines for human microbiome research, and the completed checklist is provided as Supplementary Table 1.

Endpoints and visit schedule
According to the Statistical Analysis Plan (SAP) of the ADVANZ-4 trial, the primary endpoint was the change in absolute CD4+ T-cell count from baseline to week 48. Secondary endpoints included: (i) plasma HIV-1 RNA levels and the proportion of participants achieving viral suppression (<50 copies/mL) at all scheduled visits; (ii) changes in lymphocyte subsets, clinical laboratory parameters, and immune activation, inflammation, apoptosis, and bacterial translocation markers (TNF-α and sCD14), measured at baseline, week 24, and week 48; and (iii) the occurrence of clinical outcomes including new AIDS-defining events, immune reconstitution inflammatory syndrome, adverse events, treatment discontinuations, and mortality.

Per protocol and SAP specifications, clinical assessments, blood sampling, and virological measurements were performed at baseline and weeks 4, 12, 24, 36, and 48, with extended follow-up to week 96 in this substudy. Immunological and inflammation biomarkers were assessed at baseline, week 24, and week 48, consistent with the trial protocol.

Participant exclusions
A total of 45 screened individuals were excluded prior to randomization. Most exclusions were due to medical decision (n = 31), primarily involving severe AIDS-defining conditions that compromised clinical stability (n = 11; including opportunistic infections requiring ICU admission [n = 5, of whom two died during admission], primary CNS lymphoma [n = 3], non-Hodgkin lymphoma [n = 2], and CNS toxoplasmosis [n = 1]). Additional medical reasons included severe potential drug–drug interactions (n = 8; mainly tuberculosis or psychiatric comorbidities), absence of social support (homelessness; n = 4), anticipated lack of adherence (n = 3), liver disease (n = 3), and psychiatric disorders (n = 2). Seven individuals did not meet the inclusion criteria (positive HLA-B*5701 test, n = 4; confirmatory CD4+ T cell count > 100 cells/mm3, n = 3), and none had baseline antiretroviral resistance. Finally, seven individuals declined to participate.

Sample collection and processing
Stool samples were collected either at clinical sites or at home using the OMNIgene. GUT collection kit (DNA Genotek, Ottawa, Canada). The kit stabilizes microbial DNA at ambient temperatures, minimizing degradation during transport. Upon receipt at the central laboratory (IrsiCaixa, Badalona, Spain), samples were immediately frozen at −80°C. DNA extraction was performed using the PowerSoil DNA Isolation Kit (MOBIO Laboratories, Carlsbad, CA, USA) following the Earth Microbiome Project standard protocols 56, ensuring reproducibility and minimizing technical biases. Mechanical lysis was achieved via bead-beating.

Sequencing libraries were prepared using the Nextera XT DNA Library Prep Kit (Illumina Inc., San Diego, CA, USA) targeting an average insert size of 300 bp. Paired-end sequencing (2×150 bp) was performed on the Illumina HiSeq 2500 platform. A sequencing depth target of 20 million reads per sample was established to ensure sufficient coverage for robust metagenomic profiling.

Flow cytometry
Peripheral blood mononuclear cells (PBMCs) were isolated by density-gradient centrifugation using Ficoll-Hypaque (Biomedics-Biomérieux, Madrid, Spain). Immunophenotyping included quantification of CD4+ and CD8+ T-cell subsets and characterization of naïve, effector, and memory populations based on CCR7 and CD45RA expression. Markers of immune activation (CD38, HLA-DR), senescence (CD57, CD28), and apoptosis (annexin V, CD95) were assessed in both CD4+ and CD8+ compartments.

Three-color flow cytometry was performed following previously validated protocols from the parent ADVANZ-4 trial. Data were acquired on a standard flow cytometer and analyzed using FlowJo software. For apoptosis assays, PBMCs were cultured for 18 h at 1×106 cells/mL in 24-well plates in the presence or absence of 1% phytohemagglutinin, stained with annexin V-PE (BD PharMingen, San Jose, CA) in buffer containing 5 mM CaCl2, followed by CD95 detection.

Measurement of inflammatory and immune activation markers by ELISA
Plasma concentrations of inflammatory and immune activation markers were quantified using commercially available enzyme-linked immunosorbent assays (ELISA), performed in accordance with the manufacturers’ protocols. All samples were analyzed under standardized experimental conditions, and measurements were conducted by trained personnel blinded to clinical group allocation.

Baseline plasma samples were available for all participants for all analytes included in the study, and no samples were missing at baseline. Samples were assayed in duplicate, and mean values were used for downstream analyses. Intra- and inter-assay coefficients of variation were within the ranges specified by the manufacturers.

For each analyte, concentration values below the assay-specific lower limit of detection (LLOD) were reported by the assay and retained in the dataset rather than treated as missing values. These values were included in both descriptive and inferential analyses. Apparent gaps in graphical representations therefore reflect measurements below the LLOD and not missing data. As a reference, the LLOD for IL-6 was <2 pg/mL; LLODs for the remaining analytes are provided in the corresponding kit documentation.

All raw concentration values for IL-6 and the remaining inflammatory and immune activation markers measured by ELISA are provided in Supplementary Table 2.

Microbiome data analysis
Raw reads were processed with Trimmomatic (v0.39) 57, removing Illumina adapters, trimming reads with a sliding window of 30 nucleotides, and discarding reads with average Phred scores below 20. Human reads were removed by mapping against the hg19 reference genome using Bowtie2 (v2.4.2) 58 with end-to-end alignment mode and a maximum mismatch penalty of 6. Only unmapped reads were retained.

Taxonomic profiling was performed using MetaPhlAn3 (v3.0) 59 with default parameters, based on clade-specific marker genes from the CHOCOPhlAn database (version 201901). Reads were additionally mapped to the Integrated Gene Catalog (IGC) 60 using Bowtie2 for gene-level quantification. Rarefaction curves were generated by subsampling reads to increasing coverage thresholds, and gene richness was defined at the minimum point where curves plateaued. Samples falling below the 2nd percentile of the maximum coverage distribution were excluded to mitigate depth-related biases.

Microbiome data processing
Post-taxonomic assignments, abundance matrices were imported into R as Phyloseq objects (phyloseq v1.42.0) 61 for preprocessing. Samples were filtered for completeness of metadata and microbiome data. Taxonomic abundance, taxonomy tables, and sample metadata were integrated into a TreeSummarizedExperiment (TreeSE) object (tidyTreeSummarizedExperiment v1.10.0), facilitating consistent downstream analysis.

Alpha diversity analysis
Alpha diversity was evaluated using the mia R package (v1.8.0). Metrics computed included Shannon diversity, Gini dominance, and observed richness (number of detected taxa). To account for uneven sequencing depth, a rarefaction strategy was implemented at the 2nd percentile of library sizes, repeated over 10 random iterations to ensure robustness.

Group and longitudinal differences were analyzed with linear mixed-effects models with the following formula:

diversity metric ~ treatment * timepoint + gender + age + riskgroup + (1|participant_id)

This model included fixed effects for treatment, timepoint, and their interaction, and a random intercept for participant. Multiple testing was controlled using the Benjamini–Hochberg FDR, and model-based effect sizes with 95% confidence intervals are reported. Distributions were visualized with boxplots overlaid with jittered data points and estimated marginal means for interpretability.

To characterize within-arm consistency of longitudinal trajectories, we additionally estimated a subject-level slope of change from baseline (Δ per visit) for each diversity metric and compared the dispersion of slopes between treatment arms. Equality of variances was assessed using Brown–Forsythe and Fligner–Killeen tests, and variance ratios (dolutegravir/darunavir/ritonavir) with bootstrap confidence intervals (5000 resamples) were computed to quantify differences in intra-group variability. As a confirmatory analysis within the mixed-model framework, we fitted location–scale linear mixed models allowing arm-specific residual variances (varIdent structure in the nlme package) and compared these to homoscedastic models using likelihood ratio tests.

Gene richness was separately assessed by counting the number of genes mapped above a dynamic threshold determined per sample, ensuring comparability across the cohort. This dynamic threshold was derived based on the rarefaction analysis against the IGC, selecting for each sample the coverage depth at which the rarefaction curve plateaued, thus avoiding biases due to differences in sequencing depth.

A detailed overview of all statistical analyses, including the family of tests and the number of comparisons used for FDR control in each analysis, is provided in Supplementary Table 3.

Correlation analyses
Repeated-measures correlation effects (rmcorr), which account for within-participant dependency across longitudinal samples, were calculated to explore associations between microbiome alpha diversity indices (Shannon diversity, Gini dominance, observed richness, and gene richness) and a comprehensive panel of clinical parameters. The clinical variables included circulating immune cell counts (CD4+, CD8+, CD4 nadir), immune activation markers (expression of CD38 and HLA-DR on CD4+ and CD8+ T cells), systemic inflammatory markers (C-reactive protein [CRP], interleukin-6 [IL-6], tumor necrosis factor-alpha [TNF-α], and soluble CD14 [sCD14]), and body mass index (BMI).

To normalize the distribution of cytokine levels and mitigate the influence of extreme values, cytokine data were log2(x+1) transformed prior to analysis. Repeated-measures correlations were computed both across the full cohort and within each treatment arm to capture treatment-specific host–microbiome associations.

Two-sided p values from correlation tests were adjusted using the Benjamini–Hochberg method to control for the FDR. Only associations with q < 0.05 were considered significant. Results are presented as correlation coefficients (r) and visualized in heatmaps, with color encoding the strength and direction of associations. Selected relationships were further illustrated with scatterplots showing participant-level trends and fitted regression lines.

Longitudinal evolution of clinical marker levels
Longitudinal changes in clinical markers were analyzed using linear mixed-effects models with the specification:

marker ~ treatment * timepoint + gender + age + riskgroup + (1|participant_id)

Here, time_point was modeled as a categorical factor to allow non-linear trajectories, and a random intercept for participant_id accounted for within-individual correlation. The treatment-by-time interaction enabled inference on both between-group differences at each timepoint and within-group longitudinal changes. Model coefficients (β) and 95% confidence intervals were estimated, and planned contrasts (pairwise time comparisons within treatment and between-treatment comparisons at each time point) were derived from the fitted models. Multiple testing was controlled using the Benjamini–Hochberg FDR. Visualization included boxplots of raw marker distributions by time and treatment with jittered points, forest plots summarizing model-based contrasts with FDR-adjusted p values, and longitudinal effect plots showing time on the x-axis and estimated effects (with 95% CIs) on the y-axis to highlight temporal trends.

Beta diversity analysis
Beta diversity was computed with Bray–Curtis dissimilarities from relative-abundance profiles. Community structure was visualized by NMDS (monoMDS in vegan). To test associations while accounting for repeated measurements, we fit PERMANOVA models with adonis2 using:

Distance ~ treatment * timepoint + gender + age + aidsevent + cmvnegative + riskgroup + HIVVL

Constraining permutations within participants (strata = record_id). Only complete cases for model covariates were included.

To summarize within-group dispersion, we calculated each sample’s Euclidean distance in NMDS space to its leave-one-out centroid for the corresponding treatment × time_point group (“centroid distance”). Centroid distances were analyzed with linear mixed-effects models (lmerTest) to estimate (i) longitudinal changes within treatment arms versus baseline (week 0) and (ii) between-treatment differences at each timepoint, using:

centroiddist ~ treatment * timepoint + gender + age + riskgroup + (1|participant_id)

Planned contrasts were obtained with emmeans; we report estimates, 95% CIs, and Benjamini–Hochberg FDR-adjusted p values. Results are shown as forest plots (between-treatment contrasts) and longitudinal trajectory plots (within-treatment changes vs baseline).

Associations between community variation and clinical/immunologic markers (CD4, CD8, activated T cell subsets, CRP, IL-6, TNF-α, sCD14, BMI) were assessed using repeated-measures correlation (rmcorr), reported as correlation coefficients (r) with FDR-adjusted p values, both globally and stratified by treatment. For compositional summaries, genus-level stacked barplots (top 25 genera; remainder “Other”) were ordered by NMDS1 and annotated with treatment, risk group, timepoint, and centroid distance. Figures were generated with ggplot2 and arranged with patchwork; PERMANOVA tables were rendered with gt. All analyses were performed in R.

Differential abundance analysis
We first benchmarked microbiome differential-abundance software using the dar R package to guide method selection 62. Differentially abundant taxa were then identified with ANCOMBC2 (v2.1.2) (REF) on TreeSummarizedExperiment objects at Genus and Species resolution. The model specification included fixed effects for treatment, time, and covariates, and a subject-level random intercept to account for repeated measures:

Abundance ~ treatment * timepoint + gender + age + riskgroup + (1|record_id)

Pairwise treatment contrasts at each time point were obtained by setting group = “treatment” and pairwise = TRUE. Analyses applied a prevalence filter of ≥10% across samples and used Benjamini–Hochberg FDR control. We considered taxa significant if the FDR-adjusted q value < 0.05 and the 95% confidence interval for the log fold-change excluded 0. Results are reported as effect sizes (logFC) with 95% confidence intervals and visualized as time-stratified forest plots and heatmaps of logFC (with a side annotation for prevalence).

To relate differential taxa to host immunologic and clinical markers (CD4+ and CD8+ T cell counts, activated T cell subsets, CRP, IL-6, TNF-α, sCD14, BMI), we computed repeated-measures correlations (rmcorr) between CLR-transformed taxon abundances and each marker. Cytokines were log2-transformed with a +1 offset prior to analysis. Repeated-measures correlation coefficients (r) and FDR-adjusted p values were summarized overall and by treatment arm and displayed as dotmaps and per-marker rmcorr plots.

Significance criteria were: |log2 fold-change| ≥ 0.5, FDR-adjusted p value < 0.05, and prevalence > 20%. Differential features were visualized using volcano plots and annotated barplots. Significant taxa were additionally correlated with clinical parameters using Spearman’s method.

Functional analysis
Functional profiling focused on MetaCyc pathways derived from HUMAnN3 outputs 59. We retained non–taxa-stratified pathways and excluded “UNMAPPED/UNINTEGRATED” entries. Pathway abundance tables were aligned to the study metadata and assembled into TreeSummarizedExperiment objects. Differential pathway abundance was tested with ANCOMBC2 using the model:

Abundance ~ treatment * timepoint + gender + age + riskgroup + (1|record_id)

With a subject-level random intercept (record_id) to account for repeated measures, group = “treatment”, Benjamini–Hochberg FDR control, and a prevalence filter of ≥10%. We considered effects significant when the FDR-adjusted q < 0.05 and the 95% CI for the log fold-change excluded 0. Results are reported as logFC with 95% CIs and visualized as time-stratified forest plots and heatmaps of logFC (including a prevalence side annotation).

To relate functional variation to host markers (CD4+ and CD8+ T cell counts, activated T cell subsets, CRP, IL-6, TNF-α, sCD14, BMI), we computed repeated-measures correlations (rmcorr) between CLR-transformed pathway abundances and each marker; cytokines were log2-transformed with a +1 offset. Correlations were summarized globally and by treatment arm, and multiple testing was controlled with FDR.

External cohort integration and comparative dissimilarity analysis (MetaHIV)
To contextualize ADVANZ-4 microbiome profiles against external HIV phenotypes, we integrated samples from the MetaHIV study and reprocessed both datasets with a harmonized pipeline. Short-read profiles were generated with nf-core/taxprofiler using MetaPhlAn4 to obtain species-level relative abundances. Sample metadata were standardized (common IDs, time points, and HIV phenotype labels) and merged. A TreeSummarizedExperiment object was constructed at the species rank, keeping non-missing samples with complete metadata. Beta diversity (Bray–Curtis) was computed from species-level relative abundances. For each sample i, we quantified its dissimilarity to HIV-negative microbiomes as the mean Bray–Curtis distance to all HIV-negative samples:

di = 1/|N| * Σ j∈N BC(i,j)

where N indexes HIV-negative samples. Group differences in diversity were tested with a linear mixed-effects model including a subject-level random intercept:

di ~ Profile + (1|patient_id)

followed by estimated marginal means and treatment-vs-control contrasts with HIV-negative as reference (FDR control via Benjamini–Hochberg). Results are presented as effect estimates with 95% CIs in forest plots.

For ordination, we performed NMDS (Bray–Curtis; monoMDS in vegan) on the integrated species matrix to visualize separation across HIV phenotypes and ADVANZ-4 time points. As a complementary summary, we computed each sample’s Euclidean distance in NMDS space to the HIV-negative centroid:

dist(neg) = sqrt((xi - xneg)^2 + (yi - yneg)^2)

and assessed differences versus HIV-negative with the same mixed-model framework and FDR-adjusted inference.

Network analysis
To explore genus-level co-occurrence patterns under dolutegravir (dolutegravir) versus darunavir/ritonavir (darunavir/ritonavir) therapy, we began with the pre-processed SingleCellExperiment object containing CLR- and log-transformed assays. Samples were subset by timepoint (0, 24, 48, and 96 weeks) and, within each subset, counts were agglomerated at the genus rank. Relative-abundance profiles were generated, and taxa present in fewer than 15% of samples were removed. The remaining genera were then transformed by centered log-ratio (pseudocount = 1) to stabilize variance. For each timepoint, the filtered genus-level data were split by treatment group, yielding separate dolutegravir and darunavir/ritonavir matrices for network inference.

Networks were constructed with the SPRING method in NetCoMi’s netConstruct function (NetCoMi v1.1.0) 63,64, restricting to the 200 most prevalent genera and exploring 20 values of the regularization parameter (nlambda = 20, rep.num = 10, thresh = 0.05, Rmethod = “approx”) on signed dissimilarities. The resulting microNet objects were characterized using netAnalyze, applying Louvain community detection, identifying hub taxa via eigenvector centrality, reporting unnormalized degree centrality on the LCC, and computing graphlet orbit counts together with the Graphlet Correlation Matrix (GCM) for that component. To highlight treatment-specific shifts in network micro-architecture, we calculated the difference between the dolutegravir and darunavir/ritonavir GCMs (ΔGCM = GCM_dolutegravir – GCM_darunavir/ritonavir) and visualized this differential matrix as a single heatmap, with statistical significance of Δρ values assessed by FDR-adjusted Spearman tests.

Software and reproducibility
All analyses were conducted in R (v4.5.0) using packages from the Bioconductor (v3.21) and CRAN repositories. Scripted workflows were version-controlled with Git, and data processing pipelines were documented using Quarto notebooks, ensuring transparency. Parallel computing was leveraged via the BiocParallel or future R packages. R SessionInfo and packages are provided in each Quarto notebook to facilitate reproducibility.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript presents a subanalysis of the ADVANZ-4 trial evaluating microbiome and inflammatory changes with dolutegravir versus darunavir/ritonavir in advanced HIV infection. Immunovirological endpoints and most secondary outcomes did not differ between arms in the parent ADVANZ-4 trial, which weakens the assertion here that between-group microbiome differences are clinically meaningful. In that larger cohort, among 28 T-cell and soluble biomarkers, only TNF-α and sCD14 showed greater declines with DTG at 48 weeks—effects not FDR-corrected and modest in magnitude; furthermore, the present subanalysis reports higher baseline TNF-α in the DTG arm, complicating interpretation of TNF-α changes.\nThe study addresses an important and timely question on how ART choice may influence the gut microbiome in people with advanced HIV, and the authors are to be commended for conducting a long-term, multicenter randomized trial in a highly challenging patient population. However, several methodological and interpretative issues substantially limit the strength of the conclusions: (i) the parent trial showed largely null immunovirological and inflammatory differences between arms, making the clinical relevance of observed microbiome changes uncertain; (ii) statistical approaches frequently pool repeated measures without accounting for within-subject correlation, apply inconsistent adjustments across analyses, and rely on multiple within-group significance tests that do not address between-group differences; (iii) interpretations sometimes overstate causality or benefit despite modest and exploratory findings, absence of a healthy reference group, and potential confounding from concomitant treatments; and (iv) the abstract and discussion extrapolate beyond the study population and use imprecise or overly strong language.\nTaken together, the parent trial’s largely null biomarker findings argue that the clinical relevance claimed for the microbiome shifts may be overstated.\nMETHODS AND RESULTS\n- The correlation analyses pool repeated measures within the DTG arm, violating independence assumptions, apply inconsistent restrictions by treatment group across microbiome data layers (α-/β-diversity vs. taxa/functions), and should be reanalyzed using methods that account for within-subject correlation, applied consistently across all layers. Lines 120–128, Figs. 2C–D, Fig. S2: Pooling all observations for the DTG arm violates the independence assumption of Spearman’s rank tests. This can inflate the effective sample size and underestimate standard errors, even with FDR adjustment. The rationale for restricting to DTG is unclear; absence of associations in the DRV/c arm would be equally informative. Analyses are inconsistently applied across data layers—α-diversity (Fig. 2D) and β-diversity (Fig. 3) are separated, but taxa (Fig. 4D) and functions (Fig. 6C) are not. More robust approaches include linear mixed-effects models, repeated-measures correlation, or timepoint-specific correlations.\n- Supplementary Figure S1 is difficult to interpret given the numerous within-group and between-timepoint pairwise comparisons. In an RCT, the key question is whether changes from baseline differ between arms, not whether each arm changes from its own baseline. Multiple within-group tests inflate type I error and do not address the primary comparative question. A more appropriate analysis would model change scores (or log fold-changes) with treatment arm as the main predictor, adjusting for baseline values. Presenting between-group mean (or fold) differences with 95% CIs—via line plots—would be clearer and more informative.\n- For consistency, if β-diversity and differential abundance analyses were adjusted for confounders, α-diversity comparisons should be as well.\n- Participants initiating DTG showed greater decreases in sCD14 and CRP, but given the extensive multiple testing, subsampling from the parent cohort, and exploratory nature of these endpoints, these findings should be interpreted as hypothesis-generating. Even in an RCT, small sample sizes warrant adjustment for baseline values and relevant confounders.\n- The PERMANOVA model used does not test whether β-diversity trajectories differ between groups; an interaction term (treatment × time_point) would be required.\n- ANCOM-BC models specifying treatment × time_point including main effects and interaction should account for baseline differences.\n- Concomitant antibiotics, antivirals, and antifungals are potential confounders and should be reported and, if possible, adjusted for.\n- Greater convergence of β-diversity in the DTG arm is not inherently beneficial—similar patterns could occur after broad-spectrum antibiotic exposure. This is further complicated by pooled repeated-measures correlations and negative associations with CD8 counts (a predictor of non-AIDS events), which are difficult to interpret as beneficial.\n- FDR correction: Please clarify, for each analysis, the number of comparisons included in the adjustment. This information should be reported in the Supplementary Methods or a dedicated table.\n- In line with current recommendation (e.g., ASA statements, Nature editorial policies, and the 2021 JAMA statistical guidelines), avoid overuse of “statistically significant” and focus on effect sizes with CIs and clinical relevance.\n- In Supplementary Figure S1, outliers obscure trends; consider splitting the Y-axis to improve visibility.\n- Fig. 3A is difficult to interpret due to overlapping timepoints; separating by timepoint or using trajectory plots could prevent overinterpretation of autocorrelated data.\n- Line 175: Clarify whether the Lachnospira reversal was driven by an increase in DRV/r or decrease in DTG. Correlation analyses in Figs. 4 and S4 are also affected by pooling of repeated measures.\n- Fig. 5A: The claim that DRV/r networks show weaker relationships may be overstated; visually, the differences are not marked. Visually, there appear to be three strong motif–motif correlations in the DRV/r arm and two in the DTG arm, which does not clearly support a marked contrast. Pooling repeated measures without accounting for dependence may bias correlations.\n- Lines 230–245: Consider adding a figure showing distribution of functional categories over time by group. Pooling repeated measures here has the same statistical limitations noted above.\n- Include the STORMS checklist as supplemental material.\n- In Table 1, replace variable labels with proper T-cell subset nomenclature.\n- Report frequency of CMV-negative participants by group.\n- Summarize T-cell activation and ELISA methods in the main text.\nDISCUSSION\n- Lines 281–283: The claim that DTG reduced CRP, TNF-α, and sCD14 is not supported by presented data; mean between-group changes are not shown and differences in the parent study were small. Cite original data rather than outdated reviews, adjust analyses for confounders, and emphasize effect sizes and CIs.\n- Lines 324–325: Moderate the statement that microbial metabolic outputs coincide with reduced inflammation; the ADVANZ data do not convincingly support this.\n- Lines 338–340: Avoid quantitative comparisons with other studies where populations differ substantially (ADVANZ: deeply immunosuppressed; REFRESH: stable PWH on ART with >350 CD4 cells/µL).\n- Discuss darunavir’s microbiota effects in context of recent data (PMID: 38259963), noting limitations of the cited observational work by Villanueava-Millán (unadjusted by determinants or ART choice, hence inadequate to isolate drug-specific effects).\nABSTRACT\n- Provide effect sizes and specific taxa; avoid vague terms like “beneficial taxa.”\n- Replace “restoration of functional pathways” with more neutral language; no reference group represents a “healthy” state.\n- Remove prioritization statement; DTG is already preferred in guidelines for efficacy/safety.\n- Avoid extrapolating results from advanced HIV (<100 CD4 cells/µL) to less advanced infection; note that <350 CD4 cells/µL is the typical definition.\n- Replace “architecture” with “composition.”\n- Suggested closing:\n- “In conclusion, in this cohort of individuals with advanced HIV infection, dolutegravir-based regimens were associated over two years with greater changes in gut microbiome composition and predicted metabolic functions than darunavir/ritonavir-based regimens. The clinical significance of these findings remains to be established.”", "reviewer_2": "This manuscript presents results from a multicenter, open-label, randomized clinical trial comparing the effects of dolutegravir (DTG) versus ritonavir-boosted darunavir (DRV/r), both in combination with abacavir/lamivudine, on the gut microbiome and systemic inflammation in antiretroviral therapy (ART)-naïve individuals with advanced HIV-1 infection. Over 96 weeks, the authors demonstrate that DTG significantly improves gut microbial richness, promotes the expansion of beneficial taxa, and reduces inflammatory markers more effectively than DRV/r.The study is highly relevant to both HIV treatment and microbiome research, employing advanced metagenomic and network analysis techniques. The authors make a strong case for considering microbiome restoration and systemic immune modulation when selecting ART regimens, especially in late presenters.\nStrengths\nClinical and Biological Significance\n• Focus on late-stage HIV fills an important gap in current literature.\n• Demonstrates potential anti-inflammatory advantages of DTG beyond viral suppression.\nComprehensive Longitudinal Analysis\n• 96-week follow-up with multiple timepoints and a modified ITT analysis.\n• Detailed correlation of microbiota metrics with immune activation and inflammation biomarkers (e.g., sCD14, TNF-α, CRP).\nAdvanced Microbiome Methodology\n• Whole-metagenome shotgun sequencing.\n• Functional pathway profiling (MetaCyc) and microbial network analysis (SPRING, graphlet motif evaluation) provide deep insight into microbial community architecture.\nStatistical Rigor and Transparency\n• Multiple testing corrections, multivariate models, and robust correlation analysis.\n• Open access to raw data and source code enhances reproducibility.\nClear Visual Presentation\n• Figures (esp. Figs. 2–6) are visually informative and support all major conclusions.\n• Supplementary data is thorough and well organized.\nWeaknesses\nOpen-label Design\n• While the randomized design is appropriate, the lack of blinding may introduce bias, especially for variables influenced by patient behavior (e.g., diet, adherence).\nLimited Generalizability\n• Participants are predominantly Caucasian males from Spain. Broader inclusion of global populations (e.g., African or Southeast Asian cohorts) would increase external validity.\nDietary and Environmental Confounders\n• The study does not account for dietary habits, which significantly influence gut microbiota. Even basic dietary recall data would have strengthened the analysis.\nMucosal vs. Luminal Microbiome\n• Only fecal samples were analyzed. Mucosal-associated microbiota, which may play a more direct role in gut barrier function, were not assessed.\nClinical Outcomes Not Assessed\n• Although inflammatory biomarkers were significantly reduced in the DTG arm, no data on clinical endpoints (e.g., incidence of non-AIDS comorbidities) were collected. Hence, conclusions about long-term health benefits remain speculative.\nOral Pathobionts Among DTG-Enriched Taxa\n• Some DTG-enriched taxa (e.g., Streptococcus sanguinis, Rothia mucilaginosa) are primarily oral commensals/pathobionts. Their presence in stool samples could reflect translocation or contamination—this deserves clarification.", "reviewer_3": "Thank you for the opportunity to review this manuscript. The authors present a well-designed study examining the impact of switching to dolutegravir (DTG) or boosted darunavir (DRV/r) on the gut microbiota in people with HIV. The randomized design, longitudinal sampling, and integration of microbial, inflammatory, and clinical measures are strengths. The authors show that DTG is associated with consistent increases in several alpha diversity metrics and a shift toward a more homogeneous gut microbial profile, while changes in the DRV/r group are more variable. Given the global shift to DTG-based ART regimens and the metabolic changes often observed (weight gain, blood pressure changes), these findings will be of interest to HIV and microbiome researchers. Some comments to consider are listed below.\nMethods\n- Did any participants switch to DRV/cobicistat throughout the trial?\n- In Figure 1, 31 participants are noted as being excluded for a medical reason. Specifying the nature of this exclusion would provide additional information on the generalizability of the results.\n- For IL-6, baseline measurements are missing for a substantial number of participants. Could the authors clarify whether these were due to values falling below the assay’s limit of lower detection (LLOD) or due to missing samples? If the former, providing the LLOD would assist in interpretation.\n- The manuscript presents cytokine data primarily in the context of correlations with bacterial taxa. It would be helpful to also show longitudinal cytokine trajectories over time for each treatment arm, independent of microbiota associations.\n- In Figure 2, the use of Wilcoxon rank-sum tests for repeated measures may not fully account for within-participant correlation. A paired Wilcoxon or mixed-effects model could be considered.\n- For beta-diversity analyses (Figure 3), PERMANOVA could be stratified by participant ID to avoid inflated Type I error. Clarifying this would strengthen the description of the analysis.\n- For data presented in Figure 5 (and S7), it would be helpful to state whether multiple-testing correction for cytokine–taxa correlations was applied per cytokine or across all tests. Indicating whether correlations were adjusted for repeated measures would also aid in the interpretation of this data.\nResults\n- For Figure 2, consider stating the percentage of total variance explained by the first two NMDS axes in the main text.\n- The minimum prevalence threshold of 0.2 for determining taxa significantly altered from baseline will preferentially retain common taxa, which may introduce bias, particularly where DTG appears to drive microbiome convergence. It would be helpful to include a secondary analysis reducing this filter to 0.1. Several taxa in Figure 4 appear to meet the fold-change and p-value thresholds, but are not highlighted, likely due to the prevalence filter.\n- The use of light grey to represent both baseline (week 0) and non-significant taxa in Figure 4 makes interpretation challenging. Distinguishing these categories visually could improve clarity.\n- Throughout, adding additional information to figure legends would be helpful. For example, specifying the distance metrics used (as well as in the methods section) and indicating what point size represents.\n- Please describe how missing timepoints were handled.\n- Adding the SAP-specified endpoints and visit timings to the methods section would be helpful. Although they were clearly in the protocol, adding to the main body of the paper would be helpful.\nDiscussion\n- The low proportion of female participants may limit generalizability as females have substantial differences in metabolism, gut microbiota, and immune characteristics compared to males. This should be acknowledged as a limitation (and including a statement on this would align with SAGER guidelines).\n- From Figures 2, it appears that participants in the DTG arm experience more consistent increases in alpha diversity, whereas participants who received DRV/r arm experience more heterogeneous changes. It would be helpful to highlight this difference in intra-group variability, as it provides context for why mean changes are larger and more significant in the DTG group.\n- The enrichment of oral taxa in the DTG arm is intriguing. While this could reflect re-establishment of commensals with immune recovery (as the authors suggest), in inflammatory bowel disease increased oral commensals have been linked to decreased barrier function and bacterial translocation. A brief acknowledgment of both possibilities would provide a more balanced interpretation of the findings.\n- The discussion assumes that additional Lactobacillus in the gut is beneficial; however, there are some settings (including immunocompromised individuals, individuals who have recently had antibiotic exposure, and individuals with chronic gut inflammation). The discussion would be strengthened by acknowledging that species- and strain-level information is necessary to conclude if the genus-level increase in Lactobacillus observed in DTG-treated participants is beneficial.\n- Given previous studies that have observed metabolic changes in individuals switching to a DTG-based regimen (weight gain, increased blood pressure – in some populations), it may be of interest to discuss whether the microbiota changes observed here could have metabolic implications."}, "author_response": "Answer 1: We thank the reviewer for this important point. We have reanalyzed all associations to account for within-subject dependence and to ensure methodological consistency across data layers. First, we now use repeated-measures correlation (rmcorr) for all bivariate analyses, which models a common within-participant association across visits and avoids inflating the effective sample size. Second, we added linear mixed-effects models (LMMs) as a regression-based confirmation that accommodates covariates and interaction terms; for each α-diversity metric (and, analogously, CLR-transformed taxa/functions), we fit models of the form outcome ~ treatment * time_point + covariates + (1|participant), with FDR control. Third, results are now presented for both treatment arms (DTG and DRV/r) throughout. For β-diversity, we applied PERMANOVA with participant-constrained permutations. We also include timepoint-specific (per-visit) correlations as an assumption-light sensitivity analysis. Finally, all correlation results are now presented stratified by treatment (DTG and DRV/r) in the main figures, and we also provide unstratified analyses in the Supplementary material for completeness\nAnswer 2: We thank the reviewer for this insightful comment. Following this recommendation, we have restructured the analysis and visualization of longitudinal changes in host immunologic and inflammatory markers. Instead of multiple within-arm pairwise comparisons, we now model change from baseline using linear mixed-effects models (LMMs), specifying outcome ~ treatment × time_point + covariates + (1|participant). This approach directly tests whether changes from baseline differ between treatment arms, while accounting for within-subject correlation and baseline variability.\nThe updated results are presented as estimated mean log₂ fold-changes (Δ from baseline) with 95% confidence intervals derived from the LMMs. Between-group contrasts (DTG – DRV/r) are shown alongside the trajectories to facilitate interpretation of comparative effects. These revisions are now displayed in Supplementary Figure 2, which replaces the previous Supplementary Figure 1 as the main visualization of longitudinal trends.\nSupplementary Figure 1 has been simplified to show the underlying data distributions, while statistical inference is now based solely on the LMM framework.\nAnswer 3: We agree with the reviewer and have updated the α-diversity analyses to ensure consistency with the β-diversity and differential abundance frameworks. Specifically, all α-diversity comparisons are now based on linear mixed-effects models that include the same set of covariates used elsewhere in the study (gender, age, and risk group), in addition to treatment, time point, and their interaction. This adjustment accounts for potential baseline differences and increases precision while preserving the randomized design. The updated model structure is:\n∼ treatment ∗ time_point + gender + age + risk_group + (1 ∣ participant)\nResults and figures have been revised accordingly to reflect these covariate-adjusted estimates.\nAnswer 4: We agree that, although participants initiating DTG showed larger reductions in sCD14 and CRP, the number of comparisons, subsampling, and the exploratory status of these endpoints recommend cautious interpretation. In the revised manuscript we report covariate-adjusted LMM estimates (adjusting for baseline values and relevant confounders) with FDR correction and explicitly frame these results as hypothesis-generating, noting the need for validation in larger cohorts.\nAnswer 5: We thank the reviewer for this helpful observation. We have now updated the β-diversity analysis to formally test whether temporal trajectories differed between treatment arms. Specifically, the PERMANOVA model was expanded to include an interaction term between treatment and time point (treatment × time_point) using the adonis2 function in the vegan package with participant included as a blocking stratum to account for repeated measures.\nAnswer 6: We re-fit the ANCOM-BC models including treatment, time_point, and treatment × time_point, and added patient_id as a random intercept to account for repeated measures. In this setup, the main effect captures baseline differences, and the interaction tests between-arm differences in change from baseline. Methods and Results have been updated accordingly.\nAnswer 7: We agree with the reviewer that concomitant medications can influence microbiome composition. Unfortunately, information on antibiotic or antifungal use was not available for this cohort. However, as all participants were antiretroviral-naïve at enrollment, exposure to other antivirals was uniform across treatment arms, minimizing potential confounding from antiviral therapy.\nAnswer 8: We appreciate the reviewer’s insightful comment and agree that reduced β-diversity dispersion does not necessarily imply a beneficial or “healthier” state. To address this, we have expanded our analysis by integrating data from an independent HIV cohort including uninfected controls. This allowed us to contextualize the observed convergence in DTG recipients. Notably, the DTG arm showed a progressive reduction in dissimilarity to HIV-negative individuals, whereas no such trend was observed for DRV/r. These findings suggest that the microbial community trajectories under DTG move toward a composition more similar to that of uninfected controls, supporting the interpretation of convergence as a shift toward a healthier microbiome configuration rather than a nonspecific loss of diversity.\nAnswer 9: We thank the reviewer for this valuable suggestion. We have now clarified the scope of the FDR correction for each analysis. The number of comparisons included in each adjustment is explicitly reported in the Supplementary Table 3. This table details every analysis type (α-diversity, β-diversity, differential abundance, and correlation analyses), the number of tests, adjustment level, and FDR procedure applied.\nAnswer 10: Agreed. We revised the manuscript to de-emphasize “statistically significant” terminology and instead report effect sizes with 95% CIs and adjusted p-values (q-values). We focus on magnitude, direction, and clinical relevance, including between-arm contrasts of change, and updated figures/legends to show estimates and CIs (no asterisks).\nAnswer 11: Agreed. We have split the Y-axis in the affected panels of Supplementary Figure S1 to improve visibility and clarify temporal trends without distorting the data.\nAnswer 12: Agreed. We have split the plots in Fig. 3A by time point, which improves clarity and prevents overinterpretation of autocorrelated longitudinal data.\nAnswer 13: We thank the reviewer for this observation. Because shotgun metagenomics yields relative (compositional) abundances, it is not possible to determine whether the observed pattern in Lachnospira reflects an absolute increase in the DRV/r arm, an absolute decrease under DTG, or a combination of both. What our reanalysis does show is that the relative abundance of Lachnospira rises over time in participants receiving DRV/r, while remaining comparatively stable in the DTG arm. The revised visualization makes this relative pattern more clearly appreciable, but cannot disentangle absolute abundance changes.\nRegarding the reviewer’s concern about the treatment of repeated measures, all correlation and differential-abundance analyses have now been re-run using mixed-effects models including participant ID as a random effect, which appropriately accounts for within-individual dependence.\nAnswer 14: We thank the reviewer for this comment. We have revised the description to moderate the interpretation of differences between treatment arms, acknowledging that the overall network contrast is subtle. Regarding the concern about pooling repeated measures, we clarify that all network analyses were performed independently at each time point, ensuring that repeated observations from the same participant were not pooled within the same correlation model.\nAnswer 15: We thank the reviewer for this helpful suggestion. We have now included a forest plot summarizing the distribution of functional categories over time, which provides a clearer visualization of the main trends. As noted above, all models have been re-analyzed including patient ID as a random effect, properly accounting for within-subject correlation.\nAnswer 16: The STORMS checklist has now been included as Supplementary Table 1. We have also added a statement in the Study design subsection of the Methods indicating adherence to the STORMS reporting guidelines and referring readers to the supplementary checklist.\nAnswer 17: Table 1 has been updated to use standard T-cell subset nomenclature.\nAnswer 18: Table 1 has been updated to include the frequency of CMV-negative participants by treatment group.\nAnswer 19: We thank the reviewer for the comment. We clarify that T-cell activation was not experimentally induced; rather, activated T-cell subsets were quantified by flow cytometry based on surface marker expression. Details of the ELISA assays used to measure plasma biomarkers have been added to the Methods section.\nAnswer 20: We thank the reviewer for this important comment. We agree that the original wording overstated the evidence regarding the effect of dolutegravir on systemic inflammatory markers. We have revised the text to avoid implying a clear reduction in CRP, TNF-α, and sCD14 attributable to DTG.\nAnswer 21: We agree with the reviewer and have moderated this statement accordingly.\nAnswer 22: We fully agree. Quantitative comparisons with cohorts that differ markedly in immunological status may be misleading. We have therefore removed direct quantitative comparisons between ADVANZ and REFRESH and now explicitly highlight the substantial differences in baseline CD4⁺ T cell counts, treatment history, and clinical context.\nAnswer 23: We thank the reviewer for pointing this out. We have expanded the discussion of darunavir-associated microbiota changes to incorporate recent data (PMID: 38259963).\nAnswer 24: We have replaced vague terms like “beneficial taxa” with “specific microbial taxa and functional pathways”. Due to abstract word limits, listing specific taxa with effect sizes and 95% CIs is not feasible; however, we have implemented this in the Results section, where we report named taxa, effect estimates, and 95% CIs.\nAnswer 25: We have replaced “restoration of functional pathways” with neutral wording. The text now reads: “After 2 years, dolutegravir-based regimens are associated with a gut microbiota composition more similar to that of HIV-negative individuals, compared with darunavir/ritonavir.”\nAnswer 26: We thank the reviewer for this clarification. Since reiterating priorization is indeed unnecessary, we have removed this statement from the abstract.\nAnswer 27: We thank the reviewer for this important clarification. We have revised the manuscript to ensure that our findings are not extrapolated beyond individuals with advanced HIV infection as defined in our cohort, and we explicitly acknowledge the broader clinical definition (<350 CD4 cells/µL). The text has been updated to reflect this distinction.\nAnswer 28: We thank the reviewer for this suggestion. We agree that “composition” is more precise in this context, and we have replaced the term accordingly.\nAnswer 29: We thank the reviewer for this helpful suggestion. We have revised the final paragraph of the Discussion to align more closely with the requested level of caution and to incorporate the suggested concluding sentence. The updated text now reflects the exploratory nature of the findings and avoids overinterpretation regarding clinical implications.\nAnswer 30: We acknowledge that the open-label design may introduce behavioral biases (e.g., diet, adherence). However, primary microbiome and biomarker endpoints were objectively measured in blinded laboratory analyses, minimizing the impact of participant or investigator awareness of treatment allocation.\nAnswer 31: We agree that the study population—predominantly Caucasian males from Spain—limits generalizability. We now explicitly acknowledge this limitation in the Discussion, emphasizing that further studies in more diverse global cohorts are needed to confirm these findings.\nAnswer 32: We recognize that dietary habits strongly influence the gut microbiota. Unfortunately, dietary data were not collected in this study, which we now state as a limitation. Nonetheless, the randomized design minimizes systematic dietary differences between treatment arms.\nAnswer 33: We agree that mucosal-associated microbiota may play a distinct role in gut barrier function. However, only fecal (luminal) samples were available for this cohort. This limitation is now acknowledged in the Discussion, with reference to future studies incorporating mucosal sampling.\nAnswer 34: We concur that no clinical endpoints (e.g., non-AIDS comorbidities) were measured. We now explicitly note in the Discussion that, although DTG was associated with reductions in inflammatory biomarkers, the clinical significance and long-term health implications of these findings remain to be established. In addition, after performing the updated analyses requested by the reviewers, oral pathobionts are no longer among the taxa enriched in the DTG arm. We have revised the text accordingly.\nAnswer 35: We thank the reviewer for raising this important point. Following the reviewer’s comment, we performed updated analyses to reassess taxa enriched in the DTG arm. In the revised dataset, oral commensals/pathobionts such as Streptococcus sanguinis and Rothia mucilaginosa are no longer detected among DTG-associated taxa.\nAnswer 36: No. All participants remained on their originally assigned treatment throughout the trial, and no switches to DRV/cobicistat occurred.\nAnswer 37: We thank the reviewer for this helpful suggestion. To improve transparency and strengthen the assessment of generalizability, we have added a new subsection in the Methods (“Participant exclusions”) providing a detailed description of the reasons for exclusion prior to randomization. This includes a breakdown of the medical decisions affecting 31 individuals, most commonly severe AIDS-defining conditions compromising clinical stability, as well as other medical and non-medical reasons.\nAnswer 38: We thank the reviewer for raising this point. Baseline IL-6 measurements were available for all participants, and there were no missing samples at baseline. Apparent gaps in the plots reflect values below the assay’s lower limit of detection rather than missing data. Importantly, values below the LLOD were still reported by the assay and included in the dataset.\nThe lower limit of detection for IL-6 was <2 pg/mL. We have clarified this point in the Methods section to avoid any confusion.\nAnswer 39: Thank you for this suggestion. In the initial version of the manuscript, longitudinal cytokine trajectories by treatment arm were presented in Supplementary Figure 1. We have now expanded this analysis, and the full set of cytokine trajectories can be found in the revised Supplementary Figure 1 and the new Supplementary Figure 2.\nAnswer 40: Thank you for the comment. As noted in our response to Reviewer 1, we have updated the analysis to linear mixed-effects models that include participant-level random effects to appropriately account for within-participant correlation.\nAnswer 41: Thank you for the suggestion. As clarified in our response to Reviewer 1, we implemented PERMANOVA with permutation constrained within participants (i.e., stratified by participant ID), thereby accounting for the repeated-measures structure and avoiding inflation of Type I error. We have now added this clarification to the Methods section.\nAnswer 42: Thank you for this helpful comment. As clarified in our response to Reviewer 1, all cytokine–taxa associations were estimated now using repeated-measures correlation (rmcorr), which explicitly accounts for within-participant dependence. Consistent with our analysis pipeline, multiple-testing corrections were applied per cytokine, by adjusting p-values across all taxa tested for each marker within each treatment arm.\nAnswer 43: We thank the reviewer for this suggestion. However, unlike PCA or PCoA, non-metric multidimensional scaling (NMDS) does not decompose variance along axes, so percentages of explained variance cannot be reported. Instead, the adequacy of the ordination is assessed by the stress value, which quantifies how well the reduced-dimensional configuration preserves the original dissimilarities. In our analysis, the 2D NMDS solution had a stress of 0.20, which is considered an acceptable representation of the data.\nAnswer 44: We have updated the analysis to apply a prevalence cutoff of 0.1, as suggested. In addition, taxon-level prevalence values are now reported, allowing readers to consider this factor when interpreting differential abundance results.\nAnswer 45: We have replaced the volcano plots with forest plots, which provide a clearer view of temporal changes in abundance and display effect sizes with 95% CIs rather than emphasizing significance. This avoids the visual ambiguity noted for baseline vs. non-significant taxa and aligns with Nature’s statistical reporting guidelines and the recommendations from Reviewer 1.\nAnswer 46: We have expanded the figure legends throughout the manuscript to include the distance metrics used, as well as clarifications on point size\nAnswer 47: Missing timepoints were handled using the mixed-effects modeling framework, which accommodates unbalanced longitudinal data without requiring imputation. Participants with missing visits were therefore retained, and all available observations contributed to the estimates. This clarification has been added to the Methods section.\nAnswer 48: We have added the SAP-defined endpoints and visit schedule to the Methods section so that these details are now available in the main text in addition to the protocol.\nAnswer 49: We thank the reviewer for this valuable suggestion. We agree that the limited number of female participants is an important limitation, given the known sex-related differences in metabolism, gut microbiota and immune responses. In line with the SAGER guidelines, we have added a statement in the Discussion acknowledging this limitation and emphasizing the importance of recruiting more balanced and diverse cohorts in future studies to strengthen generalizability.\nAnswer 50: We thank the reviewer for this insightful observation. In response, we quantified within-arm variability in longitudinal trajectories of α-diversity metrics by estimating subject-level slopes of change from baseline and comparing their dispersion between treatment arms. Variance ratios (DTG/DRV/r) and bootstrap confidence intervals were computed to assess differences in intra-group variability, complemented by Brown–Forsythe and Fligner–Killeen tests of equality of variances. These analyses, now presented in Supplementary Fig. 3A and described in the Methods section, show that DTG recipients displayed lower variability in slope estimates across participants, indicating more homogeneous increases in α-diversity over time. We have also incorporated a corresponding description in the Results section to contextualize this finding.\nAnswer 51: We thank the reviewer for this observation and agree that the detection of oral commensals in stool should be pointed out, as it requires careful interpretation. Therefore, we have added a paragraph in the Discussion highlighting possible explanations - including improved mucosal homeostasis, oral-gut translocation, impaired barrier function or low-level contamination, while emphasizing that our data cannot discriminate among these alternatives and should be interpreted with caution.\nAnswer 52: We thank the reviewer for the comment. We have revised the Discussion accordingly.\nAnswer 53: We thank the reviewer for this suggestion. Indeed, the microbiome shifts observed in our study could potentially contribute to metabolic reconfiguration, as suggested by the functional changes shown in Figure 6. We have expanded the Discussion to address this point in more detail.\nAnswer 54: We thank the reviewer for this constructive feedback. We agree that clear documentation is essential for reproducibility and broader community use. To address this, we have added a comprehensive README to the GitHub repository that includes environment setup, step-by-step instructions to run the analyses and information on data availability."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have done an excellent job in thoughtfully addressing the reviewer’s comments. In particular, the reanalysis using linear mixed-effects models and the replacement of multiple within-arm comparisons with model-based estimates of change from baseline represent a clear methodological improvement and substantially strengthen the longitudinal inference.\nFor Supplementary Figure 2, the legend does not fully explain how the mixed-effects model coefficients were transformed into the plotted Δ log₂ fold-changes. Please clarify (i) which outcomes were log₂-transformed prior to modeling, (ii) whether models were fitted on the log₂ scale or back-transformed post hoc, and (iii) how estimated fixed effects and contrasts were converted into Δ log₂ changes relative to baseline for visualization. Providing these details would improve transparency, reproducibility, and interpretability of an otherwise sound and well-reported analysis.", "reviewer_2": "I appreciate the authors’ comprehensive and thoughtful revisions. My previous concerns have been adequately addressed, particularly the clarification of the statistical methodology, removal of overstated claims regarding clinical implications, and reanalysis of taxa previously flagged as potential contaminants. I am satisfied with the current version and support publication.", "reviewer_3": "The authors have provided excellent rebuttal and appropriate changes in response to my prior reviewer comments. Thank you for the opportunity to review this excellent work."}, "author_response": "Point-by-point Response to Reviewers\nWe thank the reviewers for their positive assessment of our revised manuscript and their constructive feedback throughout the review process. We are pleased that the reanalysis using linear mixed-effects models has strengthened the study’s conclusions.\nReviewer #1\n1. The authors have done an excellent job in thoughtfully addressing the reviewer’s comments. In particular, the reanalysis using linear mixed-effects models and the replacement of multiple within-arm comparisons with model-based estimates of change from baseline represent a clear methodological improvement and substantially strengthen the longitudinal inference.\nResponse: We are grateful for the reviewer’s appreciation of the methodological improvements and the reanalysis performed.\n2. For Supplementary Figure 2, the legend does not fully explain how the mixed-effects model coefficients were transformed into the plotted Δ log₂ fold-changes. Please clarify (i) which outcomes were log₂-transformed prior to modeling, (ii) whether models were fitted on the log₂ scale or back-transformed post hoc, and (iii) how estimated fixed effects and contrasts were converted into Δ log₂ changes relative to baseline for visualization. Providing these details would improve transparency, reproducibility, and interpretability of an otherwise sound and well-reported analysis.\nResponse: We thank the reviewer for pointing this out. We have revised the legend of Supplementary Figure 2 to include the requested details. Specifically, we have clarified that:\n- The outcomes were log₂-transformed prior to modeling to stabilize variance.\n- The linear mixed-effects models were fitted directly on the log₂ scale.\n- The plotted values represent the estimated marginal means (EMMs) of the fixed effects (interaction terms) derived directly from the model, representing the Δ log₂ change relative to the baseline without post-hoc back-transformation.\nReviewer #2\n1. I appreciate the authors’ comprehensive and thoughtful revisions. My previous concerns have been adequately addressed, particularly the clarification of the statistical methodology, removal of overstated claims regarding clinical implications, and reanalysis of taxa previously flagged as potential contaminants. I am satisfied with the current version and support publication.\nResponse: We thank the reviewer for their time, their constructive feedback during the previous round of revision, and their support for the publication of this work.\nReviewer #3\n1. The authors have provided excellent rebuttal and appropriate changes in response to my prior reviewer comments. Thank you for the opportunity to review this excellent work.\nResponse: We sincerely thank the reviewer for their positive evaluation and their helpful comments throughout the review process.\n2. Comment (on code availability): The authors have made the code and README files available.\nResponse: We thank the reviewer for confirming the availability of the code and documentation."}]

Decision letter