SARS-CoV-2 infection reduces human nasopharyngeal commensal microbiome with inclusion of pathobionts

The microbiota of the nasopharyngeal tract (NT) play a role in host immunity against respiratory infectious diseases. However, scant information is available on interactions of SARS-CoV-2 with the nasopharyngeal microbiome. This study characterizes the effects of SARS-CoV-2 infection on human nasopharyngeal microbiomes and their relevant metabolic functions. Twenty-two (n = 22) nasopharyngeal swab samples (including COVID-19 patients = 8, recovered humans = 7, and healthy people = 7) were collected, and underwent to RNAseq-based metagenomic investigation. Our RNAseq data mapped to 2281 bacterial species (including 1477, 919 and 676 in healthy, COVID-19 and recovered metagenomes, respectively) indicating a distinct microbiome dysbiosis. The COVID-19 and recovered samples included 67% and 77% opportunistic bacterial species, respectively compared to healthy controls. Notably, 79% commensal bacterial species found in healthy controls were not detected in COVID-19 and recovered people. Similar dysbiosis was also found in viral and archaeal fraction of the nasopharyngeal microbiomes. We also detected several altered metabolic pathways and functional genes in the progression and pathophysiology of COVID-19. The nasopharyngeal microbiome dysbiosis and their genomic features determined by our RNAseq analyses shed light on early interactions of SARS-CoV-2 with the nasopharyngeal resident microbiota that might be helpful for developing microbiome-based diagnostics and therapeutics for this novel pandemic disease.


Results
The first step in managing COVID-19 is the rapid and accurate detection of SARS-CoV-2 by quantitative reverse transcription-polymerase chain reaction (RT-qPCR). In this study, patients were diagnosed positive for COVID-19 on an average 4.7 days (range 2-9 days) after the onset of clinical signs, and these patients tested negative for SARS-CoV-2 (Recovered) on an average 17.5 days (range 11-32 days) after the initial COVID-19 confirmatory www.nature.com/scientificreports/ diagnosis (Table S1). The Healthy control subjects however did not show any signs and symptoms of respiratory illness.

SARS-CoV-2 infection modulates the community composition and diversity of nasopharyngeal microbiomes.
To shed light on the effects of SARS-CoV-2 infections in the diversity and composition of human nasal microbiome, we analyzed RNAseq data of nasopharyngeal samples of randomly selected healthy human, COVID-19 and recovered patients (Fig. 1). The alpha-diversity (within sample diversity through Shannon and Simpson indices) assessment revealed significant differences in microbial species richness across three metagenomes regardless of the method used to tabulate microbial abundances i.e., either IDseq or MG-RAST (MR), showing higher diversity in the microbial niche of Recovered (p = 0.0065; Kruskal-Wallis test) followed by Healthy (p = 0.0041; Kruskal-Wallis test) and COVID-19 (p = 0.0105; Kruskal-Wallis test) samples ( Fig. 2A). The beta diversity (between sample or metagenome diversity) through principal coordinate analysis (PCoA), as measured on the Bray-Curtis distance method using IDseq and MR data, showed distinct discrimination across the metagenomes and separated samples by microbial population structure (Fig. 2B). Therefore, significant variation in microbiome diversity and composition across these three metagenomes (p = 0.0059, Kruskal-Wallis test) was evident. The diversity of microbiomes between samples (PCoA) however did not vary significantly (p = 0.651, Kruskal-Wallis rank sum test) according to the gender (male or female) of the study population (Fig. 2B). Microbiome composition at the domain level was numerically dominated by bacteria, with a relative abundance of 81.58%, followed by viruses (18.22%) and archaea (0.20%) (Fig. 1, Data S1). The unique and shared distribution of microbes found in the three participant groups are represented by comprehensive Venn diagrams (Fig. 3). In this study, we detected 532 bacterial genera, including 486, 421 and 509 in COVID-19, Recovered and Healthy nasopharyngeal samples, respectively, and of them, 72.74% genera were common in three sample groups (Fig. 3A, Data S1). Notably, we detected 2281 bacterial species through IDseq analysis, of which 919, 675 and 1476 species were found in COVID-19, Recovered, and Healthy samples, respectively (Fig. 3B, Data S1). However, compared to COVID-19 and Recovered samples, the Healthy samples had unique or sole association of 1010 (68.43%) bacterial species which underwent dysbiosis during the pathogenesis and subsequent recovery phase of COVID-19 (Data S2). Of the identified bacterial species, 50.16% and 29.04% had sole association in COVID- 19 and Recovered samples, respectively indicating their opportunistic inclusion in COVID-19 patients and re-establishment beneficial commensal flora with the recovery of SARS-CoV-2 infections (Fig. 3B, Table 1, Data S2).

