COVID-19, caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has resulted in over 780 million confirmed cases and 7 million deaths worldwide since its emergence in 2020 [1–3]. While the infection primarily targets the upper respiratory tract and often causes mild illness, it can progress to severe respiratory failure, requiring hospitalisation, and may result in death [3,4]. Approximately 10–20% of COVID-19 patients develop post-acute COVID-19 syndrome, also known as long COVID, characterised by persistent symptoms lasting beyond four weeks from symptoms onset without alternative diagnoses [5]. These symptoms affect multiple organ systems, including respiratory (cough, dyspnoea), cardiovascular (chest pain, palpitations), neuropsychiatric (anxiety, insomnia), musculoskeletal (joint pain, muscle weakness), dermatological (rashes, hair loss), and gastrointestinal (loss of appetite, diarrhoea) manifestations [5,6].
While several studies have investigated the impact of SARS-CoV-2 infection on the gut microbiome, findings regarding microbial diversity and taxonomic alterations have been inconsistent [3,7–10]. Gastrointestinal symptoms such as diarrhoea and nausea have been reported in 2–10% of COVID-19 patients [11,12]. The extrapulmonary dissemination of the virus is likely linked to the expression of its cellular receptor, angiotensin-converting enzyme II, in various susceptible human cells, including the small intestine and colon [8,13–15]. Therefore, the interaction between SARS-CoV-2 and angiotensin-converting enzyme II may disrupt gut microbial ecology and compromise intestinal homeostasis [3,8].
Currently, clinical management of long COVID is symptom-based, focusing primarily on rehabilitation and pharmacological interventions to alleviate symptoms [7,16]. Drugs targeting organ-specific syndromes approved for use by the Food and Drug Administration have been evaluated in early-phase proof-of-concept trials as a safe approach for long COVID within experimental clinical studies [17]. For instance, montelukast, an orally administered medication approved for asthma, has potential for alleviating long COVID-related respiratory symptoms due to its anti-inflammatory properties and efficacy in reducing oxidative stress in acute asthma [18,19]. However, commonly used non-antibiotic drugs may contribute to microbial gut dysbiosis by promoting or inhibiting the growth of specific bacterial taxa, thereby altering both microbial abundance and function [20,21]. Given that individuals with distinct gut microbiome profiles may exhibit different drug responses, we aimed to compare the gut microbiota composition of patients with long COVID and healthy controls, and to evaluate the effect of montelukast on the gut microbiota composition of patients with long COVID.
METHODS
Study design and sample collection
We conducted this proof-of-concept observational study at Hospital Universitari de Bellvitge using a cohort of 45 long COVID patients with persistent mild-moderate dyspnoea attended in primary healthcare without previous hospital admissions for SARS-CoV-2. These participants were previously enrolled in the clinical trial E-SPERANZA (ClinicalTrials.gov registration number: NCT04695704) [22], led by Institut Universitari d’Investigació en Atenció Primària Jordi Gol, in which they were randomised to receive treatment either with montelukast or placebo. In an independent follow-up study to investigate gut microbiome dynamics, stool samples were collected from this cohort at three time points: baseline before treatment (visit 1), after completing the 30-day treatment (visit 2), and at one year after baseline (visit 3). Patients were recruited through their primary care centres following a clinical diagnosis of long COVID, defined as a confirmed SARS-CoV-2 infection (positive RT-PCR, antigenic test, or equivalent <10 days from onset) managed in primary care, with persistent mild-to-moderate dyspnoea for more than four weeks and less than one year. All participants provided written informed consent. At enrolment, participants received a stool DNA collection & stabilisation kit (Canvax) with detailed self-collection instructions. Samples were stabilised immediately in DNA preservation buffer, transported to the microbiology laboratory and stored at -80°C until processing. A total of 106 stool samples were collected from long COVID patients (Figure S1 in the Online Supplementary Document): 45 samples at baseline (visit 1), 32 after treatment (visit 2), and 29 after one year (visit 3). Of the 45 participants, 23 received montelukast and 22 received placebo. We also included 45 stool samples from a previously published healthy control cohort [23], matched by age, sex and sample size (Table S1 in the Online Supplementary Document). Both cohorts were processed at the same institution, ensuring standardised protocols. The healthy control group, consisting of health care workers, was enrolled between November 2022 and April 2023, providing a contemporaneous geographical and temporal match to the long COVID cohort (recruited from August 2021 to March 2023). This overlap minimises potential confounding factors related to environmental changes or technical variability. The study involved two analytical approaches: a cross-sectional comparison of baseline gut microbiome composition between patients and healthy controls; and a longitudinal assessment of gut microbiome dynamics following montelukast treatment or placebo.
DNA extraction and 16S rRNA gene amplification and sequencing
Total DNA was extracted from stool samples using the KingFisher Flex system (Thermo Fisher Scientific) and the MagMAXTM Microbiome Ultra Nucleic Acid Isolation kit, which incorporates mechanical disruption with zirconia beads, following the manufacturer’s instructions. DNA concentrations were quantified using a Qubit Flex fluorometer (Invitrogen). The V3-V4 regions of the 16S rRNA gene were amplified following the Illumina’s guidelines using the following primers: forward (5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG-3′) and reverse (5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC-3′). Amplicon size and integrity were verified on a TapeStation 4200 system (Agilent Technologies) using the D1000 ScreenTape assay. Libraries were quantified with Qubit (Invitrogen), pooled in equimolar amounts using a Sciclone G3 NGSx Workstation (PerkinElmer) and sequenced on an Illumina MiSeq platform (2×300).
Bioinformatic and statistical analysis for 16S rRNA gene sequencing
Raw 16S rRNA FASTQ files were processed in QIIME2 [24] for quality control, including removal of index sequences and trimming of low-quality bases. Forward and reverse reads were merged, denoised, and clustered into amplicon sequence variants using DADA2 [25]. Taxonomic assignment was performed against the SILVA 138 reference database [26]. Amplicon sequence variants were further filtered using Phyloseq (R-Studio) [27] to remove potential contaminants (Eukaryota, Archaea, mitochondria, and chloroplasts) and aggregated at the genus level. Sequences were rarefied before alpha diversity analysis to account for uneven sequencing depth.
We verified statistical assumptions, including normality and homogeneity of variance, using the Shapiro–Wilk and Levene’s tests, respectively. For the cross-sectional analysis, we used Student’s t-test for parametric data and the Wilcoxon rank-sum used for non-parametric distributions. For longitudinal comparisons, we performed paired analyses to evaluate changes within the same individual over time, using the paired Student’s t-test or the Wilcoxon signed-rank test, as appropriate. To account for multiple comparisons in these longitudinal assessments, we adjusted P-values using the Bonferroni correction.
We assessed alpha diversity using the Shannon index. For beta diversity, we normalised microbial abundances to compositional data via total sum scaling transformation, and we evaluated community structure using Bray–Curtis, visualised via principal coordinate analysis and further analysed by permutational multivariate analysis of variance. We calculated relative abundances at all taxonomic levels, and analysed differentially abundant taxa using the linear discriminant effect size analysis (LEfSe) method [28]. We set the significance level at 0.05 (two-sided) for all statistical tests.
RESULTS
Long COVID characterisation: cross-sectional analysis
We performed a cross-sectional comparison between 45 baseline stool samples and 45 age- and sex-matched healthy control samples to investigate potential alterations in the gut microbiome of long COVID patients (Table S1 in the Online Supplementary Document). While alpha diversity analysis revealed no significant differences in the Shannon index, indicating comparable microbial richness and evenness (Figure 1, Panel A), beta diversity analysis based on Bray–Curtis dissimilarity showed a clear separation between groups (Figure 1, Panel B). Statistical significance was confirmed by PERMANOVA (ADONIS) testing, which revealed that the gut microbial composition of long COVID patients differed significantly from that of healthy controls (R2 = 0.081, P-value = 0.001). This indicates that disease status explains 8.1% of the total microbial community variation. This finding was further supported by a test of multivariate homogeneity of group dispersions, which showed no significant differences between groups (P-value >0.05), confirming that the observed separation is driven by a true compositional shift rather than differences in group dispersion.
Figure 1. Cross-sectional analysis of the gut microbiome, relative microbial abundances and differentially enriched taxa between long COVID patients and healthy subjects. Comparison includes baseline samples from long COVID patients (n = 45) and healthy controls (n = 45). Panel A. Alpha diversity (Shannon index boxplots). Panel B. Beta diversity (principal coordinate analysis plot based on Bray–Curtis distance). Relative abundance at the phylum level (Panel C) and genus level (Panel D) in long COVID patients and healthy controls. Linear discriminant effect size analysis at the phylum level (Panel E) and genus level (Panel F) showing taxa with significant differential abundance based on the linear discriminant analysis score.
At the phylum level, both groups presented a typical gut microbiota dominated by Firmicutes and Bacteroidota, with smaller proportions of Actinobacteria, Proteobacteria, and Verrucomicrobiota (Figure 1, Panel C). Discriminant analysis revealed compositional differences between groups. The LEfSe identified Verrucomicrobiota and Actinobacteriota as significantly enriched in healthy controls, whereas Firmicutes were more abundant in long COVID patients (Figure 1, Panel E). Genus-level profiling provided further resolution of these patterns. The ten most prevalent genera across all samples included Bacteroides, followed by Faecalibacterium and Blautia as dominant taxa (Figure 1, Panel D). The LEfSe further identified taxa contributing to intergroup variability: Alistipes and Akkermansia were significantly enriched in healthy controls, while Agathobacter and Faecalibacterium were more abundant in long COVID patients (Figure 1, Panel F).
Longitudinal long COVID analysis
Since long COVID patients were randomly assigned to receive either placebo or montelukast, we performed a longitudinal characterisation of the gut microbiome dynamics to explore whether treatment induced compositional or structural changes over time.
Alpha diversity analyses showed no significant temporal variation in richness and evenness, indicating that microbial diversity remained stable throughout the study period in both placebo and montelukast groups (Figure 2). Similarly, beta diversity analysis showed no significant longitudinal changes within groups, suggesting preservation of the core microbial community structure regardless of treatment (Figure S2 in the Online Supplementary Document).
Figure 2. Longitudinal analysis of alpha diversity in long COVID patients treated with montelukast or placebo. Panel A. Alpha diversity (Shannon index) at visit 1 (baseline), visit 2 (30 days), and visit 3 (one year) in the placebo and montelukast groups. Panel B. Interaction plot showing the longitudinal trend of alpha diversity measures between treatment arms across the three time points.
At the compositional level, phylum relative abundances remained largely stable across visits in both study arms. However, at the genus level, subtle temporal trends were observed among the most abundant taxa. Bacteroides dominated both groups, particularly the placebo group, where its relative abundance progressively increased from V1 to V3, accompanied by a concomitant decrease in Prevotella. In contrast, the montelukast group did not follow this trajectory, showing a transient decrease in Bacteroides abundance at V2 that interrupted the upward trend seen in the placebo group (Figure 3).
Figure 3. Longitudinal relative abundance of the top 10 genera in long COVID patients. Upper panels represent the placebo group and lower panels represent the montelukast group across the three study time points: visit 1 (baseline), visit 2 (30 days), and visit 3 (one year).
Looking into detail at these two genera, although no statistically significant differences in relative abundance of Bacteroides and Prevotella were observed between groups at any time point, visual observation of taxonomic temporal patterns highlighted subtle, yet consistent tendencies (Figure S3 in the Online Supplementary Document). In the placebo group, Bacteroides showed a continuous increase and Prevotella showed a progressive decrease between baseline and the follow-up visits, suggesting a stable directional shift in these dominant genera. In contrast, the montelukast group displayed less consistent dynamics, characterised by minor fluctuations over time, including a temporary increase in Prevotella at V2. To note, the high variability of Prevotella in the placebo group compared to the montelukast-treated group.
To further assess the effect of montelukast, we compared gut microbiome composition between groups at visit 2, corresponding to the end of the 30-day treatment. LEfSe identified Dialister as the only genus differentiating both groups, showing higher relative abundance in the montelukast-treated patients (Figure S4 in the Online Supplementary Document). No other taxa were significantly enriched in either group, and alpha and beta diversity comparisons at visit 2 revealed no significant differences, indicating that overall microbial diversity and community structure remained comparable between treatment groups.
DISCUSSION
The COVID-19 pandemic posed unprecedented global health challenges, not only during the acute phase of infection, but also through the long-term sequelae experienced by many patients with long COVID [5–7]. This clinical condition prompted an urgent search for therapeutic strategies to relieve persistent symptoms and improve quality of life [7,16,17]. While some interventions have shown partial or context-dependent benefits, others have yielded reduced efficacy. In this context, growing interest has focused on the gut microbiome as a potential modulator of persistent symptoms, although current evidence remains heterogeneous and sometimes contradictory. This study contributes to the growing body of evidence on gut microbiome alterations in long COVID by first establishing a baseline comparison with a healthy population to characterise the long COVID microbial signature, and subsequently evaluating the impact of short-term montelukast treatment on these microbial profiles.
Our findings support the growing body of evidence indicating that the gut microbiome is altered in long COVID patients [3,8–10]. In our cohort, this condition was characterised by a distinct restructuring of the microbial community composition. Specifically, our analysis revealed significant shifts in beta diversity that clearly separate long COVID patients from healthy controls, with disease status explaining 8.1% of the total community variance (R2 = 0.081). This clustering is driven by substantial differences in taxonomic profiles rather than internal group variability, as confirmed by the homogeneity of dispersions, supporting the identification of a specific microbial signature associated with long COVID in these patients. Furthermore, we observed that alpha diversity remained stable, indicating that overall microbial richness was not depleted. While SARS-CoV-2 interacts with gastrointestinal epithelial cells, these interactions do not invariably result in substantial loss of taxa, as is often reported in acute SARS-CoV-2 infection [29]; instead, these results imply that the observed dysbiosis in long COVID is qualitative rather than quantitative.
This realignment was also evident at the taxonomic level. At the phylum level, both groups had a typical gut microbiota, which was dominated by Firmicutes and Bacteroidota [30–32]. In our cohort, the three most prevalent genera were Bacteroides, Faecalibacterium, and Blautia. This reflects a core community structure that is preserved in long COVID patients [30,33]. In contrast, genera such as Clostridium, Eubacterium, and Ruminococcus were less abundant across all samples, possibly due to interindividual variability and external influences, such as diet, medication, or lifestyle [30–32,34,35]. Despite these commonalities, LEfSe discriminant analysis identified specific taxa contributing to the observed intergroup differences. Healthy controls were characterised by an enrichment of the phyla Verrucomicrobiota and Actinobacteriota, as well as the genera Akkermansia and Alistipes. Conversely, long COVID patients showed a higher abundance of Firmicutes, driven by an enrichment of genera such as Agathobacter and Faecalibacterium. These findings reinforce the idea that long COVID microbial restructuring is not defined by a generalised loss of key commensals but rather by a systemic reconfiguration of the microbial community.
Longitudinal analysis further confirmed overall stability of gut microbial diversity over time. The inverse relationship between Bacteroides and Prevotella, commonly associated with metabolic and inflammatory phenotypes [35,36], was also reflected in our cohort, suggesting group-specific trajectories after treatment. The temporal dynamics of the Bacteroides/Prevotella association diverged between groups. In the placebo group, Bacteroides abundance progressively increased while Prevotella declined, reflecting a gradual shift toward a Bacteroides-enriched state. In contrast, the montelukast group interrupted this trajectory by inducing a transient reduction in Bacteroides after the 30-day treatment. This Bacteroides reduction pattern resembles findings from both culture-based and metagenomic studies of antimicrobial exposure, where short-term disruptions to the gut microbiota, with a decline in Bacteroides alongside diversity loss are commonly observed [37]. However, in our study, montelukast caused a reduction in Bacteroides without a broad disruption of the gut microbial community. These changes are consistent with a drug-associated modulation of community dynamics rather than a simple time-dependent drift, aligning with evidence for taxon-specific effects of non-antibiotic compounds [20,21].
The LEfSe analysis combining cross-sectional and intervention data identified several differentially abundant taxa. When comparing long COVID patients vs. healthy controls, Alistipes and Akkermansia were enriched in healthy controls, while Agathobacter and Faecalibacterium were more abundant in long COVID patients. This pattern partially contrasts with acute/post-acute reports of short-chain fatty acids-producers depletion. Akkermansia muciniphila, a mucin-degrader promoting gut barrier function, is generally linked to favourable inflammatory and metabolic profiles [38] and Alistipes species exhibit context-dependent roles, with higher abundances previously associated with favourable acute COVID-19 recovery trajectories [39]. In contrast, Faecalibacterium, a major butyrate producer with anti-inflammatory activity [40], shows species-specific behaviours. For instance, F. prausnitzii enrichment has been connected to neurological long COVID symptoms such as memory loss, poor concentration, and insomnia [41]. When examining montelukast-treated vs placebo groups, Dialister was the only genus enriched in the treatment arm. This genus has been associated with pro-inflammatory conditions [42] potentially reflecting a transient inflammatory or immunomodulatory response rather than a sustained microbiome alteration. Taken together, the LEfSe findings suggest that taxon-specific fluctuations, rather than broad compositional shifts, occur both across disease states and following short-term drug exposure.
However, these results should be interpreted with caution, and several limitations should be acknowledged. Given the proof-of-concept nature of this study, the modest sample size limits statistical power, potentially obscuring the detection of subgroup-specific effects, particularly at lower taxonomic ranks. Clinical heterogeneity, including prior treatmentsand diet and lifestyle factors, may further contribute to variability. Although we applied specific exclusion criteria for major chronic comorbidities to ensure a more homogeneous cohort, we acknowledge that the lack of detailed metadata remains a potential source of bias. Additionally, the short intervention period captures only early microbial responses to montelukast, possibly missing delayed or cumulative effects. The lack of baseline samples prior to long COVID prevents confirming whether treatment or time promotes trends toward microbiome reversion.
CONCLUSIONS
These results suggest that long COVID involves significant structural reorganisation of the gut ecosystem rather than depletion of its fundamental components. The condition appears to be marked by a compositional shift and quantitative re-equilibration of key microbial taxa, most notably within the phylum Firmicutes, specifically the Faecalibacterium and Agathobacter genera, indicating an altered intestinal homeostatic balance even when overall species richness is maintained. These findings suggest a qualitative transformation of the microbiota, where the core community’s identity remains present, but its internal proportions are disrupted. Furthermore, short-term montelukast administration seems to exert transient effects at the genus level without causing broader microbial landscape destabilisation. These targeted changes highlight the complexity of gastrointestinal involvement in long COVID, where the microbial architecture undergoes a persistent realignment. Extended follow-up and integrative functional analysis are essential to determine if these shifts influence host-microbiome interactions and long-term recovery.