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

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.


INTRODUCTION
The ongoing COVID-19 pandemic caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has jeopardized global public health and safety. As of June 21, 2021, the COVID-19 pandemic has infected more than 179 million people, resulting in over 3.8 million causalities globally. Transmission of SARS-CoV-2 has commonly occurred via breathing or direct contact with viruscontaining droplets and aerosols from infected people through coughing and sneezing. When the SARS-CoV-2 particle reaches the nasal cavity of the host, it enters the host epithelial cells through the angiotensin-converting enzyme 2 (ACE2) receptor, which is prominently presented on epithelial cells lining in the respiratory and digestive tract systems 1,2 . Most infections lead to a prompt innate immune response and the virus gets eradicated or controlled quickly, manifesting none or mild symptoms. Whereas in some patients, viruses located in the upper airway replicate further toward the lower respiratory tract to activate enhanced pro-inflammatory responses, therefore resulting in severe outcomes including acute respiratory syndrome, organ malfunctioning, and even death 3 . Not limited to the respiratory tract, SARS-CoV-2 is known to target different organ systems 4 . For instance, around 55% of patients were observed with prolonged viral RNA present in the feces, even weeks after the viral clearance in their respiratory tract 5 . It has been reported that abnormal immune response, comorbid conditions, and advanced age are risk factors linked to COVID-19 severity. However, these factors are insufficient to offer a satisfactory explanation of all patients' severe disease outcomes. As global mass vaccination and proven effective therapy for COVID-19 remain unclear, explorative efforts towards new perspectives of protection and therapeutic approaches for COVID-19 are indispensable 6,7 .
The human microbiome is important for developing and maintaining immune homeostasis and it is known that microbiota imbalance or dysbiosis are highly associated with various diseases. The intestinal tract and oral cavity, with the largest and secondlargest microbiota in the human body, play significant roles in the pathogenesis of infectious disease. Previous studies have reported that oral-lung microbes can influence the outcome of many infectious diseases by regulating the host mucosal immunity [8][9][10] . Intestinal flora can affect the occurrence and progression of viral infection through the gut-lung axis [11][12][13] . Conversely, the indigenous microbiome could be disturbed by a viral infection, leading to alterations in susceptibility and disease severity through dysbiotic community structure and function 14 . The unequivocal association between influenza viral and bacterial co-infection and disease severity has been proved in early studies 15 . Likewise, evidence has suggested that SARS-CoV-2 infection could predispose patients to bacterial co-infections and superinfections, resulting in increased disease severity and mortality 14 . Further, the dysbiosis of influenza-infected individuals progresses toward microbiota homeostasis, coinciding with viral clearance and recovery, suggesting that the health status of microbial flora is likely a fine indication of disease recovery 16 . Zuo and colleagues have also suggested that fecal microbiomes from COVID-19 patients were characterized by proliferation of opportunistic pathogens and depletion of favorable commensals compared to healthy controls with the shot-gun metagenomics approach 17 ; however, this study was limited by its size of subjects used. Moreover, validation of the profiles at the genus level is currently suggested to be performed through 16S rRNA gene microbial profiling to evaluate the presence of undetected microbes in the marker gene-based profiles 18 . Accumulated evidence for the oral-gut axis has revealed its role in modulating the pathogenesis process in numerous diseases 19,20 . It is compelling to look into the oral and intestinal microbiome combined with SARS-CoV-2 infection and the crosstalk among them, which may provide an improved understanding of the initiation of viral infection and the path of disease deterioration.