SARS-CoV-2 infection associated changes in genomic potentials of the microbiomes. In our
current RNAseq data set, there was a broad variation in resistance to antibiotics and toxic compounds (RATCs) diversity and composition across three metagenomes (Fig. 8, Data S3). The categories and relative abundances of the RATCs were significantly correlated (p = 0.027, Kruskal-Wallis test) with the relative abundance of the associated microbiotas found in COVID-19, Recovered and Healthy control samples (Data S3). We detected 49 RATCs through MR analysis distributed in the microbial genomes of the three metagenomes (Fig. 8, Data S3). The COVID-19 associated microbiomes harbored the 40 different RATCs while 32 and 44 gene families were detected in the microbiomes of Recovered and Healthy controls. Of the detected RATC groups, genes associated  (Fig. 8, Data S3). By examining the correlation between the different gene families of the same KEGG pathway for COVID-19, Recovered and Healthy controls nasopharyngeal microbiomes, we found significant differences (p = 0.034, Kruskal-Wallis test) in their composition and relative abundances. Our analysis revealed that genes coding for pyruvate carboxylase (pyc) had several-fold overexpression in the microbiomes of the COVID-19 (26.31%) and Recovered (4.11%) groups compared to Healthy controls (1.80%) (Fig. 9A, Data S3). In addition, genes encoding for adherent junction (25.17%), tight junction (24.48%), environmental information processing (15.03%), carbohydrate metabolism (11.51%) and oxidative phosphorylation (7.26%) had higher relative abundances in the COVID-19 patient's nasopharyngeal microbiomes (Fig. 9A, Data S3). The nasopharyngeal microbiomes of We also sought to gain further insight into the SEED hierarchical protein functions represented by different genes, and found 37 statistically different (p = 0.013, Kruskal-Wallis test) SEED functions in COVID (COVID-19 and Recovered) and Healthy control metagenomes. Overall, the COVID-19 and Recovered samples associated microbiomes showed a higher relative abundance of these SEED functions compared to those of Healthy controls. For instance, cytokine-cytokine receptor interaction (50.0%), regulation of oxidative stress response (26.82%), phage integration and excision (24.92%), toxin-antitoxin regulation system (17.24%), protection from reactive oxygen species (14.96%), and phage regulation (6.06%) related genes had more than two-time overexpression in Healthy controls commensal nasopharyngeal microbiomes compared to COVID-19 and Recovered humans nasopharyngeal microbiomes (Fig. 9B, Data S3). The COVID-19 patients nasopharyngeal microbiomes were enriched in genes coding for cell growth and death (55.60%), Streptococcus pyogenes virulence regulators (42.44%), fratricide in Streptococcus (32.67%), neuroactive ligand receptor interaction (28.21%), epithelial cell signaling (20.35%), clustering-based subsystems (14.65%), proteolytic pathways (9.2%), prophage lysogenic conversion modules (3.69%) and bacterial invasion of epithelial cells (1.93%) compared to rest of the two metagenomes (Fig. 9B, Data S3). The Recovered humans nasopharyngeal microbiomes however had a higher abundance of SEED functions involved in glutathione: non-redox reactions (30.30%) and redox cycle (6.06%), coagulase cascade (20.0%), osmotic stress (16.07%), membrane transport (14.55%), MT1-MMP pericellular network (10.0%), BarA-UvrY(SirA) two-component regulatory system (9.29%), and oxidative stress (5.54%). In contrast, most of these SEED modules showed comparatively lower relative abundance in Healthy human nasopharyngeal microbiomes (Fig. 9B, Data S3). www.nature.com/scientificreports/ Discussion SARS-CoV-2 infection may predispose to secondary microbial infection which is associated with poor clinical outcome especially among critically ill patients. Therefore, the complex crosstalk between SARS-CoV-2 and commensal microbes of different body systems is essential for the functioning of the immune system. In the present study, we demonstrated a remarkable shift in the diversity and composition of the nasopharyngeal microbiomes (bacteria, viruses and archaea) of COVID-19 and Recovered people compared to the Healthy humans through the state-of-the-art RNAseq technology. The findings of the present study revealed that SARS-CoV-2 infection reduces commensal bacteria with inclusion of pathobionts in the nasopharyngeal tract of human (Fig. 1). Furthermore, we identified a number of microbial genomic features, altered metabolic pathways, and functional genes associated with COVID-19 pathogenesis. Although the dysbiosis of human gut microbiome by the infection of SARS-CoV-2 has been reported in several earlier studies 8,21,22 , this study for the first time determined the interactions and consequences of SARS-CoV-2 infection with resident nasopharyngeal microbiome of human. Since some of the COVID-19 patients showed persistent symptoms after recovery and/or subsequently develop multisystem inflammation 8 , we hypothesized that the altered nasopharyngeal microbiomes during SARS-CoV-2 infections could lead to abnormal inflammatory reactions to worsen the symptoms and treatment outcomes of COVID-19. The microbiome diversity (alpha and beta diversity) measures suggest that microbial dysbiosis is closely linked to SARS-CoV-2 infections. Our analysis also showed a substantial microbial disparity between COVID-19 and Healthy-controls nasal samples keeping the closest relationship between COVID-19 and Recovered people's nasopharyngeal samples, and related microbial signatures. The COVID-19-patients nasopharyngeal microbiomes had significantly lower variance in diversity than those of healthy humans, agreeing with several recent studies 23,24 . In a recent study, Zou et al. reported that loss of salutary species in COVID-19 pathogenesis in most patients, despite the clearance of SARS-CoV-2 virus, is associated with a more long-lasting detrimental effect to the gut microbiome 24 . The nasopharyngeal microbiomes predominantly identified in this study are bacteria (> 81. 0%), however, other concomitant microbial players like viruses and archaea were also detected, highlighting the novel insights of diverse microbial association with SARS-CoV-2 to exacerbating the course, common symptoms and treatment outcomes of the COVID-19 22,25,26 . In COVID-19 samples, the relative abundances of Firmicutes, Bacteroidetes and Proteobacteria phyla, and their related genera such as Streptococcus, Figure 8. Distribution of the resistance to antibiotic and toxic compounds (RATC) genes in COVID-19, recovered and healthy control metagenomes. The circular plot illustrates the diversity and relative abundance of the RATC genes detected among the microbiomes of the three metagenomes through SEED subsystems analysis. The association of the RATC genes according to metagenome is shown by different colored ribbons and the relative abundances these genes are represented by inner blue colored bars. Some of the RATC functional groups shared among microbes of the three metagenomes (COVID-19, Recovered and Healthy), and rest are effectively undetected in the microbiomes of the other metagenomes. The numerical values 1-49 represent different RATC functional groups found across the study metagenomes. More data on RATC functional groups and their relative abundances are also available in Data S2. The circular plot is built using the OmicCircos (https:// bioco nduct or. org/ packa ges/ relea se/ bioc/ html/ OmicC ircos. html), an R package.  (Table 1). For instance, Streptococcus and Veillonella had more than five-fold higher relative abundances in COVID-19 samples than the Recovered and Healthy samples. Prevotella, Veillonella, and Streptococcus are the most common genera that reside in the lungs and nasopharynx 27,28 , and recently these genera have been reported to play an opportunistic role in the progression of lung infections 28,29 . Despite having almost homogeneous genetic backgrounds and living in the same region, there were vast differences in microbiome signatures in nasal cavities of Healthy controls, COVID-19 patients and Recovered individuals. The Healthy people had higher number of bacterial genera compared to COVID-19 patients and Recovered individuals. Moreover, the inter-individual microbiome signature (bacterial genera) of participants was also observed, and was largely stable in COVID-19 patients. Notably, COVID-19 patients mainly had increased numbers of opportunistic pathogens, a part of commensal microbiota that may become pathogenic in the event of host perturbation, such as dysbiosis or immunocompromised host. The human microbiota shows a remarkable amount of diversity among different individuals, and there are ample of emerging evidences that the microbial communities can be altered to change pathophysiological state or even cure of disease is a compelling drivers of microbiome variation 30 . Moreover, co-infections in COVID-19 patients could also be attributed by different variants and strains of SARS-COV-2 5,31 . The attributes of microbiomes, for instance, pathogenicity, virulence, antibiotic resistance, and metabolic potentials are linked with the species or strain specific genomic characteristics 22,23,32 . We demonstrated that the Healthy samples possessed highest number of commensal bacterial species while SARS-CoV-2 infections and showing the average relative abundance hierarchical clustering of the predicted SEED functions in different levels among the microbiomes of four metagenomes. The color bars (column Z score) at the top represent the relative abundance of putative genes. The color codes indicate the presence and completeness of each KEGG and SEED module, expressed as a value between − 3 (lowest abundance) and 3 (highest abundance). The red color indicates the more abundant patterns, whilst green cells account for less abundant putative genes in that particular metagenome. The Heatmap is built through a stand-alone software tool called FunRich (http:// www. funri ch. org/). www.nature.com/scientificreports/ subsequent therapy reduced 37.78% and 54.23% bacterial species in COVID-19 patients and Recovered humans nasopharyngeal samples, respectively indicating the perturbations of beneficial microbes amid COVID-19. The COVID-19 patients and Recovered people had sole association of 50.16% and 29.04% bacterial species, respectively. These results are in line with several previous studies who reported that cross-talk and/or interactions between SARS-CoV-2 and oral microbiomes 26 , and SARS-CoV-2 and gut microbiomes 33 is associated with the pathophysiology of lung diseases. Therefore, majority of the bacterial species identified in COVID-19 (66.26%) and Recovered (76.78%) samples had an opportunistic inclusion by the depletion of beneficial commensal microbes during the course of SARS-CoV-2 pathogenesis and recovery phase. The composition of the bacterial communities in the nasopharynx is more diverse than any other parts of the upper respiratory tract, which are characterized by different Streptococcal species, Staphylococcus spp., Haemophilus spp., Neisseria spp., Rothia spp., and anaerobes, including Veillonella spp., Prevotella spp., and Leptotrichia spp. 18,34,35 . Despite some preliminary evidence, it is still unclear whether similar synergies exist between SARS-CoV-2 and specific bacterial infections as exist with influenza 36,37 . Surprisingly, 68.43% bacterial species had sole association with Healthy nasopharyngeal samples, which were not detected in COVID-19 patients and Recovered humans indicating the potential dysbiosis of commensal microbiomes during pathophysiology of SARS-CoV-2 infections. In this study, different species of Pedobacter, Sphingobacterium, Pseudomonas, Enterobacter, and Flavobacterium genera remained top abundant in Healthy human nasopharyngeal swabs. These microbes are ubiquitous and isolated from natural environments, such as soil and water, readily cleared by host defenses and rarely involved in human infection 10,18,34,38 . Though, the nature of nasopharyngeal commensal bacteria that exert antiviral effect in the lungs is still elusive, some members of these commensal species (desaminotyrosine producing anaerobes) are important in the priming of the pulmonary innate immune system 22 . Despite the beneficial role of these commensal bacterial species are largely unknown, reduction in the composition and relative abundances of these commensal bacteria from the nasopharyngeal tract, play an important role in microbiome formation, which can significantly affect the immune response against SARS-CoV-2 infections 39 .
The cutting-edge RNAseq approach of this study also provided an exciting opportunity for investigating the integrated cross-kingdom interactions of 'multibiome' such as viruses (other than SARS-CoV-2) and archaea which contributed about 20.0% of the total microbial sequence reads. Unlike bacteria, the diversity and composition of viruses and archaea always remained much lower in both COVID-19 and Healthy samples. In this study, of the detected viral reads (18.0%), the Alphacoronvirus (HCoV-NL63) and Betacoronovirus (SARS-CoV-2 and Pangolin coronavirus) comprises majority (> 85%) sequences indicating that both Alpha-and Beta-coronaviruses are simultaneously prevailing in Bangladesh. SARS-CoV-2 was the predominantly abundant (> 99.0% relative abundances) viral strain in COVID-19 and Healthy controls (84%) metagenomes, and rest of the viral strains concurrently detected with SARS-CoV-2 were belonged to Alphacoronavirus and Siphoviruses. On the contrary, HCoV-NL63 was the most abundant viral strain (90.89%) in Recovered humans' nasopharyngeal samples followed by SARS-CoV-2 (2.76%). Human coronaviruses (HCoVs) have only recently been shown to cause both lower and upper respiratory tract infections. The endemic HCoV-NL63 can cause co-infections, sequential infections or can be co-detected with each other or with other respiratory viruses like SARS-CoV-2 40,41 . Recent studies reported that infection with HCoV-NL63 is common worldwide, and associated with many clinical symptoms in the immunocompromised SARS-CoV-2 patients 41,42 . Other viral strains identified in the samples of COVID-19, Recovered and Healthy metagenomes were mostly represented by different strains of bacterial phages. Recently, stronger relationship among SARS-CoV-2, bacteria, fungi, and other viruses have been reported from different countries [43][44][45] . The respiratory tract microbiomes support large populations of viruses, mostly bacteriophages (phages) that infect specific bacterial species, and can appear as free phage virions (phage particles) or as dormant prophages. Since most secondary bacterial infections occur in immunocompromised or immunodeficient SARS-CoV-2 patients, phages can induce immune responses these patients 46 . The archaeal fraction of the human nasopharyngeal microbiomes across three metagenomes were predominated by different methanogenic and thermophilic genera. Though the role of these accompanying microbiomes in the pathophysiology of COVID-19 remains a mystery, it may include virus-induced airway damage, cell loss, goblet cell hyperplasia, altered mucus secretion, reduced ciliary beat frequency, function and clearance, reduced oxygen exchange, and damage to the immune system 47 .
With advancement in metagenomic sequencing techniques and downstream bioinformatics analyses, it is now possible to critically examine and analyze various microbial communities and their concurrent antimicrobial resistance in different body sites including the NT. One of the important findings of the present study is the concurrent detection of various homologs of RATCs belonging to different protein families among the microbiome of three metagenomes. Interestingly, the composition and relative abundance of the RATCs detected in this study remained significantly correlated with microbiome signature and their relative abundances in the nasal cavities of the Healthy humans, COVID-19 patients and Recovered individuals. The observed variations in the composition and relative expression of RATCs in COVID-19, Recovered and Healthy controls corroborates with the dynamic dysbiosis of microbiomes in the corresponding metagenomes, their genetic diversity, and selective pressure for the maintenance of antimicrobial resistance. The microbiomes of COVID-19 patients and Recovered humans nasopharyngeal swabs were enriched with genes coding for biofilm formation in Staphylococcus, quorum sensing: autoinducer-2 synthesis, mdtABCD, cobalt-zinc-cadmium, acriflavine, arsenic, fluoroquinolones and methicillin resistance, BlaR1 regulatory family, macA, MAR locus, pmrA, YjgK cluster, and lsrACDBFGE operon. Conversely, beta-lactamase resistance, quorum sensing in Yersinia and Pseudomonas, and teicoplanin-resistance in Staphylococcus encoding genes remained predominantly abundant in nasopharyngeal microbiomes of the Healthy controls. Remarkably, since the beginning of the COVID-19 pandemic, there has been growing concern for a potential rise in AMR secondary to increased antibiotic prescription for COVID-19 patients 48 . Moreover, severe COVID-19, which particularly affects elderly patients with multiple comorbidities, may be an important factor in determining changes in colonization pressure and multidrug resistance 49  www.nature.com/scientificreports/ The functional analysis of the microbiomes revealed that genes coding for pyruvate carboxylase (pyc), adherent junction, tight junction, environmental information processing, carbohydrates metabolism, and oxidative phosphorylation had higher relative abundances in the COVID-19 patients nasopharyngeal microbiomes compared to those of Recovered and Healthy control microbiomes. These metabolic functional changes further lead to increased cytokine-cytokine receptor interaction, regulation of oxidative stress response, phage integration and excision, toxin-antitoxin regulation, protection from reactive oxygen species, and phage regulated gene expression in COVID-19 associated microbiomes as also reported previously in other viral diseases 50,51 . The Recovered humans NT microbiomes were enriched with genes coding for focal adhesion, transport and catabolism, cell adhesion molecules, genetic information and processing, lysosome activity, endocytosis, cell-tocell communication and glycolysis. The pertinent genomic potentials of the COVID-19 patient and Recovered humans nasopharyngeal microbiomes were also evidenced by the higher expression of genes involved in glutathione: non-redox reactions, redox cycle, coagulase cascade, osmotic stress, membrane transport, MT1-MMP pericellular network, and BarA-UvrY (SirA) two-component regulatory activities. Conversely, the Healthy people nasopharyngeal microbiomes showed a lower abundance and expression of these metabolic functional genes. The metabolic health of an individual is represented by the proper functioning of organismal metabolic processes coordinated by multiple physiological systems 51 . The differentially abundant functions and pathways identified in this study corroborates with the findings from previous reports, relating to decreased lipid and glycan metabolism, increased carbohydrate metabolism, and other characteristics of the microbiome linked to COVID-19 52 . Several predicted functional pathways differed between COVID-19 and Healthy controls, perhaps reflecting metabolic changes associated with the progression of COVID-19 pathogenesis, and novel host-microbiome interactions in SARS-CoV-2 infected patients.

Conclusions
Host-microbiomes interactions are exceptionally complex especially in SARS-CoV-2 infections. Our RNAsSeq metagenomic analyses revealed that SARS-CoV-2 infection had significant effect on the diversity and composition of nasopharyngeal microbiomes of human. The identifiable changes in the microbiome diversity, composition and associated genomic features demonstrated in this study might be associated with the development, treatment, and resolution of COVID-19. The SARS-CoV-2 infection results in remarkable depletion of nasopharyngeal commensal microbiomes of Healthy humans with inclusion of different opportunistic pathogens in COVID-19 and Recovered samples. Several predicted functional pathways differed between COVID-19 patient and Recovered people nasopharyngeal samples compared to healthy individuals which reflect the roles microbial metabolic changes with the progression of SARS-CoV-2 pathogenesis. These interactions are further complicated by the common co-existence of microbiota that interact with both host and SARS-CoV-2. Therefore, findings of this study may serve as a benchmark for microbiome-based diagnostic markers and therapeutics for this pandemic disease. However, future time-course studies with a larger sample size are needed to elucidate the dynamic changes in composition and diversity of commensal microbiomes, and the inclusion of pathobionts in the whole respiratory system and gut during the progression of COVID-19.

Methods
Subject recruitment and sample collection. This study was conducted in accordance with the guidelines of the Director General of Health Services (DGHS) of Bangladesh during May to July, 2020. The study participants provided written informed consent consistent with the experiment. Twenty-two (n = 22) nasopharyngeal samples (including COVID-19 = 8, Recovered = 7, and Healthy = 7) were collected from Dhaka city of Bangladesh. The patients were diagnosed positive for COVID-19 on an average 4.7 days (range 2-9 days) after the onset of clinical signs (SARS-CoV-2 positive through RT-qPCR). On an average, 17.5 days (range 11-32 days) after the first confirmation of SARS-CoV-2 infection these patients were tested negative for COVID-19 and categorized as recovered people (Table S1). The average age of the study people was 41.86 (range 22-72) years, and of them, 68.18% and 38.82% were male and females, respectively (Data S1). The RT-qPCR was performed for ORF1ab and N genes of SARS-CoV-2 using novel Coronavirus (2019-nCoV) Nucleic Acid Diagnostic Kit (PCR-Fluorescence Probing, Sansure Biotech Inc.) according to the manufacturer's instructions. Viral RNA was extracted using a PureLink viral RNA/DNA minikit (Thermo Fisher Scientific, USA). Thermal cycling was performed at 50 °C for 30 min for reverse transcription, followed by 95 °C for 1 min, and then 45 cycles of 95 °C for 15 s, 60 °C for 30 s on an Analytik-Jena qTOWER instrument (Analytik Jena, Germany).
Nasopharyngeal samples from COVID-19 and Recovered subjects were collected on the test day (day of COVID-19 positive and COVID-19 negative identification), and subsequently sent for RNA extraction and sequencing. The Healthy control subjects were randomly selected and these people did not show any signs and symptoms of respiratory illness, and nasopharyngeal swabs from these Healthy people were collected following the same protocol for COVID-19 and Recovered humans, and immediately sent for RNA extraction and sequencing.
RNA sequencing. We utilized total RNAseq approach for the metagenomics component of the study. The www.nature.com/scientificreports/ Taxonomic mapping, classification, diversity and community analysis. The paired-end sequences of COVID-19, recovered and healthy samples (n = 22) were analyzed using both mapping-based and assembly-based hybrid methods of IDseq 53 and MG-RAST (release version 4.1) (MR) 54 . IDseq-an open-source cloud-based pipeline has been used to assign taxonomy with NT L (nucleotide alignment length in bp) >= 50 and NT %id >= 97. This pipeline used quality control, host filtering, assembly based alignment and taxonomic reporting aligning to NCBI nt (nucleotide database). The sequencing reads were filtered through BBDuk (with options k = 21, mink = 6, ktrim = r, ftm = 5, qtrim = rl, trimq = 20, minlen = 30, overwrite = true) to remove Illumina adapters, known Illumina artifacts, and phiX. Any sequence below these thresholds or reads containing more than one 'N' were discarded 32 . In the current study, 120.50 million reads were generated from the samples of three metagenomes with an average of 5.0 million reads per sample passed the quality control steps, and of these quality reads, an average of 3.22 million aligned reads per sample mapped to the ribosomal (rRNA) reference gene libraries (Data S1). Alpha diversity (diversity within samples) was estimated using the Shannon and Simpson diversity indices. To visualize differences in bacterial diversity, a PCoA (at genus level) was performed based on the Bray-Curtis distance method. Microbial taxa that were detected in one group of sample but not detected in rest of the two groups are denoted as solely (unique) associated microbiomes 32 .
Functional profiling of the microbiomes. We performed the taxonomic functional classification through mapping the reads onto the Kyoto Encyclopedia of Genes and Genomes (KEGG) database 55 , and SEED subsystem identifiers 54 on the MR server using the partially modified set parameters (e-value cutoff: 1 × 10 -30 , min. % identity cutoff: 60%, and min. alignment length cutoff: 20) 56 .

Statistical analysis.
The non-parametric test Kruskal-Wallis rank sum test was used to evaluate differences in the relative percent abundance of taxa in COVID-19, Recovered and Healthy sequence data. The gene counts were normalized by dividing the number of gene hits to individual taxa/function by total number of gene hits in each metagenome to remove bias due to differences in sequencing efforts. In addition, statistical tests were also applied with non-parametric test Kruskal-Wallis rank sum test at different KEGG and SEED subsystems levels through IBM SPSS (SPSS, Version 23.0, IBM Corp., NY, USA) 32 .
Ethical statement. The protocol for sample collection from COVID-19, recovered and healthy humans, sample processing, transport, and RNA extraction was approved by the National Institute of Laboratory Medicine and Referral Center of Bangladesh.

Data availability
All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. The sequence data reported in this article has been deposited in the National Center for Biotechnology Information (NCBI) under BioProject accession number PRJNA720904.