Altered oral and gut microbiota and its association with SARS-CoV-2 viral load in COVID-19 patients during hospitalization

Abstract

The human oral and gut commensal microbes play vital roles in the development and maintenance of immune homeostasis, while its association with susceptibility and severity of SARS-CoV-2 infection is barely understood. In this study, we investigated the dynamics of the oral and intestinal flora before and after the clearance of SARS-CoV-2 in 53 COVID-19 patients, and then examined their microbiome alterations in comparison to 76 healthy individuals. A total of 140 throat swab samples and 81 fecal samples from these COVID-19 patients during hospitalization, and 44 throat swab samples and 32 fecal samples from sex and age-matched healthy individuals were collected and then subjected to 16S rRNA sequencing and viral load inspection. We found that SARS-CoV-2 infection was associated with alterations of the microbiome community in patients as indicated by both alpha and beta diversity indexes. Several bacterial taxa were identified related to SARS-CoV-2 infection, wherein elevated Granulicatella and Rothia mucilaginosa were found in both oral and gut microbiome. The SARS-CoV-2 viral load in those samples was also calculated to identify potential dynamics between COVID-19 and the microbiome. These findings provide a meaningful baseline for microbes in the digestive tract of COVID-19 patients and will shed light on new dimensions for disease pathophysiology, potential microbial biomarkers, and treatment strategies for COVID-19.

Fig. 1. The study design and variations of the oral and gut microbiota in COVID-19 patients.

a Graphic representation of study design and sample collection. PT(+), throat swab samples from COVID-19 patients during the positive viral RNA test period; PT(−), throat swab samples from COVID-19 patients during the negative viral RNA test period; PF(+), fecal specimens from COVID-19 patients during the positive viral RNA test period; PF(−), fecal specimens from COVID-19 patients during the negative viral RNA test period. b Bacterial composition at the phylum level among groups. Sample size (n) for each group: HT, n = 44; PT(+), n = 52; PT(−), n = 88; HF, n = 32; PF(+), n = 50; PF(−), n = 31. HT, throat swab samples from the healthy controls; HF, fecal specimens from the healthy controls. c PCoA plot based on the unweighted UniFrac distance depicting differences in the bacterial community among groups.

Fig. 2. Altered oral microbiota in COVID-19 patients.

a Species diversity differences among the HT (n = 44), PT(NS) (n = 90), and PT(S) (n = 50) groups were estimated by Faith’s phylogenetic diversity index. *p < 0.05, ****p < 0.0001. p-values were obtained using one-way ANOVA followed by Tukey’s multiple comparisons test. The line in the middle of the box, bound of the box and whiskers represent the median, 25th–75th percentiles, and min-to-max values, respectively. PT(NS), throat swab samples from the non-severe patients group; PT(S), throat swab samples from the severe patients group. b PCoA plot based on unweighted UniFrac distances showing microbial structural clustering among the HT, PT(NS), and PT(S) groups. Ellipses represent 68% confidence intervals. c PCoA plot according to unweighted UniFrac distances displaying bacterial structural discrimination among the HT (n = 44), PT(abx−) (n = 49), and PT(abx+) (n = 91) groups. Ellipses represent 68% confidence intervals. PT(abx−), throat swab samples from patients without antibiotic interference; PT(abx+), throat swab samples from patients with antibiotic interference. d Bacterial taxa identified a significant difference between HT (n = 44) and PT (n = 140) groups by LEfSe. Seventeen taxa were validated using MaAsLin2 adjusting for covariates. PT throat swab samples from patients, LDA linear discriminant analysis. e Differentially abundant MetaCyc pathways identified with the functional analysis result of metagenomes. Only the pathways identified by both LEfSe and MaAsLin2 were presented. Circle size represents the relative abundance of pathways.

Fig. 3. Altered gut microbiota in COVID-19 patients.

a Species diversity differences among the HF (n = 32), PF(NS) (n = 64), and PF(S) (n = 17) groups were estimated by Faith’s phylogenetic diversity index. ns: not significant, *p < 0.05, ****p < 0.0001. p-values were obtained using one-way ANOVA followed by Tukey’s multiple comparisons test. The line in the middle of the box, bound of the box and whiskers represent the median, 25th–75th percentiles, and min-to-max values, respectively. PF(NS), fecal specimens from the non-severe patient group; PF(S), fecal specimens from the severe patient group. b PCoA plot based on unweighted UniFrac distances showing microbial structural clustering among the HF, PF(NS), and PF(S) groups. Ellipses represent 68% confidence intervals. c PCoA plot according to unweighted UniFrac distances displaying bacterial structural discrimination among the HF (n = 32), PF(abx−) (n = 30), and PF(abx+) (n = 51) groups. Ellipses represent 68% confidence intervals. PF(abx−), fecal specimens from patients without antibiotic interference; PF(abx+), fecal specimens from patients with antibiotic interference. d Bacterial taxa identified a significant difference between the HF (n = 32) and PF (n = 81) groups by LEfSe. Seventeen taxa were validated using MaAsLin adjusting for covariates. PF fecal samples from patients, LDA linear discriminant analysis. e Differentially abundant MetaCyc pathways identified with the functional analysis result of metagenomes. Only the pathways identified by both LEfSe and MaAsLin2 were presented. Circle size represents the relative abundance of pathways.

Fig. 4. Associations between oral/gut microbiota disturbance and SARS-CoV-2 viral loads in hospitalized patients.

a Longitudinal changes of viral loads in throat swab samples and fecal specimens from COVID-19 patients. Day 0 indicated the negative conversion time of SARS-CoV-2 RNA. The red line and blue line indicated longitudinal changes in viral loads represented by the abundances of the E gene and N gene, respectively. b Boxplots showing unweighted UniFrac distances within and between groups. c The longitudinal viral load represented by the abundances of the E gene significantly correlated with the timepoint-matched relative abundances of different taxa in throat swab and fecal samples separately. p-values were obtained using CCLasso.