RESULTS
Overview of microbial composition in the study subjects A total of 53 patients diagnosed with COVID-19 and 76 healthy individuals were included in this study. As indicated in Fig. 1a To characterize the microbiotas in those aforementioned samples, we sequenced the V3-V4 region of the bacterial 16S rRNA gene. As expected, our data confirmed a distinct microbiota composition between oral and gut samples 21 . The three most abundant phyla Firmicutes, Proteobacteria, and Bacteroidetes accounted for 79.3% of the community in the oral samples. The two dominant phyla Firmicutes and Bacteroidetes accounted for 83.7% of the bacterial community in the gut samples (Fig. 1b). The difference was further illustrated by the principal coordinates analysis (PCoA) plot and Permutational multivariate analysis of variance (PERMANOVA) with unweighted UniFrac distance calculated at the sequence feature level. The microbiota composition of the oral and gut samples diverged from each other along the first axis in which 31.87% of the total variance was explained, indicating that the largest source of variation was the sample types ( Fig. 1c). Similarly, the PERMANOVA result indicated 30.80% of the total variance was explained by the sample types (oral samples vs. gut samples, PERMANOVA, R 2 = 0.308, p < 0.001). More importantly, samples from patients and healthy individuals seemed to possess strong dissimilarities for both oral (PERMA-NOVA, R 2 = 0.175, p < 0.001) and gut microbiome communities (PERMANOVA, R 2 = 0.196, p < 0.001), suggesting overwhelming shifts of microbiome structure were gained after SARS-CoV-2 infection in both oral cavity and feces. Meanwhile, the microbial    Table 3). Bacterial taxa including Treponema, Aggregatibacter, and P. intermedia were identified further depleted in the PT(S) group, implying their associations with disease severity. The comparison between the controls and the PT(abx-) group was carried out to assess the impact of infection on bacteria by excluding the antibiotics interference (Supplementary Table 3). Lastly, the potential functional consequence resulted from the drastic taxa compositional alteration was evaluated. The global pattern from Procrustes analysis displayed a good-fit correlation between the oral microbial community composition and microbial function (Mantel, rho = 0.403, p = 0.001, Supplementary Fig. 1a). The significant difference in functional pathways of the oral microbial communities between COVID-19 patients and the controls was shown in Fig. 2e. Notably, the top four depleted pathways in patients were involved in the TCA cycle used by all aerobic organisms to generate energy, indicating a disturbed microbial community.
Alterations of the gut microbiota in patients Next, we explored whether dysbiosis occurred in the gut microbiome. Faith's phylogenetic diversity index showed that its diversity was significantly lower in patients with the non-   Table 4). More interestingly, two normal components of the respiratory tract flora, Granulicatella and R. mucilaginosa, were significantly increased in the patients for both oral and gut samples, and R. mucilaginosa seemed further associated with the disease severity in the fecal samples ( Supplementary Fig. 2). Additional analysis focusing on the bacteria changes associated with disease severity and antibiotics usage in fecal samples were described in Supplementary Table 5. Notably, no difference in the bacterial taxa was identified between the PF(NS) group and the PF(S) group. Similar to the oral microbiome, the global pattern from Procrustes analysis displayed a good-fit correlation between the microbial composition and predicted function profile in the gut microbiota community (Mantel, rho = 0.268, p = 0.001, Supplementary  Fig. 1b). Also, the major distinct functional pathways associated with SARS-CoV-2 infection were shown in Fig. 3e, including two depleted pyrimidine deoxyribonucleotides biosynthesis pathways involved in viral suppression through innate immunity in patients 22 .

Relationship between SARS-CoV-2 virus and bacterial species
The viral loads in the throat swab and fecal samples spanning the full collection period were compared. The median time span of the throat swabs collection was 10 (IQR: 6-12) days for P-VRTP and 10 (IQR: 3-18) days for N-VRTP. The median time span of the fecal samples collection was 8 (IQR: 3.75-14) days for P-VRTP and 8 (IQR: 2-12) days for N-VRTP. The viral copy numbers derived from the E gene showed a relatively higher abundance over the numbers derived from the N gene (Fig. 4a). As E-gene and N-gene copy numbers showed a strong correlation therefore both can accurately reflect the abundance of viral copy (throat swab samples: Spearman rho = 0.872, p < 0.001; fecal samples: Spearman rho = 0.923, p < 0.001). The overall viral loads in throat swab samples seemed equally abundant with fecal specimens collected from patients during the P-VRTP, which is consistent with the previous studies 23 . Based on unweighted UniFrac distance, it was found that the differences in beta diversity between the HT and PT (+) groups (PERMANOVA, R 2 = 0.234, p = 0.001), the HT and PT(−) groups (PERMANOVA, R 2 = 0.218, p = 0.001), were significant from each other. Likewise, unweighted UniFrac dissimilarity revealed that the microbiomes of the PF(+) group (PERMANOVA, R 2 = 0.233, p = 0.001) and the PF(−) group (PERMANOVA, R 2 = 0.250, p = 0.001) were significantly different from that of the healthy controls. During the N-VRTP, we observed a microbial community recovery towards the healthy controls, that is, unweighted UniFrac distance between PT(−) vs. HT was lower than that between PT(+) vs. HT (Fig. 4b). Compared to the oral microbiome, this recovery trend in the microbiome in the gut was less evident, which may be attributed to a relatively shorter time span of sample collection and microbiome robustness in the fecal samples. Additionally, we investigated whether the SARS-CoV-2 viral load was associated with any oral and gut bacterial species. In the oral microbiota, Pelomonas was identified as positively correlated with the viral load of SARS-CoV-2 (Fig. 4c). While in the gut microbiota, P. copri and E. dolichum were identified positively correlated with the viral load of SARS-CoV-2, and S. anginosus, Dialister, Alistipes, Ruminococcus, C. citroniae, Bifidobacterium, Haemophilus, and H. parainfluenzae were identified negatively correlated with the viral load of SARS-CoV-2.

DISCUSSION
The oral and gut microbiome offers a range of valuable properties to the host. Several most significant contributions of these microbes are to boost metabolism, improve digestive health and strengthen resistance against pathogens. Furthermore, the complex crosstalk between commensal microbes and different body systems is essential for the functioning of the immune system 24 .
Our data revealed profound alterations in both oral and gut microbiomes, which were reflected in the dramatic changes in community structure, potential bacterial marker species, and predicted functional profile. The decline in commensal bacterial diversity has been considered as a key dysbacteriosis indicator in several diseases 25 . Coherently, the oral and gut microbiome of COVID-19 patients in our study exhibited decreased diversities compared to the healthy controls. This trend of microbial imbalance in patients with severe conditions compared to nonsevere patients was also observed in both sample types, though the alpha diversity was not significant (ANOVA, p = 0.73) in the fecal samples. Moreover, the microbiomes from the severe and the non-severe patients were partitioned into two clusters in both microbial populations. Overall, diversity results of impaired microbiota suggested the strong association between the microbiome community complexity and the disease severity in COVID-19 patients. Though experimental validation is needed, our result has highlighted the possibility of personalized microbiota to affect the disease outcome of COVID-19 patients 26 .
Numerous studies have reported high occurrences of bacterial co-infection in hospitalized COVID-19 patients, and the odds of the bacterial infection get even higher for ICU patients 27,28 . Among the list of eight bacterial taxa with enlarged relative abundance in patients in Fig. 2d, Veillonella, Campylobacter, R. mucilaginosa, Granulicatella, Kingella, and Filifactor belongs to a group of periodontitis-correlated taxa, adding evidence to support the previous work denotative of a close relationship between periodontitis and SARS-CoV-2 infection 29 . Periodontalassociated cytokines may proliferate contacts of bacteria between the lungs and the mouth via driving the alteration of the respiratory epithelium and thereby promote respiratory infection in COVID-19 patients 14 . H. parainfluenzae and Kingella, two enriched taxa in the oral cavity of COVID-19 patients, are opportunistic pathogens well-known for respiratory tract infections, infective endocarditis, and meningitis 30,31 . The abrupt loss of Neisseria in patients, which is the highly abundant genus in the normal oral cavity, could raise serious damage to the oral microbiota 32,33 . In addition, all increased bacteria were previously classified as bacteremia-associated bacteria, implying a potential association between oral dysbiosis and secondary bacterial infection in COVID-19 patients 34 . Typically, mechanical ventilation supports may predispose the COVID-19 patients with severe symptoms of dyspnea to pulmonary bacterial co-infection, as the penetration of the device provides an entrance for those opportunistic infectious agents to access the lower respiratory tract from the oral cavity. Given the non-ignorable association between the oral microbiome with bacterial co-infections suggested by our and others' data 35,36 , correct and frequent oral health care measurements should be recommended by physicians to protect COVID-19 patients from secondary infections and improve survival, especially for the ones with the severe disease condition.
Elevated levels of Streptococcus, Rothia, and Actinomyces were identified in COVID-19 patients' feces, which is consistent with previous findings 37 . The host immunity symbionts beneficial bacterial species B. obeum, whose parent genus Blautia was listed as the top decreased taxon in COVID-19 patients' feces in our data, was identified to be depleted in another study 17 . Notably, two normal components of the respiratory tract flora, Granulicatella, and R. mucilaginosa, were identified enriched in patients' feces, which is likely a reflection of the migration of extra-intestinal microbes into the gut or the flourishment of potentially pathogenic bacteria. Typically, the fecal enrichment of specific oral taxa, which was suggested to be linked with the increased oral-fecal microbial transmission, has frequently been regarded as a hallmark of disease 20 . As SARS-CoV-2 viral RNA was detected in stools, the route of viral transmission from the respiratory tract to the intestinal tract can be hinted by the path of the oral flora translocation and colonization 38 . Moreover, R. mucilaginosa seemed further associated with the disease severity in the fecal samples in our study. As some severe patients with COVID-19 linked with cardiovascular disease comorbidities, this finding supports previous work wherein patients with cardiovascular disease comorbidities tended to have a higher prevalence of the Rothia ASV associated with SARS-CoV-2 39,40 . In addition, increased expression of ACE2 receptor was observed in B. longum treated mice; hence changes of B. longum might impact the individual's susceptibility to SARS-CoV-2 through modulating ACE2 expression 41 . It is worth noticing that a decreased relative abundance of B. caccae and B. coprophilus was observed in the gut flora of COVID-19 patients. A recent study by Martino et al. 42 , has demonstrated that Bacteroides may regulate viral adhesion via modifying heparan sulfate (HS), and loss of HS-modifying Bacteroides strains could predispose individuals to SARS-CoV-2 infection. Concordance of these findings with ours provides evidence for a potential relationship between B. caccae and B. coprophilus and susceptibility to SARS-CoV-2 infection. The Mantel and Procrustes analysis indicated the predicted function profiles aligned well with the bacterial taxonomy profile in both oral and fecal samples as expected. Compared to the gut microbiome, our analysis suggested a stronger impairment in the oral microbiome after SARS-CoV-2 infection, wherein a significant reduction of Neisseria, an essential oral microbial genus, was found, along with several fatal metabolic pathways involving the TCA cycle were suppressed. The result echoes a previous study that the gut microbial community possesses higher taxa-function robustness over oral community 43 , since the bacterial community in the gut generally has a higher gene and functional redundancy in comparison to communities inhabiting other body sites.
Microbiome dysbiosis persists during the COVID-19 disease course, even after the viral clearance. Moreover, previous studies have revealed that some recovered patients were re-detectable positive for SARS-CoV-2 RNA after discharge 44 . The cause of redetectable positive remains unclear. One possibility is that a prolonged detrimental effect on the microbiomes exerted by SARS-CoV-2 infection persists even after discharge, which might render convalescent patients susceptible to the residual viremia or reinfection via a more long-lasting dysfunction of the immune system. The microbiome recovery at the taxonomy level seemed to be slowly healed in the oral samples as the PT(−) group shifts towards the healthy controls group based on the dissimilarity distance. In contrast, in the gut samples, the PF(−) group shifts even away from the healthy controls (Fig. 4b). The inconsistency may be due to the period span of sample collection or microbiome characteristics within the oral and gut samples. In a previous metagenomics study, Zuo and colleagues suggested four negatively associated Bacteroides species could be involved in the downregulation of ACE2 expression, and positively regulated E. bacterium 2_2_44A may promote augmenting SARS-CoV-2 infection in the gut 17 . However, a similar analysis method in our study did not identify the same agents to be associated with the viral load of SARS-CoV-2 (Fig. 4c). Thus, we urge validation of the results that should be taken before recognizing the results.
It is well accepted that unnecessary use of antibiotics should be avoided in the episode of acute respiratory infections as antibiotics have no role in treatment. In reality, 67.9% of patients in our study were treated with antibiotics to prevent potential secondary bacterial co-infections at the early phase of the pandemic. Besides, antibiotics prescribing was elevated in patients with severe conditions as compared to non-severe patients. In a meta-analysis study, Langford has reported that 74.6% of COVID-19 patients received antimicrobial treatment 28 . Antibiotics treatment can not only eliminate pathogens but commensal microorganisms indiscriminately, which may lead to microbiota dysbiosis and antimicrobial resistance. To circumvent and minimize the effect of antibiotics on the microbiome analysis, we employed a method to count antibiotics treatment as an adjusting variable in the model or to stratify the data for only patients without antibiotics treatment. In addition, unsupervised dimensionality reduction results demonstrated that COVID-19 patients with and without antibiotics treatment shared similar bacterial community structures with a declined diversity in both oral and gut microbiome, shifting away from healthy individuals. It seems that the impacts derived from infection are overwhelming over that of the antibiotic treatment that may apply to the microbiome. The taxa differential analysis also supported this finding, that the differential taxa between the controls and the antibiotics-free patients were largely overlapped with those identified between the controls and all patients (Supplementary Tables 3 and 5).
Collectively, we reported the alterations in both oral and gut microbiomes of SARS-CoV-2 infected patients during hospitalization and made comprehensive analyses to evaluate their potential consequences and implication. The associations between microbial species with disease severity and viral load in patients have suggested the potential of microbiome-based intervention in the prevention and treatment of COVID-19. We believe that the data provides new knowledge with innovative perspectives for tackling and managing the ongoing COVID-19 pandemic.

Study subject and sample collection
A total of 53 COVID-19 patients and 76 healthy individuals were included in this study. Patients with suspected SARS-CoV-2 infection were confirmed after two sequential positive respiratory tract sample real-time RT-PCR results. Patients were kept hospitalized and under strict observation until the virus was completely eliminated in both respirational and intestinal territories by real-time RT-PCR results. Depending on the sample availability, serial samples were collected from a patient throughout his/ her hospitalization period. More specifically, throat swab and fecal samples were collected in both the positive viral RNA test period (P-VRTP, defined as the period of positive nucleic acid tests until the first day of continuous negative tests) and the negative viral RNA test period (N-VRTP, defined as the interval between the first day of negative nucleic acid test until the hospital discharge) for both sample types. Fecal samples were collected from patients who were ever detected with viral RNA in their feces. Only one sample, either throat swab or fecal specimen, was collected from healthy individuals during their physical examination. None of the COVID-19 patients was received antibiotics nor probiotics within 8 weeks before the infection, and none of the healthy individuals was either before this study recruitment. Patients were categorized into two groups based on disease severity: the non-severe group (mild/moderate) and the severe group (severe/critical) following the instruction of the New Coronavirus Pneumonia Prevention and Control Program (7th edition) published by the National Health Commission of China. The demographic information, underlying diseases, clinical indexes, and treatments were summarized from official patients' medical records. This study was reviewed and approved by the Medical Ethical Committee of the Fifth Affiliated Hospital of Sun Yat-Sen University (approval # K162-1). Written informed consent was obtained from each enrolled subject.

Sample library preparation and sequencing
Samples were inactivated at 56°C for 30 min before DNA extraction. Extraction of nucleic acids was performed with the CFDA approved nucleic acid extraction kits (QIAamp Viral RNA Mini Kit, Catalog #: 52904, QIAGEN). The concentration and the purity were measured using the NanoDrop One (ThermoFisher Scientific, MA, USA). Sequencing libraries were generated using NEBNext® UltraTM II DNA Library Prep Kit for Illumina (New England Biolabs, MA, USA) following manufacturer's recommendations, and index codes were added. The universal primers 338F (5′-ACTCCTACGGGAGG-CAGCA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) were used for amplification of the V3-V4 region of the bacterial 16S rRNA gene. DNA libraries were generated from PCR amplicons targeting the hypervariable regions V3-V4 of the bacterial 16S rRNA gene. After quality assessment with the Qubit 2.0 Fluorometer (ThermoFisher Scientific, MA, USA), the library was then sequenced on an Illumina NovaSeq 6000 platform, and a minimum of 50,000 of 250 bp paired-end reads was generated for the samples.

Bioinformatics analysis
Raw FASTQ files were demultiplexed using the QIIME 2 demux plugin based on their unique barcodes 45 . Demultiplexed sequences from each sample were stitched, quality filtered, trimmed, de-noised, and then the chimeric sequences were identified and removed using the QIIME 2 dada2 Y. Wu et al. plugin to obtain the feature table 46 . The QIIME 2 feature-classifier plugin was then used to align feature sequences to a pre-trained GREENGENES 13_8 99% database to generate the taxonomy table 47 . The data was rarefied prior to alpha and beta diversity analysis using a depth of 13,111 reads. Any contaminating mitochondrial and chloroplast sequences were filtered using the QIIME 2 feature-table plugin. Diversity metrics were calculated and plotted using the core-diversity plugin and the emperor plugin within QIIME 2 48 . The beta diversity significance among groups was examined with PERMANOVA and PERMDISP tests using QIIME 2 plugins and Vegan package in R (version 4.0.2). The differences in the relative abundance of taxa between the patients and healthy control groups were identified in both sample types separately using the linear discriminant analysis effect size (LEfSe) 49 . Given the possible confounding impact, the output from LEfSe was further validated with Multivariate Association with Linear Models (MaAsLin2) by adjusting confounding factors (age, sex, antibiotic usage, PCR detection result, and patient ID), so only the bacteria taxa agreed by both methods were presented 50 . PICRUSt2 was used to predict the microbial metabolic pathways to assess the potential functional implication 51 . QIIME 2 Procrustes plugin was used to examine the fitness of the functional properties and the bacterial composition in the microbiome community with the Procrustes plot, and the correlation between functional properties (Bray Curtis distance) and the bacterial composition (unweighted UniFrac distance) in the microbiome community was examined by two-sided Mantel test.

Detection of SARS-CoV-2 viral load
SARS-CoV-2 viral loads in the throat and fecal swabs were measured using a real-time RT-PCR assay. Viral RNA from throat swabs and fecal samples were extracted using QIAamp Viral RNA Mini Kit (QIAamp Viral RNA Mini Kit, Catalog #: 52904, QIAGEN). Up to 0.1 g of stool or throat swab was suspended in a 2 mL viral transport medium (in 1:10 dilution), followed by centrifugation at 3000 × g for 30 min. The aliquot of the filtrate was used as the starting material. The real-time RT-PCR was carried out with the Novel Coronavirus (2019-nCoV) real-time RT-PCR kit from LifeRiver Ltd. (Catalog #: RR-0479-02). Nucleocapsid gene (N), membrane gene (E), and RNA dependent RNA polymerase gene (RdRp) were the three targeted genes simultaneously amplified and tested. According to the manufacturer's instructions, a combined result of the three SARS-CoV-2 viral gene targets was used to yield a positive result. Samples were considered negative if the cycle threshold values exceeded 43 cycles. Plasmids containing the full N gene and E gene were obtained to assess SARS-CoV-2 viral copy in the samples (PCDNA6B-SARS-CoV-2-N and PCDNA6B-SARS-CoV-2-E, gifts from Peihui Wang, Cheeloo College of Medicine, Shandong University). Serial 10fold dilutions of known copies of these plasmids were prepared separately for generating the standard curve. The Ct values of real-time RT-PCR from patient samples were converted into viral RNA copies based on a standard curve.

Statistical analysis
Categorical variables were presented as numbers and percentages (n/N, %) whereas continuous variables were reported as median and interquartile ranges (IQR). According to the distribution of data sets, the one-way analysis of variance (ANOVA) followed by Tukey's multiple comparisons test or the Kruskal-Wallis test followed by Dunn's multiple comparisons test were performed using Prism 8.3.0 (GraphPad Software). Given the nature of compositional data for the 16S microbial relative abundance information, we implemented CCLasso to define the correlation between the longitudinal viral load and the timepoint-matched relative abundances of different taxa in the throat swab and fecal samples separately 52 . The Spearman rank test was used to calculate the correlation between the viral load represented by the abundances of the N gene and the E gene, as their values did not follow the normal distribution according to the Kolmogorov-Smirnov test. All statistical tests were two-sided, and a p-value, of <0.05 was considered significant. Statistical analysis was performed using SPSS version 25.0 (SPSS Inc). Additional custom R scripts were used to making the plots.

Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.

DATA AVAILABILITY
Raw FASTQ files of 16S rRNA gene sequences were archived in the Sequence Read Archive (SRA) under Bioproject accession number PRJNA684070. All supporting data are publicly available at https://github.com/XiaominCheng/COVID19-16S.