Respiratory Microbiome Dynamics in COVID-19: A Comprehensive Multi-Omics Study

Respiratory Microbiome Dynamics in COVID-19: A Comprehensive Multi-Omics Study | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Article Respiratory Microbiome Dynamics in COVID-19: A Comprehensive Multi-Omics Study Deyvid Amgarten, Fernanda de Mello Malta, Raquel Riyuzo de Almeida Franco, and 19 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-8013295/v1 This work is licensed under a CC BY 4.0 License Status: Under Review Version 1 posted You are reading this latest preprint version Abstract The outcomes of SARS-CoV-2 infection exhibit significant heterogeneity, and the role of the airway microbiota remains insufficiently understood. In a cohort comprising 561 participants, an integrated multi-omics approach was employed. This approach combined 16S/ITS amplicons (n = 542), total RNA metatranscriptomics (n = 495), and total DNA metagenomics (n = 113), with 101 samples analyzed using all three omics. Metatranscriptomes resolved lineage dynamics across major infection waves and identified microbial coinfections; shotgun metagenomes recovered 190 MAGs, including 13 putative novel species, expanding the catalog of respiratory taxa. Cross‑omics testing showed species‑level biomarkers and resistome burden were more linked to clinical outcomes than genus‑level signals, highlighting the value of high‑resolution profiling. Functional analyses revealed outcome‑based shifts in microbial processes and host metatranscriptomic signatures enriched for interferon programs and tissue remodeling pathways. This work delineates mechanisms associated with COVID‑19 severity, uncovers candidate biomarkers, and provides a multi‑omics framework for integrated pathogen–microbiome–host surveillance in future respiratory pandemics. Biological sciences/Computational biology and bioinformatics/High-throughput screening Health sciences/Diseases/Infectious diseases/Viral infection Health sciences/Biomarkers/Prognostic markers Biological sciences/Microbiology/Microbial communities/Clinical microbiology Biological sciences/Microbiology/Microbial communities/Metagenomics Respiratory Microbiome COVID-19 SARS-CoV-2 Multi-Omics 16S/ITS Amplicon Sequencing RNA Sequencing Metagenomics Disease Severity Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 1. INTRODUCTION The COVID-19 pandemic, driven by SARS-CoV-2, has created an unprecedented global health crisis, characterized by clinical outcomes ranging from asymptomatic infection to severe disease. While significant research has focused on viral evolution and host immune responses, the influence of the respiratory microbiome on COVID-19 severity remains understudied 1 – 3 . The human respiratory tract harbors a diverse microbiome that is fundamental for maintaining health. In healthy individuals, the upper respiratory tract is predominantly colonized by oral commensal microbes, whereas the lower respiratory tract generally maintains a lower biomass 4 . The respiratory microbiome is continuously shaped by microbial migration from the upper tract and is balanced by clearance mechanisms involving mucociliary activity and innate immune responses. Disruption of this delicate equilibrium can lead to pathogenic overgrowth, contributing to respiratory illnesses, such as pneumonia and chronic lung diseases 4 , 5 . The respiratory microbiome plays a pivotal role in modulating local and systemic immune responses. Studies have demonstrated that upper respiratory tract microbiota can significantly influence immune regulation in the lungs, a relationship that is particularly relevant in respiratory infections such as COVID-19 6,7 . Microbial dysbiosis, characterized by an imbalance in microbial community composition, correlates with increased inflammation and is implicated in exacerbating conditions such as asthma and chronic obstructive pulmonary disease (COPD). Factors such as infections, environmental exposures, and therapeutic interventions can perturb the microbiome, disrupt respiratory homeostasis, and potentially influence disease progression 8 – 10 . In COVID-19, the relationship between respiratory microbiome composition and disease severity has become an important research area. Emerging evidence suggests that specific microbiome profiles may modulate the immune responses critical for controlling SARS-CoV-2 infection 6 , 11 . Certain oropharyngeal taxa have been linked to immune activation patterns associated with varying disease severities. Additionally, SARS-CoV-2 infection may alter the respiratory microbiome, leading to dysbiosis that contributes to complications, such as bacterial pneumonia. Conversely, specific microbial communities may confer protective effects, mitigating severe COVID-19 outcomes 12 , 13 . Clarifying the role of the respiratory microbiome could provide insights into COVID-19 pathogenesis and potentially guide therapeutic strategies. To address this gap, we conducted a large-scale, multi-omics study of nasopharyngeal samples collected from 561 individuals, including 544 COVID-19 patients and 17 uninfected controls. By integrating 16S/ITS amplicon sequencing, shotgun metatranscriptomics sequencing, and shotgun metagenomics, our approach enabled us to track SARS-CoV-2 lineage dynamics, uncover co-infections, and delineate shifts in microbial diversity associated with the disease severity. In doing so, we not only recovered high-quality metagenome-assembled genomes (MAGs) and characterized functional gene profiles, but we also identified potential microbial biomarkers linked to disease progression. Our findings reveal key host–microbe interactions and metabolic alterations in the respiratory tract during COVID-19, suggesting that the microbiome may play a role in modulating clinical outcomes of the disease. This study provides a comprehensive resource that lays the groundwork for future therapeutic and prognostic strategies targeting the respiratory microbiome in viral infections. 2. RESULTS 2.1. Cohort and experimental design Swab samples were collected from 561 participants between July 2020 and July 2022. Seventeen participants formed the control group (reverse transcription – polymerase chain reaction (RT-PCR) negative, no respiratory symptoms), while 544 had confirmed SARS-CoV-2 infection. Age ranged from 1–98 years (mean of 39.3), increasing significantly with disease severity (p < 0.05): asymptomatic (33.9) vs. severe (57.5). Body Mass Index (BMI) ranged 17–44, with higher means in severe (29.1) and moderate (29.3) groups than in mild (26.2) or asymptomatic (28.4) (p < 0.05). Gender distribution differed significantly, with severe cases enriched in males (p < 0.05). See Extended Data Table 1. Three omics approaches were applied depending on sample suitability: 16S/ITS amplicon sequencing (542 samples), shotgun RNA metagenomics (495), and shotgun DNA metagenomics (113). In total, 101 samples underwent all three methods, and 488 at least two (Fig. 1 ). Amplicon sequencing assessed bacterial/fungal diversity; RNA metatranscriptomics recovered viral genomes, functional profiles, and host partial transcriptomes; DNA metagenomics reconstructed bacterial metagenome-assembled genome (MAG) and profiled virulence, and antimicrobial resistance (AMR). 2.2. SARS-CoV-2 diversity assessed by metatranscriptomics RNA metatranscriptomics generated on average 7.35 million paired-end reads/sample (Q30 scores > 88%). Complete/partial viral genomes were recovered from 366/495 samples (73.9%), enabling lineage assignment: Delta (138), Gamma (110), and Omicron (83). Recruitment covered the three major Brazilian waves (Fig. 2 a). Mutation analysis revealed high prevalence in ORF1ab and Spike proteins. Nsp3 accumulated mutations in Alpha/Gamma; nsp12 in Alpha/Delta. Omicron had a mean 29.4 Spike mutations versus 3.8 in Wuhan (Fig. 2 b). Hotspots localized to N-terminal domain (NTD), receptor-binding domain (RBD), heptad repeat regions (HR1 and HR2) (Fig. 2 c). Gamma carried characteristic mutations at positions 417, 484, and 50, later seen in Omicron (Fig. 2 d). Principal component analysis (PCA) showed lineage-specific clustering (Fig. 2 e). 2.3. Co-infections assessed by DNA and RNA Metagenomics Metatranscriptomics also detected co-infections 72/495 (13.1%). Common agents included Influenza A (n = 9), Haemophilus parainfluenzae (n = 7), and Streptococcus mitis (n = 6). Rare pathogens such as Influenza C, Vientovirus, and Human Polyomavirus 5 were also observed (Extended Data Table 2). 2.4. Naso-oropharingeal microbiome diversity Amplicon sequencing (n = 542) provided Amplicon Sequence Variants (ASVs) and relative abundances. Alpha diversity indices (Shannon, evenness, observed features, and Faith's Phylogenetic Diversity) showed no significant differences (Fig. 3 a). Beta-diversity (Bray-Curtis) revealed significant compositional differences (p < 0.05). Controls diverged from mild cases; mild vs. severe also differed, while asymptomatic comparisons were not significant. Grouping by variants of concern (VOCs) showed significant differences between non-VOCs vs. VOCs, Gamma vs. Delta/Omicron, and Delta vs. Omicron (all p < 0.05). A list of all comparisons and their associated p-values is provided Extended Data Table 3. Core microbiome analysis showed a single ASV ( Streptococcus spp .) shared by controls only. No outcome-associated core ASVs were identified (Extended Data Table 4). Shotgun DNA metagenomic composition analyses revealed a high prevalence of the genera Cutibacterium , Streptococcus , Schaalia , Mogibacterium , Gamella , and Prevotella , with slight variations observed among different outcome groups and antibiotic usage (Fig. 3 b). 2.5. Bacterial MAGs recovered by DNA metagenomics Shotgun DNA metagenomic (n = 113) yielded on average 21.36 million paired-end reads/sample (Q30 > 89%). From 64 samples, 190 metagenome-assembled genomes (MAGs) were recovered: 55 species, 31 genera, and 20 families. Genome sizes ranged 0.58–6.9 Mbp (median 1.6 Mbp). Figure 3 c shows average nucleotide identity (ANI) vs. genome length with family assignment. MAGs completeness ranged 40%–100% (median 68%); 59 genomes exceeded 90% completeness. Contamination was low (median 1.75%, all < 5%). Median ANI to GTDB genomes was 96.6%, indicating close similarity. Thirteen high-quality MAGs (< 96% ANI, ≥ 90% completeness, ≤ 5% contamination) may represent novel species (e.g., Abiotrophia, Actinomyces, Corynebacterium, Granulicatella, Porphyromonas, Rothia, and Streptococcus). A phylogenetic tree of MAGs is shown in Fig. 3 d. Complete description of each MAG is provided as Supplementary Table S1 . 2.6. Functional Landscape of the Microbiome Functional annotation of RNA data indicated a diverse representation of the COG categories (Fig. 4 a). COG analysis revealed significant differences in normalized counts for ten out of nineteen functional categories, with category I (Lipid transport and metabolism) demonstrating significant results concerning outcomes (mild vs. moderate) and hospitalization status (Fig. 4 b). Categories I, S, and Q are depicted in Fig. 4 c. A detailed comparison and effect size is provided in Extended Data Table 5. Contigs were also annotated for Protein Families (Pfam), Enzyme Commission (E.C.), Carbohydrate-Active enzymes (CAZy), and Kyoto Encyclopedia of Genes and Genomes (KEGG), with the top ten identifiers for each source annotation shown in Fig. 4 d. Pfam-conserved domains reflected two main trends: 1. Transport and Nutrient Uptake (ABC_tran, BPD_transp_1, SBP_bac_5) and 2. Protein Synthesis and RNA Processing (GTP_EFTU, MMR_HSR1, S1, tRNA-synt_1). Enzymes classification showed predominance of transferases (506,980 counts), hydrolases (299,026), and oxidoreductases (246,892). CAZy annotation highlighted Glycoside Hydrolase family 13 (GH13), which breaks down starches and glycogen, and carbohydrate-binding module family 48, typically found in enzymes that modify or degrade glycopolymers. KEGG counts and abundance analysis revealed 93 orthologs groups significantly differing across mild and moderate outcomes (Fig. 4 e). Noteworthy, top three KOs are show in Fig. 4 f. See Supplementary Table S2 for datailed comparisons and effect size. Additionally, the presence of known toxins as classified by the Database of Bacterial Exotoxin for Human (DBETH) was investigated, yelding interesting findings. Severe cases did not present exotoxins as expected, whereas mild cases presented a cluster of cytolysin produced by Streptococcus spp . Common respiratory microbiota species, such as Vellonella spp., Neisseria spp., and Staphylococcus aureus , were also found to present exotoxin hits (Fig. 5 ). 2.7. Antibiotics therapy and resistance genes A total of 2,514 matches met the criteria for inclusion in any tier, with metagenomic data generating three times more matches (1904) than metatranscriptomic data (610 matches). Henceforth, only matches included in moderate/high-confidence tiers were considered, which amounted to 808 ARG identifications (630 in metagenomics and 178 in metatranscriptomics). A total of 183 samples (36%) exhibited at least one ARG. Of these, 101 were unique to metatranscriptomic data, 68 to metagenomic data, and 14 were identified by both techniques. Notably, the same gene was identified in half of the samples analyzed. The results suggest that utilizing both metagenomic and metatranscriptomic techniques may be advantageous for studies, as our analysis indicates that these techniques are complementary. Staphylococcus accounted for the most ARGs (251), followed by Streptococcus (65) and Corynebacterium (37). Species-level analysis highlighted Staphylococcus aureus (138) and Staphylococcus epidermis (40). Metatranscriptomic ARGs were most frequent in Streptococcus mitis (16) and Haemophilus parainfluenzae (14). ARGs were also founded in commensal genera ( Rothia , Veillonella, Prevotella ). Figure 6 a presents a heatmap of absolute counts for specific ARGs (right), along with the highest taxonomic group of the contig where the ORF was identified (bottom). To the left, the dataset originating from the contig, clinical outcome of the patient, and antibiotic usage are also provided as annotations. In the figure, a cluster of high counts of Staphylococcus and Streptococcus species is observed. Notably, two samples linked to a phage contig ( Salasvirus genus) presented the same 100% identity match to a tet(C) ARG, which alone can confer resistance to tetracycline. Antibiotic-treated patients had significantly higher ARG counts (mean 10.8) than untreated (5.7; p < 0.05). 2.8. Host Gene Expression Metatranscriptomics analysis identified 274 differentially expressed genes (DEGs; p < 0.05), being 226 for control cases, 38 for asymptomatic, 19 for moderate, and 6 for severe outcomes (Fig. 6 b). Severe-specific DEGs included SEMA3E, TET2, SPART, and DAB1. Many DEGS were linked to immune response (CXCL9/10/11, TNFSF4/8, IFIT3), antiviral pathways (RSAD2, SAMD9L, DDX60), cell cycle/apoptosis (BCL2, NEK7, CENPU), protein synthesis (EIF4B, DDX46, RPL12), and signaling (MAGI2, ROCK1P1, ARHGEF7). Matrix remodeling genes (MMP8/16, ADAMTS6) and neural-related-genes (SCN1A, NLGN1, GRID2) also varied. Several uncharacterized LINC/LOC genes emerged. Overall, gene expression patterns reflect interplay of immune, metabolic signaling, and tissue responses in COVID-19. For a comprehensive listing of all 226 DEGs, including pairwise comparisons, and log2 fold-change values with associated statistical measures, refer to Supplementary Table S3 . 2.9. Microbial Transcripts and Clinical Progression Metatranscriptomic profiling identified 67 species differing across COVID-19 outcomes (Kruskal-Wallis, robust λ factor). Severe cases were enriched for Haemophilus spp., Acinetobacter spp., and Pseudomonas spp. Gut microbiota associated taxa ( Campylobacter concisus , Fusobacterium nucleatum , and Clostridium spp.) also varied, suggesting systemic dysbiosis. Respiratory commensals ( Veillonella spp., Gemella spp.) were associated with specific outcomes (Fig. 6 c). Additional comparisons revealed 130 species differing in COVID-19 vs. controls, 69 in hospitalization, 110 in intensive care unit (ICU) admission, and 34 in severity grouping. Genus level results were weaker but still highlighted Haemophilus spp. Acinetobacter spp. Veillonella spp. and Gemella spp., as key taxa. 2.10. Taxonomic Markers via 16S/ITS Genera were initially assessed for differential abundance across clinical outcomes using 16S/ITS abundance and composition tables, which yielded no significant results. However, when outcomes were re-categorized as benign (asymptomatic plus mild) versus critical (moderate plus severe), three genera exhibited significant differential abundance (p < 0.05, effect sizes ranging from − 0.194 to 0.267): Diaphorobacter , Gemella , and Haemophilus . Specifically, Gemella and Haemophilus were more abundant in benign cases, whereas Diaphorobacter was more abundant in critical cases (Fig. 6 d). The Wilcoxon rank-sum test produced well-calibrated p-values, as demonstrated by the QQ plots (Fig. 6 e). These genera may serve as taxonomic markers of COVID-19 prognosis. 3. DISCUSSION This study presents a well-structured cohort and experimental framework for exploring the clinical and molecular aspects of COVID-19. By including 561 participants, we ensured robust comparisons to assess infection-related trends. Age was critical to influence disease severity, which is consistent with reports linking advanced age to worse outcomes, likely due to diminished adaptive immunity and enhanced inflammatory responses 14 , 15 . Similarly, BMI displayed notable disparities, with moderate and severe cases showing higher averages. This corroborates prior research indicating that obesity exacerbates systemic inflammation and metabolic dysregulation in the context of SARS-CoV-2 infection 16 . Gender differences across clinical groups and hospitalization statuses also highlighted possible demographic influences on COVID-19 outcomes. However, viral load estimations did not reveal significant variations, suggesting that other host and viral dynamics may play a more decisive role. The multi-omics approach of this study, which included 16S/ITS amplicon sequencing, RNA metatranscriptomics, and DNA metagenomics, provided complementary insights into host-microbiome interactions during SARS-CoV-2 infection. Microbial diversity, viral genome recovery, functional analyses, and antimicrobial resistance profiling have collectively advanced our understanding of these interactions. Notably, the overlapping analyses of the 101 samples underscored the value of combining omics techniques to reveal intricate host-pathogen and microbe-microbe interactions that single-omics approaches alone might overlook. For example, metatranscriptomics facilitated robust reconstruction of SARS-CoV-2 lineages and identification of co-infecting pathogens, while shotgun metagenomics allowed the recovery of high-quality bacterial MAGs and detailed resistome profiling, and 16S/ITS amplicon sequencing yielded insights into bacterial and fungal community shifts correlated with disease severity. Collectively, these findings demonstrate that leveraging the strengths of diverse omics platforms in an integrated framework is invaluable for unraveling the complex biological networks modulating COVID-19 pathogenesis and offers a rich resource for identifying potential biomarkers and therapeutic targets. This study thereby sets a precedent for future investigations aiming to dissect multi-layered host-microbiome dynamics in infectious diseases through a comprehensive multi-omics lens The successful recovery of SARS-CoV-2 genomes and lineage determination through our metatranscriptomics approach highlights the temporal dynamics and prevalence of specific VOCs, aligning with global epidemiological patterns during the study period 17 , 18 . The observed mutation accumulation in Omicron, particularly in the Spike protein, reflects evolutionary adaptations for immune evasion and enhanced transmissibility, as previously reported 19 , 20 . Lineage-specific mutation patterns in ORF1ab and Spike proteins further underscore the importance of genomic hotspots, such as the RBD and NTD domains, in adaptive evolution, consistent with the selective pressures driving SARS-CoV-2 evolution 21 . The shared spike mutations at positions 417, 484, and 501 between Gamma and Omicron suggest convergent evolution, reinforcing the functional significance of these residues in viral fitness and immune evasion. The phylogenetic clustering observed in the PCA analyses confirmed the genetic divergence among VOCs, in alignment with broader genomic surveillance studies 22 . These findings emphasize the importance of genomic monitoring in anticipating shifts in variant behavior and their implications for public health. The detection of co-infections in 13.1% of the samples underscores the utility of a broad-spectrum metatranscriptomics approach in capturing the diverse etiological landscape of respiratory infections. Co-infection with pathogens such as Influenza A , Haemophilus parainfluenzae , and Streptococcus mitis was common among COVID-19 cases, consistent with studies showing that viral respiratory infections predispose patients to bacterial infections, which can exacerbate disease outcomes 23 . The presence of opportunistic pathogens suggests that immune dysregulation resulting from SARS-CoV-2 infection may facilitate such co-infections, as reported in previous studies 24 . Moreover, the identification of less common pathogens, including Influenza C , Vientovirus , and Human Polyomavirus 5, highlights the added value of metatranscriptomics in detecting underrecognized contributors to respiratory disease. This suggests that comprehensive pathogen detection strategies could significantly impact diagnosis and treatment in both epidemic and endemic settings in the future. Our analysis of bacterial and fungal communities revealed significant differences in beta diversity across clinical outcomes and SARS-CoV-2 variant classifications, although no significant differences in alpha diversity were observed. The beta diversity findings are consistent with those of studies demonstrating that disease severity can influence the microbial composition of the respiratory tract 25 , 26 . Specifically, differences in microbial communities between mild and severe cases support the hypothesis that microbial dysbiosis correlates with disease severity, potentially influencing disease progression via host-microbiome interactions and immune modulation 27 – 29 . The significant beta diversity differences based on SARS-CoV-2 variant classification suggest that different viral lineages may exert distinct pressures on microbial communities, as the virus adapts to the human host, as observed in previous studies on the impact of SARS-CoV-2 variants on the respiratory microbiome 29 . The recovery of high-quality metagenome-assembled genomes (MAGs) from the respiratory microbiome has demonstrated the power of shotgun metagenomics in elucidating microbial diversity. The observed recovery of 190 MAGs spanning 55 species and 31 genera corroborates prior studies describing microbial heterogeneity in both healthy and diseased respiratory states 25 , 30 . The identification of MAGs from clinically relevant genera, such as Streptococcus and Porphyromonas , aligns with the well-documented associations between these taxa and respiratory infections 31 . High completeness and low contamination exceeded the criteria for high-quality MAGs defined by the MIMAG standards 32 , thus ensuring the reliability of downstream taxonomic and functional inferences. The relationship between ANI values (median 96.6%) and genome completeness suggests a significant overlap with previously described species in the Genome Taxonomy Database (GTDB). However, 13 MAGs with ANI values ≤ 96%, combined with high completeness and low contamination, suggest the presence of potentially novel species within the Granulicatella and Corynebacterium genera. These findings underscore the potential of shotgun metagenomics to expand our understanding of the human microbiome, as evidenced by studies that have uncovered novel taxa in diverse environments 33 . Collectively, these results underscore the complex and largely unexplored diversity of the respiratory microbiome, paving the way for targeted studies to investigate the roles of both novel and known taxa in respiratory health and diseases. The observed functional diversity in the naso-oropharyngeal microbiome highlights its critical role in health and disease, particularly in the context of respiratory infections, such as COVID-19. Functional categories Lipid transport and metabolism (I), and Secondary metabolite biosynthesis, transport, and catabolism (Q) were more prominent in critical patients, suggesting that the naso-oropharyngeal microbiome in these individuals is undergoing metabolic reprogramming to adapt to the inflammatory and resource-limited environment of the respiratory tract during severe COVID-19. The increased activity in lipid metabolism pathways may reflect microbial competition for essential nutrients, as well as the need to exploit alternative energy sources when the host immune response alters metabolite availability and tissue homeostasis. Similar patterns have been observed in other respiratory illnesses, where increased lipid remodeling and secondary metabolite production by opportunistic pathogens correlate with more severe disease phenotypes and worse clinical outcomes 34 . The distribution and transcriptional activity of antibiotic resistance genes (ARGs) identified in our cohort underscore the complexity and ubiquity of the respiratory resistome in COVID-19 patients. Employing both metagenomics and metatranscriptomics, we observed that these approaches are mutually informative, with each uncovering unique as well as overlapping ARGs, highlighting their complementarity and supporting previous findings suggesting that integrated multi-omics is essential for a comprehensive assessment of the resistome and its potential clinical implications 35 , 36 . Importantly, the substantial prevalence of ARGs among Staphylococcus spp., particularly S. aureus , within metagenomic data is in line with prior studies identifying this genus as a significant reservoir of antibiotic resistance in the airway and hospital settings. Our ability to detect active (transcribed) ARGs in genera such as Streptococcus, Haemophilus, Rothia, Veillonella , and Prevotella points to widespread resistance mechanisms persisting across both canonical pathogens and commensals of the respiratory microbiome, echoing recent large-scale metagenomic surveys of healthy and diseased populations 37 . A statistically significant increase in ARG numbers among patients treated with antibiotics, compared to those not exposed, further reinforces clinical concerns that antibiotic therapy exerts strong selective pressure on the respiratory microbiota, expanding the resistome. The reclassification of patient outcomes into benign versus critical categories in our 16S/ITS differential abundance analysis revealed distinct taxonomic markers associated with COVID-19 prognosis. Notably, the genera Gemella and Haemophilus were significantly enriched in benign cases, while Diaphorobacter exhibited greater abundance in critical cases, highlighting a potentially protective or risk-related role for these taxa. The robust calibration of p-values and effect sizes supports the reliability of these findings, suggesting that subtle shifts in respiratory microbiota composition can serve as informative prognostic indicators. These results align with prior observations that beneficial commensals such as Gemella and Haemophilus are associated with milder respiratory disease, potentially through immune modulation or competitive exclusion of pathogenic species, whereas enrichment of genera like Diaphorobacter may reflect, or contribute to, a disrupted airway microenvironment predisposing to severe outcomes 27 , 38 . Although the overall number of differentially abundant taxa was limited, the identification of these genera underscores the importance of high-resolution microbiome profiling and refined clinical categorization strategies for biomarker discovery. Collectively, our multi-omics approach illuminates the complex interplay between host, microbiota, and disease trajectory, and suggests that defined microbial signatures when validated in larger cohorts may assist in risk stratification and guide targeted interventions in COVID-19. 4. METHODS (Online) 4.1. Study cohort, ethics and sample collection A total of 561 nasopharyngeal swab specimens were prospectively collected from patients at nine healthcare institutions across four states and three regions of Brazil, spanning the period from July 2020 to July 2022. 544 subjects were RT-PCR-confirmed positive for SARS-CoV-2, while 17 asymptomatic, SARS-CoV-2-negative individuals were included as controls. Demographic and clinical metadata (age, sex, BMI, comorbidities, COVID-19 severity, hospitalization, ICU status, antibiotic exposure, and RT-PCR cycle threshold) were systematically recorded for all participants. Informed consent was obtained in accordance with the Declaration of Helsinki, and study protocols received approval from the institutional ethics comitee by CAAE 48098621.4.0000.0071 report 4.986.409. For sample collection, two flocked swabs were used, one in each nostril, and the other in the oropharyngeal cavity, using a swab suitable for sample collection. The samples were stored in sterile tubes without any liquid until they were transported to the laboratory, where 1 mL of sterile saline was added to each tube and the samples were tested for the presence of SARS-CoV-2 by RT-PCR. The tubes containing the swabs in saline solution were then kept at -80°C for long-term storage. Prior to further processing the samples, the tubes were vortexed, and the DNA/RNA extraction was performed according to the experiment (Total RNA metatranscriptomic, metagenomic 16S and shotgun DNA metagenomic). 4.2. Nucleic acid extraction and library preparation 4.2.1. Total RNA metatranscriptomcis To extract total RNA, 140µl of nasopharyngeal/oropharyngeal swab samples were processed using the QIAamp® Viral RNA Mini kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. As an internal control, 2.5 µl of MS2 RNA phage (10 7 copies/µl; Roche Diagnostics, Mannheim, Germany) was added to the swab samples before extraction. The RNA was eluted to a final volume of 50 µl. Subsequently, 40µl of the total RNA underwent DNase I treatment and concentration using the RNA Clean & Concentrator kit (Zymo Research, Irvine, CA, USA). The concentrated RNA was then eluted in 17µl of nuclease-free water. For rRNA depletion, 1µl of QIAseq FastSelect -rRNA HMR (Qiagen) and 1µl of QIAseq FastSelect − 5S/16S/23S rRNA (Qiagen) were added to 15µl of concentrated RNA and then submitted to temperature variations according to the manufacturer's instructions. First-strand cDNA was synthesized using Superscript IV (SSIV) reverse transcriptase (Invitrogen™, Thermo Fisher Scientific Brand, Carlsbad, CA, USA). 1 µl of Sol-PrimerA (40 pmol/µl; 5’-GTTTCCCACTGGAGGATA-N9-3’) and 4 µl of depleted RNA were mixed and heated at 65°C for 5 min, followed by cooling at 4°C for at least 1 min (primer annealing step). Subsequently, 5 µl of reverse transcriptase Master Mix (comprising 2 µl 5× SSIV Buffer, 0.65 µL nuclease-free water, 1.25 µL 10 mM dNTP mix, 0.5 µL 100 mM DTT, and 0.6 µL SSIV RT) was added to the RNA and incubated at 50°C for 10 min, followed by an inactivation step at 94°C for 2 min. The second-strand DNA synthesis was carried out using Sequenase Version 2.0 DNA polymerase Master Mix (Invitrogen™). This mixture was then incubated at 37°C for 16 min, followed by an inactivation step at 94°C for 2 min (adapted from Greninger et al. 2015) 39 . For the subsequent PCR reaction, Sol-PrimerB (5’-GTTTCCCACTGGAGGATA-3’) and Platinum Taq DNA Polymerase Master Mix (Invitrogen™, Thermo Fisher Scientific Brand, Carlsbad, CA, USA) were employed. This mixture was then subjected to an incubation at 94°C for 2 min, followed by 25 cycles of 94°C for 30 s, 50°C for 45 s, and 72°C for 1 min, concluding with a final step at 72°C for 5 min (adapted from Greninger et al. 2015). The concentration of PCR product was determined using the 1X dsDNA HS Assay kit on a Qubit 4.0 fluorometer (Invitrogen™). Next-generation sequencing (NGS) libraries were generated from 100 ng of DNA (PCR product). This process utilized the Illumina DNA Prep kit (Illumina, San Diego, CA, USA), following the manufacturer's protocol. For library clean-up we used a 1X bead ratio and 30 ul of resuspension buffer to elute the library. To assess the average fragment sizes of the DNA libraries, the High Sensitivity D1000 ScreenTape Assay (Agilent, Santa Clara, CA, USA) was employed, and the analysis was performed on a 4150 TapeStation System (Agilent). The libraries were quantified utilizing the 1X dsDNA HS Assay kit on a Qubit 4.0 fluorometer (Invitrogen™), and their concentrations were adjusted to 2nM. Subsequently, the sequencing was performed on a NextSeq 550 platform (Illumina), employing a 2 x 150 bp sequencing configuration, resulting in the generation of a minimum of 3 million paired-end reads per sample. 4.2.2. Amplicon 16S/ITS metagenomics To isolate microbial genomic DNA from 250µl of nasal swab samples, we employed the DNeasy® PowerSoil® Pro Kit (Qiagen, Hilden, Germany). Once quantified, the DNA was diluted to a concentration of 1ng/mL. The 16S rRNA variable (V3 and V4) and ITS regions (ITS1 and ITS2) were enriched in the initial PCR stage using phased primers sourced from the QIAseq 16S/ITS Panel, following a round of bead purification, according to the manufacturer's instructions. Library adapters (QIAseq 16S/ITS Index kit) were incorporated during the second-stage PCR. Subsequently, a second round of bead cleanup was performed, and the resulting libraries were quantified. The quantification of these libraries was carried out using the QIAseq™ Library Quant Assay Kit for illumina (Qiagen), employing real-time PCR methodology. To assess the average fragment sizes of the 16S/ITS libraries, the D1000 ScreenTape Assay (Agilent, Santa Clara, CA, USA) was employed, and the analysis was performed on a 4150 TapeStation System (Agilent). The final denatured library concentration was 10 pM on a MiSeq (V3 kit), sequencing was carried out on a MiSeq System (Illumina) using a 2 x 300 bp sequencing configuration, a 150k paired-end reads generated per sample. 4.2.3. Shotgun DNA metagenomics Microbial genomic DNA extraction from 200µl nasal swab samples was performed using the HostZERO™ Microbial DNA Kit (Zymo Research, South America) following the manufacturer's protocol. This kit employs a methodology specifically designed to selectively lyse eukaryotic cells and degrade their DNA, thereby diminishing the presence of host DNA contaminants prior to total DNA purification. Subsequently, quantification of the extracted DNA was conducted using the 1X dsDNA HS Assay kit on a Qubit 4.0 fluorometer (Invitrogen™). The samples were then submitted to Multiple Displacement Amplification (MDA) using the REPLI-g® UltraFast Mini kit (Qiagen) to enhance the DNA quantity. Next-generation sequencing (NGS) libraries were generated from 2ng of DNA. This process utilized the Nextera XT DNA Library Preparation kit (Illumina, San Diego, CA, USA), following the manufacturer's protocol. To optimize reagent usage, the protocol employed halved reagent quantities. For library clean-up, Agencourt AMPure XP beads (Beckman Coulter, Indianapolis, IN, USA) were used at a 1X bead ratio. To assess the average fragment sizes of the DNA libraries, the High Sensitivity D1000 ScreenTape Assay (Agilent, Santa Clara, CA, USA) was employed, and the analysis was performed on a 4150 TapeStation System (Agilent). The libraries were quantified utilizing the 1X dsDNA HS Assay kit on a Qubit 4.0 fluorometer (Invitrogen™), and their concentrations were adjusted to 2nM. The sequencing was performed on the SP flow cell type using the NovaSeq 6000 platform (Illumina), employing a 2 x 150 bp sequencing configuration, which resulted in generating a minimum of 10 million paired-end reads per sample. 4.3. Bioinformatics and statistical analyses 4.3.1. Quality Control and Host Read Removal All raw data underwent quality filtering using FastQC ( https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ ) and Cutadapt 40 . Host-derived reads were subtracted by mapping to the GRCh38 human genome (with non-canonical contigs removed) using BWA 41 , followed by removal of matching reads. Low complexity reads were filtered using Prinseq 42 for DNA metagenomics and RNA metatranscriptomics, while ribosomal reads were removed using sortMeRNA 43 . Quality metrics were inspected using Varsmetagen, an in-house online platform for analysis and interpretation of NGS data from microbial experiments. 4.3.2. Microbial Profiling Amplicon data were processed with DADA2 and QIIME2 for ASV inference 44 , 45 , taxonomic assignment using SILVA (v138) and NCBI databases 46 , and abundance tables constructed. For shotgun data, Kraken2 and Bracken were used for taxonomic assignment 47 . Alpha and beta diversity metrics (Shannon, Chao1, Bray-Curtis distances) were computed using the vegan R package, with statistical significance assessed via Kruskal-Wallis and PERMANOVA as appropriate. 4.3.3. Metagenome-Assembled Genomes (MAGs) Contigs from metagenomic data were assembled with metaSPAdes 48 . MAGs were binned using Metawrap, which is a wrapper consensus tools that uses third-party binning tools 49 . MetaBAT2, CONCOCT and MaxBin were set in Metawrap and all bins were considered 50 – 52 . Quality for the bins were assessed by CheckM 53 . Taxonomic identification was performed with GTDB-Tk 54 , and completeness/contamination thresholds set at ≥ 90% and ≤ 5%, respectively, for high-quality bins. Average nucleotide identity (ANI) was calculated to infer novelty relative to GTDB genomes. A de novo maximum likelihood phylogeny tree that merges query and reference genomes to visualise evolutionary context was build using GTDB toolkit. 4.3.4. Functional and Resistome Annotation Open reading frames were predicted using Prodigal 55 , and functional annotation performed against COG, KEGG, Pfam, E.C. and CAZy databases 56 – 60 via Prokka and EggNOG-mapper 61 , 62 . Antimicrobial resistance genes were identified by searching against the CARD database 63 . Open Reading Frame (ORFs) sequences predicted in assembled contigs generated from both metagenomics and metatranscriptomics data were searched against the CARD database to elucidate the distribution of Antibiotic Resistance Genes (ARGs). Positive matches were categorized into three tiers: high-confidence identity (id) > 90% and percentage length of reference sequence (plrs) > 80%); moderate-confidence (30% 50%); and low-confidence (30% ≤ id < 60% and plrs ≥ 30. Identified ARGs were further annotated using resistance gene identifiers and linked to contig taxonomy. Metagenomic and metatranscriptomic contigs were screened using BLASTx searches against the DBETH database 64 – 66 . Predicted toxins were identified at the species and genus level by applying thresholds of > = 70% sequence identity and alignment coverage > = 50% of the database protein length. The resulting toxin-genus/species networks were visualized using Gephi. 4.3.5. Pathogen, Co-infection, and Viral Genome Analysis Complete or near-complete SARS-CoV-2 genomes were reconstructed from RNA-seq data using two approaches: Consensus mapping to NC_045512.2 reference and denovo assembly with metaSPADES. Lineage assignment was performed using PANGOLIN 67 . Mutational analyses, including entropy calculations for the spike protein, were conducted using custom scripts leveraging Biopython and the scikit-bio package. Coinfecting viruses and bacteria were detected by Kraken2/Bracken in metagenomic and metatranscriptomic data and called as true co-infections only if confirmed by clinical assessment through manual curation in Varsmetagen. Phylogenetic analyses were performed using IQ-TREE, with support assessed by 1,000 ultrafast bootstraps 68 . 4.3.6. Differential Abundance and Expression Microbiome and metatranscriptome differential analyses were performed from centered log-ratio (CLR)-transformed abundance tables using the DESeq2 R packages 69 , controlling for FDR with Benjamini-Hochberg correction. Functional and taxonomic associations with clinical covariates were tested by Mann-Whitney U or Kruskal-Wallis as appropriate, with further FDR multiple testing corrections. 4.3.7. Host Transcriptome Analysis Human reads from RNA-seq were mapped to the GRCh38 transcriptome with STAR 70 , and gene-level counts generated with featureCounts 71 . Differential expression analysis was performed in DESeq2 69 , and gene set enrichment for up/downregulated genes assessed using Enrichr (Reactome pathways), with FDR < 0.05 as the significance threshold. 4.3.8. Data Visualization and manual curation Quality control and organism detection tables and default visualizations were generated by Varsmetagen for inspection by a Physician or Laboratory analyst. Analysis and interpretation of findings was conducted following the College of American Pathologist guidelines, with two independent analysis and one review 72 . Custom visualizations (heatmaps, NMDS, PCA plots) were created with ggplot2 and seaborn 73 , 74 . 4.4. Statistical analyses Metadata variables were analyzed using ANOVA or the Chi-square test, depending on the type and distribution of the variables. Alpha diversity was compared across clinical outcome groups using the Kruskal-Wallis test, while beta diversity was assessed using PERMANOVA with multiple test corrections applied via the Benjamini-Hochberg method. Functional variables were evaluated using the Kruskal-Wallis test with Bonferroni correction. The comparative analysis of ARG between critical and benign outcome groups was conducted using the Mann-Whitney U test. Differential gene expression analysis was performed using DESEQ with Benjamini-Hochberg correction. Differential abundance of microbial transcripts was assessed using either the Kruskal-Wallis or Mann-Whitney test, as appropriate, with Benjamini-Hochberg correction. To address potential inflation of results, the λ factor was calculated, and only tests with 0.95 < λ < 1.05 were considered. Microbial taxonomic markers distinguishing benign and critical groups were identified using the Wilcoxon rank-sum test with Benjamini-Hochberg correction, with effect sizes ranging from − 0.194 to 0.267. Effect sizes were calculated to quantify the magnitude of group differences, thereby complementing statistical significance testing and facilitating interpretation of biological relevance. The prevalence test using metatranscriptomic data was conducted with Fisher’s test and Benjamini-Hochberg correction. In all analyses, a p-value or adjusted p-value threshold of 0.05 was used to determine significance. Where applicable, the robustness of the tests was evaluated using QQ-plots and other visual resources. 4.5. Data availability The sequencing data underpinning this study has been deposited in the NCBI SRA under Bioproject PRJNA1329314. The SARS-CoV-2 genomes have been submitted to GISAID, with identification details provided in Supplementary Table S4 . Processed tables and structured data have been consolidated in a multiomics data lake, with access potentially available upon request. Declarations Acknowledgements This study was supported by the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) in the call MCTIC/CNPq/FNDCT/MS/SCTIE/Decit Nº 07/2020 process 402669/2020-7 and by the MCTI Corona Omica network. We gratefully acknowledge their financial support, which made this research possible. The authors would like to express their sincere gratitude to the PROADI team for their critical support throughout the course of this study. The contributions from the staff at Laboratório clínico, Hospital Israelita Albert Einstein were instrumental in sample processing and analysis. We are also deeply appreciative of André Mario Doi, Cristóvão Luis Pitangueira Mangueira, Rubia Anita Ferraz Santana and Maiara Antonio for their technical assistance and unwavering commitment to the project. Appreciation is also extended to the technical IT and data teams of Genesis Genomics for their vital role in developing and maintaining the cloud infrastructure, which enabled seamless data integration and advanced multiomics analysis. Their expertise in infrastructure provisioning and scalable analytics greatly enhanced the quality and scope of this research. Author contributions D.A., F.M.M., T.F.A., M.C.C., and J.R.R.P. conceptualized and designed the study. F.M.M., N.M.Z.T., A.L.B., and A.P.M.S. conducted the wet-lab experiments. D.A., R.R., A.C.S.S., R.H., T.M., and E.D. carried out the bioinformatics analyses. P.H.R.S. and B.M.F. performed the statistical analyses with contributions from W.S. and D.A. The architecture for omics data integration was conceived and assembled by W.S., R.S.R., P.H.R.S., B.M.F., and D.A. F.S., A.S.P., K.G.O, R.S., A.R.O.S. and J.R.R.P. were responsible for articulating patient inclusion and sample collection. A.R.M., F.M.M., J.R.R.P and D.A. carefully discussed results and provided insights to this manuscript. The manuscript was authored by D.A., F.M.M., A.R.M., and J.R.R.P., with input from all co-authors, and all authors have reviewed and approved the final version. Competing interests The authors declare no competing interests. References Yamamoto, S. et al. The human microbiome and COVID-19: A systematic review. PLoS One 16 , e0253293 (2021). Li, Z. et al. Impact of SARS-CoV-2 infection on respiratory and gut microbiome stability: a metagenomic investigation in long-term-hospitalized COVID-19 patients. NPJ Biofilms Microbiomes 10 , 1–11 (2024). Xie, L. et al. The clinical outcome of COVID-19 is strongly associated with microbiome dynamics in the upper respiratory tract. Journal of Infection 88 , (2024). Natalini, J. G., Singh, S. & Segal, L. N. The dynamic lung microbiome in health and disease. Nat Rev Microbiol 21 , 222 (2022). Pérez-Cobas, A. E., Ginevra, C., Rusniok, C., Jarraud, S. & Buchrieser, C. The respiratory tract microbiome, the pathogen load, and clinical interventions define severity of bacterial pneumonia. Cell Rep Med 4 , 101167 (2023). von Ameln Lovison, O. et al. Unveiling the role of the upper respiratory tract microbiome in susceptibility and severity to COVID-19. Front Cell Infect Microbiol 15 , 1531084 (2025). Man, W. H., De Steenhuijsen Piters, W. A. A. & Bogaert, D. The microbiota of the respiratory tract: gatekeeper to respiratory health. Nat Rev Microbiol 15 , 259 (2017). Azim, A. et al. Peripheral airways type 2 inflammation, neutrophilia and microbial dysbiosis in severe asthma. Allergy 76 , 2070–2078 (2021). Liu, C. et al. Microbial dysbiosis and childhood asthma development: Integrated role of the airway and gut microbiome, environmental exposures, and host metabolic and immune response. Front Immunol 13 , 1028209 (2022). Song, W., Yue, Y. & Zhang, Q. Imbalance of gut microbiota is involved in the development of chronic obstructive pulmonary disease: A review. Biomedicine & Pharmacotherapy 165 , 115150 (2023). Zhu, T., Jin, J., Chen, M. & Chen, Y. The impact of infection with COVID-19 on the respiratory microbiome: A narrative review. Virulence 13 , 1076–1087 (2022). Ma, S. et al. Metagenomic analysis reveals oropharyngeal microbiota alterations in patients with COVID-19. Signal Transduct Target Ther 6 , (2021). Galeana-Cadena, D. et al. Microbiome in the nasopharynx: Insights into the impact of COVID-19 severity. Heliyon 10 , e31562 (2024). Williamson, E. J. et al. Factors associated with COVID-19-related death using OpenSAFELY. Nature 2020 584:7821 584 , 430–436 (2020). Yang, J. et al. Prevalence of comorbidities and its effects in patients infected with SARS-CoV-2: a systematic review and meta-analysis. International Journal of Infectious Diseases 94 , 91–95 (2020). Popkin, B. M. et al. Individuals with obesity and COVID-19: A global perspective on the epidemiology and biological relationships. Obesity Reviews 21 , e13128 (2020). Tegally, H. et al. Detection of a SARS-CoV-2 variant of concern in South Africa. Nature 2021 592:7854 592 , 438–443 (2021). Callaway, E. Delta coronavirus variant: scientists brace for impact. Nature 595 , 17–18 (2021). Cao, Y. et al. Omicron escapes the majority of existing SARS-CoV-2 neutralizing antibodies. Nature 2021 602:7898 602 , 657–663 (2021). Harvey, W. T. et al. SARS-CoV-2 variants, spike mutations and immune escape. Nature Reviews Microbiology 2021 19:7 19 , 409–424 (2021). Mlcochova, P. et al. SARS-CoV-2 B.1.617.2 Delta variant replication and immune evasion. Nature 599 , 114–119 (2021). O’Toole, Á. et al. Tracking the international spread of SARS-CoV-2 lineages B.1.1.7 and B.1.351/501Y-V2 with grinch. Wellcome Open Res 6 , 121 (2021). Musuuza, J. S. et al. Prevalence and outcomes of co-infection and superinfection with SARS-CoV-2 and other pathogens: A systematic review and meta-analysis. PLoS One 16 , e0251170 (2021). Manna, S., Baindara, P. & Mandal, S. M. Molecular pathogenesis of secondary bacterial infection associated to viral infections including SARS-CoV-2. J Infect Public Health 13 , 1397–1404 (2020). Merenstein, C., Bushman, F. D. & Collman, R. G. Alterations in the respiratory tract microbiome in COVID-19: current observations and potential significance. Microbiome 10 , 1–14 (2022). Xie, L. et al. The clinical outcome of COVID-19 is strongly associated with microbiome dynamics in the upper respiratory tract. Journal of Infection 88 , 106118 (2024). De Maio, F. et al. Nasopharyngeal Microbiota Profiling of SARS-CoV-2 Infected Patients. Biol Proced Online 22 , 1–4 (2020). Hoque, M. N. et al. SARS-CoV-2 infection reduces human nasopharyngeal commensal microbiome with inclusion of pathobionts. Sci Rep 11 , 1–17 (2021). Nath, S. et al. Upper respiratory tract microbiome profiles in SARS-CoV-2 Delta and Omicron infected patients exhibit variant specific patterns and robust prediction of disease groups. Microbiol Spectr 11 , (2023). Ancona, G. et al. Gut and airway microbiota dysbiosis and their role in COVID-19 and long-COVID. Front Immunol 14 , 1080043 (2023). Bellocchio, L. et al. COVID-19 on Oral Health: A New Bilateral Connection for the Pandemic. Biomedicines 2024, Vol. 12, Page 60 12 , 60 (2023). Bowers, R. M. et al. Minimum information about a single amplified genome (MISAG) and a metagenome-assembled genome (MIMAG) of bacteria and archaea. Nat Biotechnol 35 , 725–731 (2017). Pasolli, E. et al. Extensive Unexplored Human Microbiome Diversity Revealed by Over 150,000 Genomes from Metagenomes Spanning Age, Geography, and Lifestyle. Cell 176 , 649-662.e20 (2019). Maras, J. S. et al. Multi-omics analysis of respiratory specimen characterizes baseline molecular determinants associated with SARS-CoV-2 outcome. iScience 24 , (2021). Wu, J. et al. Dysbiosis of oropharyngeal microbiome and antibiotic resistance in hospitalized COVID-19 patients. J Med Virol 95 , e28727 (2023). Xu, H. et al. Pulmonary microbial spectrum of Burkholderia multivorans infection identified by metagenomic sequencing. Front Med (Lausanne) 12 , 1577363 (2025). Siew, S. W., Musa, S. M., Sabri, N. ‘Azyyati, Farida Asras, M. F. & Ahmad, H. F. The Microbiome and Metabolome Analyses of Pre-Treated Healthcare Wastes During Covid-19 Pandemic Reveal Potent Pathogens, Antibiotics Residues, and Antibiotic Resistance Genes Against Beta-Lactams. SSRN Electronic Journal (2022) doi:10.2139/SSRN.4240492. Galeeva, J. S. et al. Microbial communities of the upper respiratory tract in mild and severe COVID-19 patients: a possible link with the disease course. Frontiers in Microbiomes 2 , 1067019 (2023). Greninger, A. L. et al. Rapid metagenomic identification of viral pathogens in clinical samples by real-time nanopore sequencing analysis. Genome Med 7 , 1–13 (2015). Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J 17 , 10–12 (2011). Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics 25 , 1754–1760 (2009). Schmieder, R. & Edwards, R. Quality control and preprocessing of metagenomic datasets. Bioinformatics 27 , 863–864 (2011). Kopylova, E., Noé, L. & Touzet, H. SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data. Bioinformatics 28 , 3211–3217 (2012). Callahan, B. J. et al. DADA2: High-resolution sample inference from Illumina amplicon data. Nat Methods 13 , 581–583 (2016). Estaki, M. et al. QIIME 2 Enables Comprehensive End-to-End Analysis of Diverse Microbiome Data and Comparative Studies with Publicly Available Data. Curr Protoc Bioinformatics 70 , e100 (2020). Quast, C. et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res 41 , D590–D596 (2013). Lu, J. et al. Metagenome analysis using the Kraken software suite. Nat Protoc 17 , 2815–2839 (2022). Nurk, S., Meleshko, D., Korobeynikov, A. & Pevzner, P. A. metaSPAdes: a new versatile metagenomic assembler. Genome Res 27 , 824–834 (2017). Uritskiy, G. V., Diruggiero, J. & Taylor, J. MetaWRAP - A flexible pipeline for genome-resolved metagenomic data analysis. Microbiome 6 , 1–13 (2018). Kang, D. D. et al. MetaBAT 2: An adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies. PeerJ 2019 , e7359 (2019). Alneberg, J. et al. Binning metagenomic contigs by coverage and composition. Nature Methods 2014 11:11 11 , 1144–1146 (2014). Wu, Y. W., Simmons, B. A. & Singer, S. W. MaxBin 2.0: an automated binning algorithm to recover genomes from multiple metagenomic datasets. Bioinformatics 32 , 605–607 (2016). Parks, D. H., Imelfort, M., Skennerton, C. T., Hugenholtz, P. & Tyson, G. W. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res 25 , 1043–1055 (2015). Chaumeil, P. A., Mussig, A. J., Hugenholtz, P. & Parks, D. H. GTDB-Tk v2: memory friendly classification with the genome taxonomy database. Bioinformatics 38 , 5315–5316 (2022). Hyatt, D. et al. Prodigal: Prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics 11 , 1–11 (2010). Tatusov, R. L. et al. The COG database: An updated vesion includes eukaryotes. BMC Bioinformatics 4 , 1–14 (2003). Kanehisa, M. & Goto, S. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res 28 , 27–30 (2000). Bateman, A. et al. The Pfam protein families database. Nucleic Acids Res 32 , D138–D141 (2004). Bairoch, A. The ENZYME data bank. Nucleic Acids Res 22 , 3626–3627 (1994). Cantarel, B. I. et al. The Carbohydrate-Active EnZymes database (CAZy): an expert resource for Glycogenomics. Nucleic Acids Res 37 , D233–D238 (2009). Seemann, T. Prokka: rapid prokaryotic genome annotation. Bioinformatics 30 , 2068–2069 (2014). Cantalapiedra, C. P., Hern̗andez-Plaza, A., Letunic, I., Bork, P. & Huerta-Cepas, J. eggNOG-mapper v2: Functional Annotation, Orthology Assignments, and Domain Prediction at the Metagenomic Scale. Mol Biol Evol 38 , 5825–5829 (2021). Alcock, B. P. et al. CARD 2023: expanded curation, support for machine learning, and resistome prediction at the Comprehensive Antibiotic Resistance Database. Nucleic Acids Res 51 , D690–D699 (2023). Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. J Mol Biol 215 , 403–410 (1990). Liu, B., Zheng, D., Zhou, S., Chen, L. & Yang, J. VFDB 2022: a general classification scheme for bacterial virulence factors. Nucleic Acids Res 50 , D912–D917 (2022). Chakraborty, A., Ghosh, S., Chowdhary, G., Maulik, U. & Chakrabarti, S. DBETH: A Database of Bacterial Exotoxins for Human. Nucleic Acids Res 40 , D615–D620 (2012). O’Toole, Á. et al. Assignment of epidemiological lineages in an emerging pandemic using the pangolin tool. Virus Evol 7 , (2021). Minh, B. Q. et al. IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era. Mol Biol Evol 37 , 1530–1534 (2020). Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15 , 1–21 (2014). Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29 , 15–21 (2013). Liao, Y., Smyth, G. K. & Shi, W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30 , 923–930 (2014). Hardwick, S. A., Deveson, I. W. & Mercer, T. R. Reference standards for next-generation sequencing. Nature Reviews Genetics 2017 18:8 18 , 473–484 (2017). Wickham, H. ggplot2. Wiley Interdiscip Rev Comput Stat 3 , 180–185 (2011). Waskom, M. L. seaborn: statistical data visualization. doi:10.21105/joss.03021. Additional Declarations There is NO Competing Interest. Supplementary Files ExtendedDataTable1.docx Extended Data Table 1 ExtendedDataTable4.docx Extended data table 4 ExtendeddataTable2.docx Extended data table 2 ExtendedDataTable5.docx Extended data table 5 ExtendedDataTable3.docx Extended data table 3 SupplementaryTableS4.xlsx Supplementary Table S4 SupplementaryTableS3.xlsx Supplementary Table S3 SupplementaryTableS2.xlsx Supplementary Table S2 SupplementaryTableS1.xlsx Supplementary Table S1 Cite Share Download PDF Status: Under Review Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-8013295","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":542769730,"identity":"c90aee3b-7e34-4e5e-ad68-8ad89d40862e","order_by":0,"name":"Deyvid Amgarten","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABDUlEQVRIiWNgGAWjYDCCA0CcwAbn2gAxY+MBYrUwNjAwpIFpwloYEFoOIwRxAb7jPWYSD8rsouXdzz5/8HHHebu17YeBttTYROPSInnmjJlEwrnk3I1n0g0bZ565nbztTCJQy7G03AYcWgxu5BgbJLYx525sSGNs5m27nWx2AKiFseEwbi3334C01Odu7H/G2Py37Vyy2fmHBLTc4DF8kNh2OHe+BNAWxrYDdmY3CNgieSat8EHCueO5GySeMc7sbUtOMLsBtCUBj1/4jh/ecPBHWXXu/P40hg8/2+zszc6nP3zwocYGpxYGBg4DiAsPQLiJYJUJOJWDAPsDMCUPNdQer+JRMApGwSgYkQAAx/puZDDa5HAAAAAASUVORK5CYII=","orcid":"https://orcid.org/0000-0002-2612-5990","institution":"Faculdade Israelita de Ciências da Saúde Albert Einstein, Hospital Israelita Albert Einstein","correspondingAuthor":true,"prefix":"","firstName":"Deyvid","middleName":"","lastName":"Amgarten","suffix":""},{"id":542769731,"identity":"7384b892-b3ac-4e19-a4dd-6f5cc12caa81","order_by":1,"name":"Fernanda de Mello Malta","email":"","orcid":"https://orcid.org/0000-0001-8887-5060","institution":"Hospital Israelita Albert Einstein","correspondingAuthor":false,"prefix":"","firstName":"Fernanda","middleName":"de Mello","lastName":"Malta","suffix":""},{"id":542769732,"identity":"643ae2a2-1d95-4c27-9290-8951a9703fb4","order_by":2,"name":"Raquel Riyuzo de Almeida Franco","email":"","orcid":"","institution":"Genesis Genomics","correspondingAuthor":false,"prefix":"","firstName":"Raquel","middleName":"Riyuzo de Almeida","lastName":"Franco","suffix":""},{"id":542769733,"identity":"7e72a372-7236-458a-9747-3237fc04240d","order_by":3,"name":"Ana Laura Boechat","email":"","orcid":"","institution":"Hospital Israelita Albert Einstein","correspondingAuthor":false,"prefix":"","firstName":"Ana","middleName":"Laura","lastName":"Boechat","suffix":""},{"id":542769734,"identity":"2183f9e6-6ec9-45ae-ac68-f863156bdc3e","order_by":4,"name":"Ana Carolina dos Santos Soares","email":"","orcid":"","institution":"Genesis Genomics","correspondingAuthor":false,"prefix":"","firstName":"Ana","middleName":"Carolina dos Santos","lastName":"Soares","suffix":""},{"id":542769735,"identity":"6e53f89a-5848-481c-8159-31b2760f2c1b","order_by":5,"name":"Erick Dorlass","email":"","orcid":"","institution":"Hospital Israelita Albert Einstein","correspondingAuthor":false,"prefix":"","firstName":"Erick","middleName":"","lastName":"Dorlass","suffix":""},{"id":542769736,"identity":"2dae1119-431f-4dcc-ae9e-63468e612d32","order_by":6,"name":"Ana Moreira Salles","email":"","orcid":"","institution":"Hospital Israelita Albert Einstein","correspondingAuthor":false,"prefix":"","firstName":"Ana","middleName":"Moreira","lastName":"Salles","suffix":""},{"id":542769737,"identity":"f77b853b-6916-4bdb-8ae5-084c89b20f6b","order_by":7,"name":"Pedro Sebe Rodrigues","email":"","orcid":"","institution":"Genesis Genomics","correspondingAuthor":false,"prefix":"","firstName":"Pedro","middleName":"Sebe","lastName":"Rodrigues","suffix":""},{"id":542769738,"identity":"453c0b6e-f013-4fa9-88b2-0b47b5042f88","order_by":8,"name":"Noelly Zimpel Trecenti","email":"","orcid":"","institution":"Hospital Israelita Albert Einstein","correspondingAuthor":false,"prefix":"","firstName":"Noelly","middleName":"Zimpel","lastName":"Trecenti","suffix":""},{"id":542769739,"identity":"797e6f71-4fab-4b24-bbbb-049c301e7d94","order_by":9,"name":"Tania Mangolini","email":"","orcid":"","institution":"Genesis Genomics","correspondingAuthor":false,"prefix":"","firstName":"Tania","middleName":"","lastName":"Mangolini","suffix":""},{"id":542769740,"identity":"e9fbae01-e39b-434f-a58f-ab13a5e6ff3f","order_by":10,"name":"Bruno Martinez de Farias","email":"","orcid":"","institution":"Genesis Genomics","correspondingAuthor":false,"prefix":"","firstName":"Bruno","middleName":"Martinez","lastName":"de Farias","suffix":""},{"id":542769741,"identity":"16629a36-f39f-4306-8124-18234686f665","order_by":11,"name":"Raquel Hurtado Castillo","email":"","orcid":"","institution":"Hospital Israelita Albert Einstein","correspondingAuthor":false,"prefix":"","firstName":"Raquel","middleName":"Hurtado","lastName":"Castillo","suffix":""},{"id":542769742,"identity":"84ad0424-c060-4906-b9a6-3020f753e266","order_by":12,"name":"Rodrigo de Souza Reis","email":"","orcid":"","institution":"Genesis Genomics","correspondingAuthor":false,"prefix":"","firstName":"Rodrigo","middleName":"de Souza","lastName":"Reis","suffix":""},{"id":542769743,"identity":"f787d7c5-4965-4a51-864f-07bd092af78b","order_by":13,"name":"Welliton de Souza","email":"","orcid":"","institution":"Genesis Genomics","correspondingAuthor":false,"prefix":"","firstName":"Welliton","middleName":"","lastName":"de Souza","suffix":""},{"id":542769744,"identity":"0cd3b4fd-a386-434a-933c-2f857805cea7","order_by":14,"name":"Ketti Gleyzer de Oliveira","email":"","orcid":"","institution":"Hospital Israelita Albert Einstein","correspondingAuthor":false,"prefix":"","firstName":"Ketti","middleName":"Gleyzer","lastName":"de Oliveira","suffix":""},{"id":542769745,"identity":"f2ce7b81-b3fd-436a-a187-e3ce9d72ba01","order_by":15,"name":"Roberta Sitnik","email":"","orcid":"","institution":"Hospital Israelita Albert Einstein","correspondingAuthor":false,"prefix":"","firstName":"Roberta","middleName":"","lastName":"Sitnik","suffix":""},{"id":542769746,"identity":"ccf7cda9-0135-4646-b158-c135699f0f36","order_by":16,"name":"Fernando Spilki","email":"","orcid":"","institution":"Feevale University","correspondingAuthor":false,"prefix":"","firstName":"Fernando","middleName":"","lastName":"Spilki","suffix":""},{"id":542769747,"identity":"6813a2cd-681d-4e84-b0c1-8840ffccd08e","order_by":17,"name":"Arlene dos Santos Pinto","email":"","orcid":"","institution":"Fundação de Medicina Tropical Dr. Heitor Vieira Dourado","correspondingAuthor":false,"prefix":"","firstName":"Arlene","middleName":"dos Santos","lastName":"Pinto","suffix":""},{"id":542769748,"identity":"2a0aa5fe-a113-4b66-8136-59385eab65a0","order_by":18,"name":"Alexandre Rodrigues Marra","email":"","orcid":"","institution":"Hospital Israelita Albert Einstein","correspondingAuthor":false,"prefix":"","firstName":"Alexandre","middleName":"Rodrigues","lastName":"Marra","suffix":""},{"id":542769749,"identity":"c4b95a97-eb7b-46ae-928d-a42edff0afc4","order_by":19,"name":"Tatiana Ferreira de Almeida","email":"","orcid":"","institution":"Hospital Israelita Albert Einstein","correspondingAuthor":false,"prefix":"","firstName":"Tatiana","middleName":"Ferreira","lastName":"de Almeida","suffix":""},{"id":542769750,"identity":"f2aebd0f-80e2-4bc6-b95d-e8133d90ba19","order_by":20,"name":"Murilo Castro Cervato","email":"","orcid":"","institution":"Hospital Israelita Albert Einstein","correspondingAuthor":false,"prefix":"","firstName":"Murilo","middleName":"Castro","lastName":"Cervato","suffix":""},{"id":542769751,"identity":"cb40b647-3c3c-403b-9be3-c099c10c2665","order_by":21,"name":"Joao Renato Rebello Pinho","email":"","orcid":"https://orcid.org/0000-0003-3999-0489","institution":"Hospital Albert Einstein","correspondingAuthor":false,"prefix":"","firstName":"Joao","middleName":"Renato Rebello","lastName":"Pinho","suffix":""}],"badges":[],"createdAt":"2025-11-02 22:00:10","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-8013295/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-8013295/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":96805433,"identity":"519bd2b8-9af5-4ec3-b779-9587d0eca84b","added_by":"auto","created_at":"2025-11-26 09:09:57","extension":"jpg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":319243,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eExperimental design\u003c/strong\u003e \u003cstrong\u003eof the study\u003c/strong\u003e. Nasopharyngeal swab samples were collected from the study participants and subjected to RT-PCR testing for SARS-CoV-2. Subsequently, the samples were processed for total RNA metagenomics, total DNA metagenomics, and 16S/ITS amplicon microbiome sequencing targeting the V3 and V4 regions. The Venn diagram illustrates the overlap among the samples in which the techniques were employed. Subsequent analyses encompassed lineage and variant tracking, microbial functional profiling, and human transcriptome analyses from RNA metagenomics; MAG recovery, virulence, and multidrug resistance analyses from DNA metagenomics; and microbial diversity and abundance profiling from 16S/ITS amplicon sequencing. This multi-omics approach integrates diverse data types to provide comprehensive insights into the microbial and host landscapes.\u003c/p\u003e","description":"","filename":"Fig1.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8013295/v1/de7e799fb1e2e677753df529.jpg"},{"id":96805424,"identity":"ff8a1046-0f02-431f-9b1c-e18580d5a85b","added_by":"auto","created_at":"2025-11-26 09:09:52","extension":"jpg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":973630,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eEmergence, mutational landscapes, and evolutionary patterns of SARS-CoV-2 lineages.\u003c/strong\u003e(\u003cstrong\u003ea\u003c/strong\u003e) Temporal dynamics of major SARS-CoV-2 lineages based on the sampling date (July 2020–July 2022). Each colored trace represents the monthly counts of a given lineage. (\u003cstrong\u003eb\u003c/strong\u003e) Heatmap showing the average number of mutations per sample across key viral genes (left) and details of ORF1ab corresponding to nonstructural mature proteins (right). Rows correspond to specific genes or nsp, and columns correspond to the indicated lineages. The color intensity reflects the magnitude of the average number of mutations. (\u003cstrong\u003ec\u003c/strong\u003e) Shannon entropy (bits) for each amino acid position of the spike protein, highlighting regions of higher sequence variability. Major functional domains, including the N-terminal domain (NTD), receptor-binding domain (RBD), fusion peptide (FP), HR1, HR2, cytoplasmic tail (CT), and transmembrane domain (TM), are indicated along the x-axis. (\u003cstrong\u003ed\u003c/strong\u003e) Proportion of samples harboring mutations at each spike position within each lineage. Bar colors indicate the frequencies of position-specific changes for Wuhan, B.1.1.28, B.1.1.33, Alpha, Gamma, Delta, Omicron, and other variants. (\u003cstrong\u003ee\u003c/strong\u003e) Principal component analysis of SARS-CoV-2 whole-genome sequence variations (top) colored by lineage. PC1 and PC2, accounting for 28.4% and 22.7% of the variance, respectively, separated the major variants of the study. The scree plot (bottom) shows the explained variance of the first eight principal components.\u003c/p\u003e","description":"","filename":"Fig2.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8013295/v1/2823d4c1a7a98e1aeb4b1bfa.jpg"},{"id":96805432,"identity":"f58b39cb-47f6-4962-a923-bb7b0beaab18","added_by":"auto","created_at":"2025-11-26 09:09:56","extension":"jpg","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":1606690,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eDiversity and composition of the naso-oropharyngeal microbiome across clinical outcomes of COVID-19. \u003c/strong\u003e(\u003cstrong\u003ea\u003c/strong\u003e) \u003cstrong\u003eAlpha diversity and community structure assessed using 16S/ITS data.\u003c/strong\u003e The top panels show the Shannon (left) and Chao1 (right) diversity indices in individuals classified as Control, Asymptomatic, Mild, Moderate, or Severe; each dot represents one sample. The bottom panels show non-metric multidimensional scaling (NMDS) ordinations based on Bray–Curtis dissimilarities, colored by clinical outcome (left) or COVID-19 status (right). Dashed ellipses represent the 95% confidence intervals for each group. (\u003cstrong\u003eb\u003c/strong\u003e) \u003cstrong\u003eShotgun metagenomics–based composition heatmap.\u003c/strong\u003e Rows represent the top 50 bacterial genera (ranked by information gain in distinguishing COVID-19 outcomes), and the columns represent individual samples. Colors indicate relative abundance on a log10 scale (light to dark gold). A side color bar denotes antibiotic use (True, False, or NA).(\u003cstrong\u003ec\u003c/strong\u003e) \u003cstrong\u003eMAG genome characteristics.\u003c/strong\u003e Scatter plot of metagenome-assembled genomes (MAGs), with genome length on the x-axis and average nucleotide identity (ANI) on the y-axis. Each point represents a MAG, sized by estimated completeness and colored by the bacterial family. \u003cstrong\u003e(d) Phylogenetic analysis of recovered MAGs.\u003c/strong\u003e Maximum-likelihood phylogenetic tree of MAGs recovered in this study, combined with reference MAGs from GTDB. The outer ring highlights the major phyla (\u003cem\u003eActinobacteriota, Firmicutes, Bacteroidota, and Proteobacteria\u003c/em\u003e). The scale bar indicates one amino acid substitution per site.\u003c/p\u003e","description":"","filename":"Fig3.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8013295/v1/c47004da87bb335498cb031a.jpg"},{"id":96805467,"identity":"427b282a-6b84-4fe4-afbf-0959501f0844","added_by":"auto","created_at":"2025-11-26 09:10:02","extension":"jpg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":2143529,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFunctional gene profiles in metatranscriptomics across different COVID-19 outcomes.\u003c/strong\u003e (\u003cstrong\u003ea\u003c/strong\u003e) Heatmap of annotated gene functions using COG categories (columns) in individual samples (rows) grouped by clinical outcome. The color scale reflects the relative abundance of each function, from low (light) to high (dark red). (\u003cstrong\u003eb\u003c/strong\u003e) The volcano plot illustrates a comparison of CLR normalized counts of COG categories between mild and moderate outcomes (above), as well as hospitalization status (below). It highlights significantly differentially abundant COG categories as red dots positioned above the blue dashed line on the y-axis, with fold-change referenced against mild outcomes on the x-axis.\u003cstrong\u003e(c)\u003c/strong\u003e Box plots of CLR normalized counts are presented to compare the top three COG categories across outcomes (left) and hospitalization status (right): Lipid transport and metabolism (I), Function unknown (S), and Secondary metabolite biosynthesis, transport, and catabolism (Q). Each dot represents an individual sample; the boxes illustrate the interquartile ranges, horizontal lines denote medians, and whiskers extend to 1.5 times the interquartile range (IQR).(\u003cstrong\u003ed\u003c/strong\u003e) Bar chart of the most frequently identified functional annotations, colored by annotation database source (PFAM, EC, CAZy, and KEGG). The x-axis denotes the total number of occurrences in general, whereas the y-axis lists the annotation identifiers. (\u003cstrong\u003ee\u003c/strong\u003e) Volcano plot comparing mild versus moderate outcomes, showing significant differentially abundant KOs as red dots above the blue dashed line in the y-axis, and fold-change referenced against mild outcome in the x-axis. \u003cstrong\u003e(f) \u003c/strong\u003eBox plots of CLR normalized counts for top 3 orthologous groups: K11735, K07644 and K03325.\u003c/p\u003e","description":"","filename":"Fig4.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8013295/v1/f3095dce7fbee2188c336da7.jpg"},{"id":96917121,"identity":"0d46a7ef-5ed2-478e-96c5-9222b9f68439","added_by":"auto","created_at":"2025-11-27 14:09:17","extension":"jpg","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":774504,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eMetatranscriptomic and metagenomic detection of toxins and related virulence factors in COVID-19 patient samples across clinical outcomes. \u003c/strong\u003eToxin-taxa networks derived from metatranscriptomic and metagenomic data of COVID-19 patient samples across different clinical categories. Nodes represent toxin genes and their corresponding microbial species or genus. Node colors indicate toxin categories (e.g., cytolysins, hemolysins, cytotoxins, superantigens, and metal-binding toxins), while node size is proportional to the number of connections with other nodes.\u003c/p\u003e","description":"","filename":"Fig5.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8013295/v1/e253b2fb3efc7c4060d809cc.jpg"},{"id":96805413,"identity":"f9b0d86f-1bac-4488-bc9d-11a904609e41","added_by":"auto","created_at":"2025-11-26 09:09:47","extension":"jpg","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":302278,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eOverview of resistance gene profiles, differential expression analyses, and microbial associations with clinical outcomes\u003c/strong\u003e. (a) Heatmap of antimicrobial resistance gene (ARG) identification across multiple taxa obtained from metagenomic and metatranscriptomic data, colored by the absolute count of each ARG. The color bars on the left indicate the dataset where ARGs were identified, clinical outcomes, and antibiotic usage (see legend). (b) QQ and volcano plots of differential human gene expression in metatranscriptomics data for mild outcomes compared with other outcomes. The QQ plots (left panels) compare observed versus expected –log10(p) values to assess p-value distribution, and the volcano plots (right panels) highlight significantly up- or downregulated features (red/teal) based on fold-change and false discovery rate (FDR) thresholds. (c) Bubble chart depicting the prevalence of selected microbial species in the metatranscriptomics data across five outcome groups (Control, Asymptomatic, Mild, Moderate, and Severe). The bubble size is proportional to the prevalence; darker shades indicate higher values. All 13 taxa yielded significant prevalence as assessed by Fisher’s test (p \u0026lt; 0.05). (d) Swarm plots are presented for the three genera that demonstrated significant results in the comparative analysis of the microbiome 16S/ITS abundance tables. Positional measurements are depicted as boxplots, and CLR-normalized abundance values are plotted on the y-axis. (e) QQ plots of observed versus expected p-values from Fisher’s tests on metatranscriptomic prevalence (related to panel c) and DESeq-based differential abundance analyses (related to panel d). The diagonal line indicates the theoretical null distribution; deviations from this line suggest p-value inflation or deflation. Statistical significance in all panels was assessed using the FDR correction for multiple comparisons.\u003c/p\u003e","description":"","filename":"Fig6.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8013295/v1/d320f0608c27113e7c68892d.jpg"},{"id":96922884,"identity":"616a153f-f850-4e04-b756-17d2d4fb5d98","added_by":"auto","created_at":"2025-11-27 14:20:06","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":7831085,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-8013295/v1/ae9f1df3-2610-4efc-8abc-73311d4bcf8b.pdf"},{"id":96805436,"identity":"e27d2ae3-492a-47de-86d9-69a57fa6f984","added_by":"auto","created_at":"2025-11-26 09:09:59","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":16050,"visible":true,"origin":"","legend":"Extended Data Table 1","description":"","filename":"ExtendedDataTable1.docx","url":"https://assets-eu.researchsquare.com/files/rs-8013295/v1/d1809807da72c3231e69f003.docx"},{"id":96805469,"identity":"5f7f5fb6-b501-4021-b2ee-481bb817e82f","added_by":"auto","created_at":"2025-11-26 09:10:04","extension":"docx","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":14589,"visible":true,"origin":"","legend":"Extended data table 4","description":"","filename":"ExtendedDataTable4.docx","url":"https://assets-eu.researchsquare.com/files/rs-8013295/v1/a444ef488b0771b8ce918c35.docx"},{"id":96805463,"identity":"4a5ba19c-5c24-4c96-9ead-856d60996671","added_by":"auto","created_at":"2025-11-26 09:10:00","extension":"docx","order_by":3,"title":"","display":"","copyAsset":false,"role":"supplement","size":16817,"visible":true,"origin":"","legend":"Extended data table 2","description":"","filename":"ExtendeddataTable2.docx","url":"https://assets-eu.researchsquare.com/files/rs-8013295/v1/162bbf360208007ebb459475.docx"},{"id":96805485,"identity":"d409d097-81b8-4a01-8c7a-1f093f92883a","added_by":"auto","created_at":"2025-11-26 09:10:06","extension":"docx","order_by":4,"title":"","display":"","copyAsset":false,"role":"supplement","size":19572,"visible":true,"origin":"","legend":"Extended data table 5","description":"","filename":"ExtendedDataTable5.docx","url":"https://assets-eu.researchsquare.com/files/rs-8013295/v1/f1c61653d37e2816a9a086ba.docx"},{"id":96805428,"identity":"45a2105f-41f9-4c61-b869-0f55c171d38b","added_by":"auto","created_at":"2025-11-26 09:09:54","extension":"docx","order_by":5,"title":"","display":"","copyAsset":false,"role":"supplement","size":17348,"visible":true,"origin":"","legend":"Extended data table 3","description":"","filename":"ExtendedDataTable3.docx","url":"https://assets-eu.researchsquare.com/files/rs-8013295/v1/d5c41a71eccd701f526571e0.docx"},{"id":96805430,"identity":"bbfba8b8-218b-461c-a5b4-6363910cd7e6","added_by":"auto","created_at":"2025-11-26 09:09:55","extension":"xlsx","order_by":6,"title":"","display":"","copyAsset":false,"role":"supplement","size":19172,"visible":true,"origin":"","legend":"Supplementary Table S4","description":"","filename":"SupplementaryTableS4.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-8013295/v1/cfbb23b8af0d7f492dcfa327.xlsx"},{"id":96805434,"identity":"3e2ea926-b9be-455d-a340-6bfa9b0884b1","added_by":"auto","created_at":"2025-11-26 09:09:58","extension":"xlsx","order_by":7,"title":"","display":"","copyAsset":false,"role":"supplement","size":37341,"visible":true,"origin":"","legend":"Supplementary Table S3","description":"","filename":"SupplementaryTableS3.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-8013295/v1/37fcab6dd08462d6db01f468.xlsx"},{"id":96805465,"identity":"685654e5-9310-45ad-a795-504a26621c7d","added_by":"auto","created_at":"2025-11-26 09:10:01","extension":"xlsx","order_by":8,"title":"","display":"","copyAsset":false,"role":"supplement","size":470311,"visible":true,"origin":"","legend":"Supplementary Table S2","description":"","filename":"SupplementaryTableS2.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-8013295/v1/f42135d5f3d082500f713bb8.xlsx"},{"id":96805429,"identity":"de36145e-eb57-42e0-8192-ecd5ad40486e","added_by":"auto","created_at":"2025-11-26 09:09:55","extension":"xlsx","order_by":9,"title":"","display":"","copyAsset":false,"role":"supplement","size":115080,"visible":true,"origin":"","legend":"Supplementary Table S1","description":"","filename":"SupplementaryTableS1.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-8013295/v1/e714aef411ddcb8a0f24d4e1.xlsx"}],"financialInterests":"There is \u003cb\u003eNO\u003c/b\u003e Competing Interest.","formattedTitle":"Respiratory Microbiome Dynamics in COVID-19: A Comprehensive Multi-Omics Study","fulltext":[{"header":"1. INTRODUCTION","content":"\u003cp\u003eThe COVID-19 pandemic, driven by SARS-CoV-2, has created an unprecedented global health crisis, characterized by clinical outcomes ranging from asymptomatic infection to severe disease. While significant research has focused on viral evolution and host immune responses, the influence of the respiratory microbiome on COVID-19 severity remains understudied\u003csup\u003e\u003cspan additionalcitationids=\"CR2\" citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u003c/sup\u003e. The human respiratory tract harbors a diverse microbiome that is fundamental for maintaining health. In healthy individuals, the upper respiratory tract is predominantly colonized by oral commensal microbes, whereas the lower respiratory tract generally maintains a lower biomass\u003csup\u003e\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e\u003c/sup\u003e. The respiratory microbiome is continuously shaped by microbial migration from the upper tract and is balanced by clearance mechanisms involving mucociliary activity and innate immune responses. Disruption of this delicate equilibrium can lead to pathogenic overgrowth, contributing to respiratory illnesses, such as pneumonia and chronic lung diseases\u003csup\u003e\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e,\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eThe respiratory microbiome plays a pivotal role in modulating local and systemic immune responses. Studies have demonstrated that upper respiratory tract microbiota can significantly influence immune regulation in the lungs, a relationship that is particularly relevant in respiratory infections such as COVID-19\u003csup\u003e6,7\u003c/sup\u003e. Microbial dysbiosis, characterized by an imbalance in microbial community composition, correlates with increased inflammation and is implicated in exacerbating conditions such as asthma and chronic obstructive pulmonary disease (COPD). Factors such as infections, environmental exposures, and therapeutic interventions can perturb the microbiome, disrupt respiratory homeostasis, and potentially influence disease progression\u003csup\u003e\u003cspan additionalcitationids=\"CR9\" citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eIn COVID-19, the relationship between respiratory microbiome composition and disease severity has become an important research area. Emerging evidence suggests that specific microbiome profiles may modulate the immune responses critical for controlling SARS-CoV-2 infection\u003csup\u003e\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e,\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e\u003c/sup\u003e. Certain oropharyngeal taxa have been linked to immune activation patterns associated with varying disease severities. Additionally, SARS-CoV-2 infection may alter the respiratory microbiome, leading to dysbiosis that contributes to complications, such as bacterial pneumonia. Conversely, specific microbial communities may confer protective effects, mitigating severe COVID-19 outcomes\u003csup\u003e\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e,\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e\u003c/sup\u003e. Clarifying the role of the respiratory microbiome could provide insights into COVID-19 pathogenesis and potentially guide therapeutic strategies.\u003c/p\u003e\u003cp\u003eTo address this gap, we conducted a large-scale, multi-omics study of nasopharyngeal samples collected from 561 individuals, including 544 COVID-19 patients and 17 uninfected controls. By integrating 16S/ITS amplicon sequencing, shotgun metatranscriptomics sequencing, and shotgun metagenomics, our approach enabled us to track SARS-CoV-2 lineage dynamics, uncover co-infections, and delineate shifts in microbial diversity associated with the disease severity. In doing so, we not only recovered high-quality metagenome-assembled genomes (MAGs) and characterized functional gene profiles, but we also identified potential microbial biomarkers linked to disease progression. Our findings reveal key host\u0026ndash;microbe interactions and metabolic alterations in the respiratory tract during COVID-19, suggesting that the microbiome may play a role in modulating clinical outcomes of the disease. This study provides a comprehensive resource that lays the groundwork for future therapeutic and prognostic strategies targeting the respiratory microbiome in viral infections.\u003c/p\u003e"},{"header":"2. RESULTS","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e\u003ch2\u003e2.1. Cohort and experimental design\u003c/h2\u003e\u003cp\u003eSwab samples were collected from 561 participants between July 2020 and July 2022. Seventeen participants formed the control group (reverse transcription \u0026ndash; polymerase chain reaction (RT-PCR) negative, no respiratory symptoms), while 544 had confirmed SARS-CoV-2 infection. Age ranged from 1\u0026ndash;98 years (mean of 39.3), increasing significantly with disease severity (p\u0026thinsp;\u0026lt;\u0026thinsp;0.05): asymptomatic (33.9) vs. severe (57.5). Body Mass Index (BMI) ranged 17\u0026ndash;44, with higher means in severe (29.1) and moderate (29.3) groups than in mild (26.2) or asymptomatic (28.4) (p\u0026thinsp;\u0026lt;\u0026thinsp;0.05). Gender distribution differed significantly, with severe cases enriched in males (p\u0026thinsp;\u0026lt;\u0026thinsp;0.05). See Extended Data Table\u0026nbsp;1.\u003c/p\u003e\u003cp\u003eThree omics approaches were applied depending on sample suitability: 16S/ITS amplicon sequencing (542 samples), shotgun RNA metagenomics (495), and shotgun DNA metagenomics (113). In total, 101 samples underwent all three methods, and 488 at least two (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). Amplicon sequencing assessed bacterial/fungal diversity; RNA metatranscriptomics recovered viral genomes, functional profiles, and host partial transcriptomes; DNA metagenomics reconstructed bacterial metagenome-assembled genome (MAG) and profiled virulence, and antimicrobial resistance (AMR).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec4\" class=\"Section2\"\u003e\u003ch2\u003e2.2. SARS-CoV-2 diversity assessed by metatranscriptomics\u003c/h2\u003e\u003cp\u003eRNA metatranscriptomics generated on average 7.35\u0026nbsp;million paired-end reads/sample (Q30 scores\u0026thinsp;\u0026gt;\u0026thinsp;88%). Complete/partial viral genomes were recovered from 366/495 samples (73.9%), enabling lineage assignment: Delta (138), Gamma (110), and Omicron (83). Recruitment covered the three major Brazilian waves (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea).\u003c/p\u003e\u003cp\u003eMutation analysis revealed high prevalence in ORF1ab and Spike proteins. Nsp3 accumulated mutations in Alpha/Gamma; nsp12 in Alpha/Delta. Omicron had a mean 29.4 Spike mutations versus 3.8 in Wuhan (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eb). Hotspots localized to N-terminal domain (NTD), receptor-binding domain (RBD), heptad repeat regions (HR1 and HR2) (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ec). Gamma carried characteristic mutations at positions 417, 484, and 50, later seen in Omicron (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ed). Principal component analysis (PCA) showed lineage-specific clustering (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ee).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec5\" class=\"Section2\"\u003e\u003ch2\u003e2.3. Co-infections assessed by DNA and RNA Metagenomics\u003c/h2\u003e\u003cp\u003eMetatranscriptomics also detected co-infections 72/495 (13.1%). Common agents included Influenza A (n\u0026thinsp;=\u0026thinsp;9), \u003cem\u003eHaemophilus parainfluenzae\u003c/em\u003e (n\u0026thinsp;=\u0026thinsp;7), and \u003cem\u003eStreptococcus mitis\u003c/em\u003e (n\u0026thinsp;=\u0026thinsp;6). Rare pathogens such as Influenza C, Vientovirus, and Human Polyomavirus 5 were also observed (Extended Data Table\u0026nbsp;2).\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec6\" class=\"Section2\"\u003e\u003ch2\u003e2.4. Naso-oropharingeal microbiome diversity\u003c/h2\u003e\u003cp\u003eAmplicon sequencing (n\u0026thinsp;=\u0026thinsp;542) provided Amplicon Sequence Variants (ASVs) and relative abundances. Alpha diversity indices (Shannon, evenness, observed features, and Faith's Phylogenetic Diversity) showed no significant differences (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ea).\u003c/p\u003e\u003cp\u003eBeta-diversity (Bray-Curtis) revealed significant compositional differences (p\u0026thinsp;\u0026lt;\u0026thinsp;0.05). Controls diverged from mild cases; mild vs. severe also differed, while asymptomatic comparisons were not significant. Grouping by variants of concern (VOCs) showed significant differences between non-VOCs vs. VOCs, Gamma vs. Delta/Omicron, and Delta vs. Omicron (all p\u0026thinsp;\u0026lt;\u0026thinsp;0.05). A list of all comparisons and their associated p-values is provided Extended Data Table\u0026nbsp;3. Core microbiome analysis showed a single ASV (\u003cem\u003eStreptococcus spp\u003c/em\u003e.) shared by controls only. No outcome-associated core ASVs were identified (Extended Data Table\u0026nbsp;4).\u003c/p\u003e\u003cp\u003eShotgun DNA metagenomic composition analyses revealed a high prevalence of the genera \u003cem\u003eCutibacterium\u003c/em\u003e, \u003cem\u003eStreptococcus\u003c/em\u003e, \u003cem\u003eSchaalia\u003c/em\u003e, \u003cem\u003eMogibacterium\u003c/em\u003e, \u003cem\u003eGamella\u003c/em\u003e, and \u003cem\u003ePrevotella\u003c/em\u003e, with slight variations observed among different outcome groups and antibiotic usage (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eb).\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec7\" class=\"Section2\"\u003e\u003ch2\u003e2.5. Bacterial MAGs recovered by DNA metagenomics\u003c/h2\u003e\u003cp\u003eShotgun DNA metagenomic (n\u0026thinsp;=\u0026thinsp;113) yielded on average 21.36\u0026nbsp;million paired-end reads/sample (Q30\u0026thinsp;\u0026gt;\u0026thinsp;89%). From 64 samples, 190 metagenome-assembled genomes (MAGs) were recovered: 55 species, 31 genera, and 20 families. Genome sizes ranged 0.58\u0026ndash;6.9 Mbp (median 1.6 Mbp). Figure\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ec shows average nucleotide identity (ANI) vs. genome length with family assignment.\u003c/p\u003e\u003cp\u003eMAGs completeness ranged 40%\u0026ndash;100% (median 68%); 59 genomes exceeded 90% completeness. Contamination was low (median 1.75%, all \u0026lt;\u0026thinsp;5%). Median ANI to GTDB genomes was 96.6%, indicating close similarity. Thirteen high-quality MAGs (\u0026lt;\u0026thinsp;96% ANI, \u0026ge;\u0026thinsp;90% completeness, \u0026le;\u0026thinsp;5% contamination) may represent novel species (e.g., \u003cem\u003eAbiotrophia, Actinomyces, Corynebacterium, Granulicatella, Porphyromonas, Rothia, and Streptococcus).\u003c/em\u003e A phylogenetic tree of MAGs is shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ed. Complete description of each MAG is provided as Supplementary Table \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec8\" class=\"Section2\"\u003e\u003ch2\u003e2.6. Functional Landscape of the Microbiome\u003c/h2\u003e\u003cp\u003eFunctional annotation of RNA data indicated a diverse representation of the COG categories (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea). COG analysis revealed significant differences in normalized counts for ten out of nineteen functional categories, with category I (Lipid transport and metabolism) demonstrating significant results concerning outcomes (mild vs. moderate) and hospitalization status (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eb). Categories I, S, and Q are depicted in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ec. A detailed comparison and effect size is provided in Extended Data Table\u0026nbsp;5.\u003c/p\u003e\u003cp\u003eContigs were also annotated for Protein Families (Pfam), Enzyme Commission (E.C.), Carbohydrate-Active enzymes (CAZy), and Kyoto Encyclopedia of Genes and Genomes (KEGG), with the top ten identifiers for each source annotation shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ed. Pfam-conserved domains reflected two main trends: 1. Transport and Nutrient Uptake (ABC_tran, BPD_transp_1, SBP_bac_5) and 2. Protein Synthesis and RNA Processing (GTP_EFTU, MMR_HSR1, S1, tRNA-synt_1). Enzymes classification showed predominance of transferases (506,980 counts), hydrolases (299,026), and oxidoreductases (246,892). CAZy annotation highlighted Glycoside Hydrolase family 13 (GH13), which breaks down starches and glycogen, and carbohydrate-binding module family 48, typically found in enzymes that modify or degrade glycopolymers. KEGG counts and abundance analysis revealed 93 orthologs groups significantly differing across mild and moderate outcomes (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ee). Noteworthy, top three KOs are show in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ef. See Supplementary Table \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e for datailed comparisons and effect size. Additionally, the presence of known toxins as classified by the Database of Bacterial Exotoxin for Human (DBETH) was investigated, yelding interesting findings. Severe cases did not present exotoxins as expected, whereas mild cases presented a cluster of cytolysin produced by \u003cem\u003eStreptococcus spp\u003c/em\u003e. Common respiratory microbiota species, such as \u003cem\u003eVellonella spp., Neisseria spp., and Staphylococcus aureus\u003c/em\u003e, were also found to present exotoxin hits (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec9\" class=\"Section2\"\u003e\u003ch2\u003e2.7. Antibiotics therapy and resistance genes\u003c/h2\u003e\u003cp\u003eA total of 2,514 matches met the criteria for inclusion in any tier, with metagenomic data generating three times more matches (1904) than metatranscriptomic data (610 matches). Henceforth, only matches included in moderate/high-confidence tiers were considered, which amounted to 808 ARG identifications (630 in metagenomics and 178 in metatranscriptomics). A total of 183 samples (36%) exhibited at least one ARG. Of these, 101 were unique to metatranscriptomic data, 68 to metagenomic data, and 14 were identified by both techniques. Notably, the same gene was identified in half of the samples analyzed. The results suggest that utilizing both metagenomic and metatranscriptomic techniques may be advantageous for studies, as our analysis indicates that these techniques are complementary.\u003c/p\u003e\u003cp\u003e\u003cem\u003eStaphylococcus\u003c/em\u003e accounted for the most ARGs (251), followed by \u003cem\u003eStreptococcus\u003c/em\u003e (65) and \u003cem\u003eCorynebacterium\u003c/em\u003e (37). Species-level analysis highlighted \u003cem\u003eStaphylococcus aureus\u003c/em\u003e (138) and \u003cem\u003eStaphylococcus epidermis\u003c/em\u003e (40). Metatranscriptomic ARGs were most frequent in \u003cem\u003eStreptococcus mitis\u003c/em\u003e (16) and \u003cem\u003eHaemophilus parainfluenzae\u003c/em\u003e (14). ARGs were also founded in commensal genera (\u003cem\u003eRothia\u003c/em\u003e, \u003cem\u003eVeillonella, Prevotella\u003c/em\u003e). Figure\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ea presents a heatmap of absolute counts for specific ARGs (right), along with the highest taxonomic group of the contig where the ORF was identified (bottom). To the left, the dataset originating from the contig, clinical outcome of the patient, and antibiotic usage are also provided as annotations. In the figure, a cluster of high counts of \u003cem\u003eStaphylococcus\u003c/em\u003e and \u003cem\u003eStreptococcus\u003c/em\u003e species is observed. Notably, two samples linked to a phage contig (\u003cem\u003eSalasvirus\u003c/em\u003e genus) presented the same 100% identity match to a tet(C) ARG, which alone can confer resistance to tetracycline. Antibiotic-treated patients had significantly higher ARG counts (mean 10.8) than untreated (5.7; p\u0026thinsp;\u0026lt;\u0026thinsp;0.05).\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec10\" class=\"Section2\"\u003e\u003ch2\u003e2.8. Host Gene Expression\u003c/h2\u003e\u003cp\u003eMetatranscriptomics analysis identified 274 differentially expressed genes (DEGs; p\u0026thinsp;\u0026lt;\u0026thinsp;0.05), being 226 for control cases, 38 for asymptomatic, 19 for moderate, and 6 for severe outcomes (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eb). Severe-specific DEGs included SEMA3E, TET2, SPART, and DAB1.\u003c/p\u003e\u003cp\u003eMany DEGS were linked to immune response (CXCL9/10/11, TNFSF4/8, IFIT3), antiviral pathways (RSAD2, SAMD9L, DDX60), cell cycle/apoptosis (BCL2, NEK7, CENPU), protein synthesis (EIF4B, DDX46, RPL12), and signaling (MAGI2, ROCK1P1, ARHGEF7). Matrix remodeling genes (MMP8/16, ADAMTS6) and neural-related-genes (SCN1A, NLGN1, GRID2) also varied. Several uncharacterized LINC/LOC genes emerged. Overall, gene expression patterns reflect interplay of immune, metabolic signaling, and tissue responses in COVID-19. For a comprehensive listing of all 226 DEGs, including pairwise comparisons, and log2 fold-change values with associated statistical measures, refer to Supplementary Table \u003cspan refid=\"MOESM3\" class=\"InternalRef\"\u003eS3\u003c/span\u003e.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec11\" class=\"Section2\"\u003e\u003ch2\u003e2.9. Microbial Transcripts and Clinical Progression\u003c/h2\u003e\u003cp\u003eMetatranscriptomic profiling identified 67 species differing across COVID-19 outcomes (Kruskal-Wallis, robust λ factor). Severe cases were enriched for \u003cem\u003eHaemophilus\u003c/em\u003e spp., \u003cem\u003eAcinetobacter\u003c/em\u003e spp., and \u003cem\u003ePseudomonas\u003c/em\u003e spp. Gut microbiota associated taxa (\u003cem\u003eCampylobacter concisus\u003c/em\u003e, \u003cem\u003eFusobacterium nucleatum\u003c/em\u003e, and \u003cem\u003eClostridium\u003c/em\u003e spp.) also varied, suggesting systemic dysbiosis. Respiratory commensals (\u003cem\u003eVeillonella\u003c/em\u003e spp., \u003cem\u003eGemella\u003c/em\u003e spp.) were associated with specific outcomes (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ec).\u003c/p\u003e\u003cp\u003eAdditional comparisons revealed 130 species differing in COVID-19 vs. controls, 69 in hospitalization, 110 in intensive care unit (ICU) admission, and 34 in severity grouping. Genus level results were weaker but still highlighted \u003cem\u003eHaemophilus\u003c/em\u003e spp. \u003cem\u003eAcinetobacter\u003c/em\u003e spp. \u003cem\u003eVeillonella\u003c/em\u003e spp. and \u003cem\u003eGemella\u003c/em\u003e spp., as key taxa.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec12\" class=\"Section2\"\u003e\u003ch2\u003e2.10. Taxonomic Markers via 16S/ITS\u003c/h2\u003e\u003cp\u003eGenera were initially assessed for differential abundance across clinical outcomes using 16S/ITS abundance and composition tables, which yielded no significant results. However, when outcomes were re-categorized as benign (asymptomatic plus mild) versus critical (moderate plus severe), three genera exhibited significant differential abundance (p\u0026thinsp;\u0026lt;\u0026thinsp;0.05, effect sizes ranging from \u0026minus;\u0026thinsp;0.194 to 0.267): \u003cem\u003eDiaphorobacter\u003c/em\u003e, \u003cem\u003eGemella\u003c/em\u003e, and \u003cem\u003eHaemophilus\u003c/em\u003e. Specifically, \u003cem\u003eGemella\u003c/em\u003e and \u003cem\u003eHaemophilus\u003c/em\u003e were more abundant in benign cases, whereas \u003cem\u003eDiaphorobacter\u003c/em\u003e was more abundant in critical cases (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ed). The Wilcoxon rank-sum test produced well-calibrated p-values, as demonstrated by the QQ plots (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ee). These genera may serve as taxonomic markers of COVID-19 prognosis.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e"},{"header":"3. DISCUSSION","content":"\u003cp\u003eThis study presents a well-structured cohort and experimental framework for exploring the clinical and molecular aspects of COVID-19. By including 561 participants, we ensured robust comparisons to assess infection-related trends. Age was critical to influence disease severity, which is consistent with reports linking advanced age to worse outcomes, likely due to diminished adaptive immunity and enhanced inflammatory responses\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e,\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e\u003c/sup\u003e. Similarly, BMI displayed notable disparities, with moderate and severe cases showing higher averages. This corroborates prior research indicating that obesity exacerbates systemic inflammation and metabolic dysregulation in the context of SARS-CoV-2 infection\u003csup\u003e\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e\u003c/sup\u003e. Gender differences across clinical groups and hospitalization statuses also highlighted possible demographic influences on COVID-19 outcomes. However, viral load estimations did not reveal significant variations, suggesting that other host and viral dynamics may play a more decisive role.\u003c/p\u003e\u003cp\u003eThe multi-omics approach of this study, which included 16S/ITS amplicon sequencing, RNA metatranscriptomics, and DNA metagenomics, provided complementary insights into host-microbiome interactions during SARS-CoV-2 infection. Microbial diversity, viral genome recovery, functional analyses, and antimicrobial resistance profiling have collectively advanced our understanding of these interactions. Notably, the overlapping analyses of the 101 samples underscored the value of combining omics techniques to reveal intricate host-pathogen and microbe-microbe interactions that single-omics approaches alone might overlook. For example, metatranscriptomics facilitated robust reconstruction of SARS-CoV-2 lineages and identification of co-infecting pathogens, while shotgun metagenomics allowed the recovery of high-quality bacterial MAGs and detailed resistome profiling, and 16S/ITS amplicon sequencing yielded insights into bacterial and fungal community shifts correlated with disease severity. Collectively, these findings demonstrate that leveraging the strengths of diverse omics platforms in an integrated framework is invaluable for unraveling the complex biological networks modulating COVID-19 pathogenesis and offers a rich resource for identifying potential biomarkers and therapeutic targets. This study thereby sets a precedent for future investigations aiming to dissect multi-layered host-microbiome dynamics in infectious diseases through a comprehensive multi-omics lens\u003c/p\u003e\u003cp\u003eThe successful recovery of SARS-CoV-2 genomes and lineage determination through our metatranscriptomics approach highlights the temporal dynamics and prevalence of specific VOCs, aligning with global epidemiological patterns during the study period\u003csup\u003e\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e,\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e\u003c/sup\u003e. The observed mutation accumulation in Omicron, particularly in the Spike protein, reflects evolutionary adaptations for immune evasion and enhanced transmissibility, as previously reported\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e,\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e\u003c/sup\u003e. Lineage-specific mutation patterns in ORF1ab and Spike proteins further underscore the importance of genomic hotspots, such as the RBD and NTD domains, in adaptive evolution, consistent with the selective pressures driving SARS-CoV-2 evolution\u003csup\u003e\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e\u003c/sup\u003e. The shared spike mutations at positions 417, 484, and 501 between Gamma and Omicron suggest convergent evolution, reinforcing the functional significance of these residues in viral fitness and immune evasion. The phylogenetic clustering observed in the PCA analyses confirmed the genetic divergence among VOCs, in alignment with broader genomic surveillance studies\u003csup\u003e\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e\u003c/sup\u003e. These findings emphasize the importance of genomic monitoring in anticipating shifts in variant behavior and their implications for public health.\u003c/p\u003e\u003cp\u003eThe detection of co-infections in 13.1% of the samples underscores the utility of a broad-spectrum metatranscriptomics approach in capturing the diverse etiological landscape of respiratory infections. Co-infection with pathogens such as \u003cem\u003eInfluenza A\u003c/em\u003e, \u003cem\u003eHaemophilus parainfluenzae\u003c/em\u003e, and \u003cem\u003eStreptococcus mitis\u003c/em\u003e was common among COVID-19 cases, consistent with studies showing that viral respiratory infections predispose patients to bacterial infections, which can exacerbate disease outcomes\u003csup\u003e\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e\u003c/sup\u003e. The presence of opportunistic pathogens suggests that immune dysregulation resulting from SARS-CoV-2 infection may facilitate such co-infections, as reported in previous studies\u003csup\u003e\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e\u003c/sup\u003e. Moreover, the identification of less common pathogens, including \u003cem\u003eInfluenza C\u003c/em\u003e, \u003cem\u003eVientovirus\u003c/em\u003e, and Human Polyomavirus 5, highlights the added value of metatranscriptomics in detecting underrecognized contributors to respiratory disease. This suggests that comprehensive pathogen detection strategies could significantly impact diagnosis and treatment in both epidemic and endemic settings in the future.\u003c/p\u003e\u003cp\u003eOur analysis of bacterial and fungal communities revealed significant differences in beta diversity across clinical outcomes and SARS-CoV-2 variant classifications, although no significant differences in alpha diversity were observed. The beta diversity findings are consistent with those of studies demonstrating that disease severity can influence the microbial composition of the respiratory tract\u003csup\u003e\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e,\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e\u003c/sup\u003e. Specifically, differences in microbial communities between mild and severe cases support the hypothesis that microbial dysbiosis correlates with disease severity, potentially influencing disease progression via host-microbiome interactions and immune modulation\u003csup\u003e\u003cspan additionalcitationids=\"CR28\" citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e\u003c/sup\u003e. The significant beta diversity differences based on SARS-CoV-2 variant classification suggest that different viral lineages may exert distinct pressures on microbial communities, as the virus adapts to the human host, as observed in previous studies on the impact of SARS-CoV-2 variants on the respiratory microbiome\u003csup\u003e\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eThe recovery of high-quality metagenome-assembled genomes (MAGs) from the respiratory microbiome has demonstrated the power of shotgun metagenomics in elucidating microbial diversity. The observed recovery of 190 MAGs spanning 55 species and 31 genera corroborates prior studies describing microbial heterogeneity in both healthy and diseased respiratory states\u003csup\u003e\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e,\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e\u003c/sup\u003e. The identification of MAGs from clinically relevant genera, such as \u003cem\u003eStreptococcus\u003c/em\u003e and \u003cem\u003ePorphyromonas\u003c/em\u003e, aligns with the well-documented associations between these taxa and respiratory infections\u003csup\u003e\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e\u003c/sup\u003e. High completeness and low contamination exceeded the criteria for high-quality MAGs defined by the MIMAG standards\u003csup\u003e\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e\u003c/sup\u003e, thus ensuring the reliability of downstream taxonomic and functional inferences.\u003c/p\u003e\u003cp\u003eThe relationship between ANI values (median 96.6%) and genome completeness suggests a significant overlap with previously described species in the Genome Taxonomy Database (GTDB). However, 13 MAGs with ANI values\u0026thinsp;\u0026le;\u0026thinsp;96%, combined with high completeness and low contamination, suggest the presence of potentially novel species within the \u003cem\u003eGranulicatella\u003c/em\u003e and \u003cem\u003eCorynebacterium\u003c/em\u003e genera. These findings underscore the potential of shotgun metagenomics to expand our understanding of the human microbiome, as evidenced by studies that have uncovered novel taxa in diverse environments\u003csup\u003e\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e\u003c/sup\u003e. Collectively, these results underscore the complex and largely unexplored diversity of the respiratory microbiome, paving the way for targeted studies to investigate the roles of both novel and known taxa in respiratory health and diseases.\u003c/p\u003e\u003cp\u003eThe observed functional diversity in the naso-oropharyngeal microbiome highlights its critical role in health and disease, particularly in the context of respiratory infections, such as COVID-19. Functional categories Lipid transport and metabolism (I), and Secondary metabolite biosynthesis, transport, and catabolism (Q) were more prominent in critical patients, suggesting that the naso-oropharyngeal microbiome in these individuals is undergoing metabolic reprogramming to adapt to the inflammatory and resource-limited environment of the respiratory tract during severe COVID-19. The increased activity in lipid metabolism pathways may reflect microbial competition for essential nutrients, as well as the need to exploit alternative energy sources when the host immune response alters metabolite availability and tissue homeostasis. Similar patterns have been observed in other respiratory illnesses, where increased lipid remodeling and secondary metabolite production by opportunistic pathogens correlate with more severe disease phenotypes and worse clinical outcomes\u003csup\u003e\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eThe distribution and transcriptional activity of antibiotic resistance genes (ARGs) identified in our cohort underscore the complexity and ubiquity of the respiratory resistome in COVID-19 patients. Employing both metagenomics and metatranscriptomics, we observed that these approaches are mutually informative, with each uncovering unique as well as overlapping ARGs, highlighting their complementarity and supporting previous findings suggesting that integrated multi-omics is essential for a comprehensive assessment of the resistome and its potential clinical implications\u003csup\u003e\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e,\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e\u003c/sup\u003e. Importantly, the substantial prevalence of ARGs among \u003cem\u003eStaphylococcus\u003c/em\u003e spp., particularly \u003cem\u003eS. aureus\u003c/em\u003e, within metagenomic data is in line with prior studies identifying this genus as a significant reservoir of antibiotic resistance in the airway and hospital settings. Our ability to detect active (transcribed) ARGs in genera such as \u003cem\u003eStreptococcus, Haemophilus, Rothia, Veillonella\u003c/em\u003e, and \u003cem\u003ePrevotella\u003c/em\u003e points to widespread resistance mechanisms persisting across both canonical pathogens and commensals of the respiratory microbiome, echoing recent large-scale metagenomic surveys of healthy and diseased populations\u003csup\u003e\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e\u003c/sup\u003e. A statistically significant increase in ARG numbers among patients treated with antibiotics, compared to those not exposed, further reinforces clinical concerns that antibiotic therapy exerts strong selective pressure on the respiratory microbiota, expanding the resistome.\u003c/p\u003e\u003cp\u003eThe reclassification of patient outcomes into benign versus critical categories in our 16S/ITS differential abundance analysis revealed distinct taxonomic markers associated with COVID-19 prognosis. Notably, the genera \u003cem\u003eGemella\u003c/em\u003e and \u003cem\u003eHaemophilus\u003c/em\u003e were significantly enriched in benign cases, while \u003cem\u003eDiaphorobacter\u003c/em\u003e exhibited greater abundance in critical cases, highlighting a potentially protective or risk-related role for these taxa. The robust calibration of p-values and effect sizes supports the reliability of these findings, suggesting that subtle shifts in respiratory microbiota composition can serve as informative prognostic indicators. These results align with prior observations that beneficial commensals such as \u003cem\u003eGemella\u003c/em\u003e and \u003cem\u003eHaemophilus\u003c/em\u003e are associated with milder respiratory disease, potentially through immune modulation or competitive exclusion of pathogenic species, whereas enrichment of genera like \u003cem\u003eDiaphorobacter\u003c/em\u003e may reflect, or contribute to, a disrupted airway microenvironment predisposing to severe outcomes\u003csup\u003e\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e,\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e\u003c/sup\u003e. Although the overall number of differentially abundant taxa was limited, the identification of these genera underscores the importance of high-resolution microbiome profiling and refined clinical categorization strategies for biomarker discovery. Collectively, our multi-omics approach illuminates the complex interplay between host, microbiota, and disease trajectory, and suggests that defined microbial signatures when validated in larger cohorts may assist in risk stratification and guide targeted interventions in COVID-19.\u003c/p\u003e"},{"header":"4. METHODS (Online)","content":"\u003cdiv id=\"Sec15\" class=\"Section2\"\u003e\u003ch2\u003e4.1. Study cohort, ethics and sample collection\u003c/h2\u003e\u003cp\u003eA total of 561 nasopharyngeal swab specimens were prospectively collected from patients at nine healthcare institutions across four states and three regions of Brazil, spanning the period from July 2020 to July 2022. 544 subjects were RT-PCR-confirmed positive for SARS-CoV-2, while 17 asymptomatic, SARS-CoV-2-negative individuals were included as controls. Demographic and clinical metadata (age, sex, BMI, comorbidities, COVID-19 severity, hospitalization, ICU status, antibiotic exposure, and RT-PCR cycle threshold) were systematically recorded for all participants. Informed consent was obtained in accordance with the Declaration of Helsinki, and study protocols received approval from the institutional ethics comitee by CAAE 48098621.4.0000.0071 report 4.986.409.\u003c/p\u003e\u003cp\u003eFor sample collection, two flocked swabs were used, one in each nostril, and the other in the oropharyngeal cavity, using a swab suitable for sample collection. The samples were stored in sterile tubes without any liquid until they were transported to the laboratory, where 1 mL of sterile saline was added to each tube and the samples were tested for the presence of SARS-CoV-2 by RT-PCR. The tubes containing the swabs in saline solution were then kept at -80\u0026deg;C for long-term storage. Prior to further processing the samples, the tubes were vortexed, and the DNA/RNA extraction was performed according to the experiment (Total RNA metatranscriptomic, metagenomic 16S and shotgun DNA metagenomic).\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec16\" class=\"Section2\"\u003e\u003ch2\u003e4.2. Nucleic acid extraction and library preparation\u003c/h2\u003e\u003cdiv id=\"Sec17\" class=\"Section3\"\u003e\u003ch2\u003e4.2.1. Total RNA metatranscriptomcis\u003c/h2\u003e\u003cp\u003eTo extract total RNA, 140\u0026micro;l of nasopharyngeal/oropharyngeal swab samples were processed using the QIAamp\u0026reg; Viral RNA Mini kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. As an internal control, 2.5 \u0026micro;l of MS2 RNA phage (10\u003csup\u003e7\u003c/sup\u003e copies/\u0026micro;l; Roche Diagnostics, Mannheim, Germany) was added to the swab samples before extraction. The RNA was eluted to a final volume of 50 \u0026micro;l. Subsequently, 40\u0026micro;l of the total RNA underwent DNase I treatment and concentration using the RNA Clean \u0026amp; Concentrator kit (Zymo Research, Irvine, CA, USA). The concentrated RNA was then eluted in 17\u0026micro;l of nuclease-free water. For rRNA depletion, 1\u0026micro;l of QIAseq FastSelect -rRNA HMR (Qiagen) and 1\u0026micro;l of QIAseq FastSelect \u0026minus;\u0026thinsp;5S/16S/23S rRNA (Qiagen) were added to 15\u0026micro;l of concentrated RNA and then submitted to temperature variations according to the manufacturer's instructions.\u003c/p\u003e\u003cp\u003eFirst-strand cDNA was synthesized using Superscript IV (SSIV) reverse transcriptase (Invitrogen\u0026trade;, Thermo Fisher Scientific Brand, Carlsbad, CA, USA). 1 \u0026micro;l of Sol-PrimerA (40 pmol/\u0026micro;l; 5\u0026rsquo;-GTTTCCCACTGGAGGATA-N9-3\u0026rsquo;) and 4 \u0026micro;l of depleted RNA were mixed and heated at 65\u0026deg;C for 5 min, followed by cooling at 4\u0026deg;C for at least 1 min (primer annealing step). Subsequently, 5 \u0026micro;l of reverse transcriptase Master Mix (comprising 2 \u0026micro;l 5\u0026times; SSIV Buffer, 0.65 \u0026micro;L nuclease-free water, 1.25 \u0026micro;L 10 mM dNTP mix, 0.5 \u0026micro;L 100 mM DTT, and 0.6 \u0026micro;L SSIV RT) was added to the RNA and incubated at 50\u0026deg;C for 10 min, followed by an inactivation step at 94\u0026deg;C for 2 min.\u003c/p\u003e\u003cp\u003eThe second-strand DNA synthesis was carried out using Sequenase Version 2.0 DNA polymerase Master Mix (Invitrogen\u0026trade;). This mixture was then incubated at 37\u0026deg;C for 16 min, followed by an inactivation step at 94\u0026deg;C for 2 min (adapted from Greninger et al. 2015)\u003csup\u003e\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eFor the subsequent PCR reaction, Sol-PrimerB (5\u0026rsquo;-GTTTCCCACTGGAGGATA-3\u0026rsquo;) and Platinum Taq DNA Polymerase Master Mix (Invitrogen\u0026trade;, Thermo Fisher Scientific Brand, Carlsbad, CA, USA) were employed. This mixture was then subjected to an incubation at 94\u0026deg;C for 2 min, followed by 25 cycles of 94\u0026deg;C for 30 s, 50\u0026deg;C for 45 s, and 72\u0026deg;C for 1 min, concluding with a final step at 72\u0026deg;C for 5 min (adapted from Greninger et al. 2015). The concentration of PCR product was determined using the 1X dsDNA HS Assay kit on a Qubit 4.0 fluorometer (Invitrogen\u0026trade;).\u003c/p\u003e\u003cp\u003eNext-generation sequencing (NGS) libraries were generated from 100 ng of DNA (PCR product). This process utilized the Illumina DNA Prep kit (Illumina, San Diego, CA, USA), following the manufacturer's protocol. For library clean-up we used a 1X bead ratio and 30 ul of resuspension buffer to elute the library.\u003c/p\u003e\u003cp\u003eTo assess the average fragment sizes of the DNA libraries, the High Sensitivity D1000 ScreenTape Assay (Agilent, Santa Clara, CA, USA) was employed, and the analysis was performed on a 4150 TapeStation System (Agilent). The libraries were quantified utilizing the 1X dsDNA HS Assay kit on a Qubit 4.0 fluorometer (Invitrogen\u0026trade;), and their concentrations were adjusted to 2nM. Subsequently, the sequencing was performed on a NextSeq 550 platform (Illumina), employing a 2 x 150 bp sequencing configuration, resulting in the generation of a minimum of 3\u0026nbsp;million paired-end reads per sample.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec18\" class=\"Section3\"\u003e\u003ch2\u003e4.2.2. Amplicon 16S/ITS metagenomics\u003c/h2\u003e\u003cp\u003eTo isolate microbial genomic DNA from 250\u0026micro;l of nasal swab samples, we employed the DNeasy\u0026reg; PowerSoil\u0026reg; Pro Kit (Qiagen, Hilden, Germany). Once quantified, the DNA was diluted to a concentration of 1ng/mL. The 16S rRNA variable (V3 and V4) and ITS regions (ITS1 and ITS2) were enriched in the initial PCR stage using phased primers sourced from the QIAseq 16S/ITS Panel, following a round of bead purification, according to the manufacturer's instructions.\u003c/p\u003e\u003cp\u003eLibrary adapters (QIAseq 16S/ITS Index kit) were incorporated during the second-stage PCR. Subsequently, a second round of bead cleanup was performed, and the resulting libraries were quantified. The quantification of these libraries was carried out using the QIAseq\u0026trade; Library Quant Assay Kit for illumina (Qiagen), employing real-time PCR methodology. To assess the average fragment sizes of the 16S/ITS libraries, the D1000 ScreenTape Assay (Agilent, Santa Clara, CA, USA) was employed, and the analysis was performed on a 4150 TapeStation System (Agilent). The final denatured library concentration was 10 pM on a MiSeq (V3 kit), sequencing was carried out on a MiSeq System (Illumina) using a 2 x 300 bp sequencing configuration, a 150k paired-end reads generated per sample.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec19\" class=\"Section3\"\u003e\u003ch2\u003e4.2.3. Shotgun DNA metagenomics\u003c/h2\u003e\u003cp\u003eMicrobial genomic DNA extraction from 200\u0026micro;l nasal swab samples was performed using the HostZERO\u0026trade; Microbial DNA Kit (Zymo Research, South America) following the manufacturer's protocol. This kit employs a methodology specifically designed to selectively lyse eukaryotic cells and degrade their DNA, thereby diminishing the presence of host DNA contaminants prior to total DNA purification. Subsequently, quantification of the extracted DNA was conducted using the 1X dsDNA HS Assay kit on a Qubit 4.0 fluorometer (Invitrogen\u0026trade;). The samples were then submitted to Multiple Displacement Amplification (MDA) using the REPLI-g\u0026reg; UltraFast Mini kit (Qiagen) to enhance the DNA quantity.\u003c/p\u003e\u003cp\u003eNext-generation sequencing (NGS) libraries were generated from 2ng of DNA. This process utilized the Nextera XT DNA Library Preparation kit (Illumina, San Diego, CA, USA), following the manufacturer's protocol. To optimize reagent usage, the protocol employed halved reagent quantities. For library clean-up, Agencourt AMPure XP beads (Beckman Coulter, Indianapolis, IN, USA) were used at a 1X bead ratio.\u003c/p\u003e\u003cp\u003eTo assess the average fragment sizes of the DNA libraries, the High Sensitivity D1000 ScreenTape Assay (Agilent, Santa Clara, CA, USA) was employed, and the analysis was performed on a 4150 TapeStation System (Agilent). The libraries were quantified utilizing the 1X dsDNA HS Assay kit on a Qubit 4.0 fluorometer (Invitrogen\u0026trade;), and their concentrations were adjusted to 2nM. The sequencing was performed on the SP flow cell type using the NovaSeq 6000 platform (Illumina), employing a 2 x 150 bp sequencing configuration, which resulted in generating a minimum of 10\u0026nbsp;million paired-end reads per sample.\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv id=\"Sec20\" class=\"Section2\"\u003e\u003ch2\u003e4.3. Bioinformatics and statistical analyses\u003c/h2\u003e\u003cdiv id=\"Sec21\" class=\"Section3\"\u003e\u003ch2\u003e4.3.1. Quality Control and Host Read Removal\u003c/h2\u003e\u003cp\u003eAll raw data underwent quality filtering using FastQC (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://www.bioinformatics.babraham.ac.uk/projects/fastqc/\u003c/span\u003e\u003cspan address=\"https://www.bioinformatics.babraham.ac.uk/projects/fastqc/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) and Cutadapt\u003csup\u003e\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e\u003c/sup\u003e. Host-derived reads were subtracted by mapping to the GRCh38 human genome (with non-canonical contigs removed) using BWA\u003csup\u003e\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e\u003c/sup\u003e, followed by removal of matching reads. Low complexity reads were filtered using Prinseq\u003csup\u003e\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e\u003c/sup\u003e for DNA metagenomics and RNA metatranscriptomics, while ribosomal reads were removed using sortMeRNA\u003csup\u003e\u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e\u003c/sup\u003e. Quality metrics were inspected using Varsmetagen, an \u003cem\u003ein-house\u003c/em\u003e online platform for analysis and interpretation of NGS data from microbial experiments.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec22\" class=\"Section3\"\u003e\u003ch2\u003e4.3.2. Microbial Profiling\u003c/h2\u003e\u003cp\u003eAmplicon data were processed with DADA2 and QIIME2 for ASV inference\u003csup\u003e\u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e,\u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e45\u003c/span\u003e\u003c/sup\u003e, taxonomic assignment using SILVA (v138) and NCBI databases\u003csup\u003e\u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e\u003c/sup\u003e, and abundance tables constructed. For shotgun data, Kraken2 and Bracken were used for taxonomic assignment\u003csup\u003e\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e\u003c/sup\u003e. Alpha and beta diversity metrics (Shannon, Chao1, Bray-Curtis distances) were computed using the vegan R package, with statistical significance assessed via Kruskal-Wallis and PERMANOVA as appropriate.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec23\" class=\"Section3\"\u003e\u003ch2\u003e4.3.3. Metagenome-Assembled Genomes (MAGs)\u003c/h2\u003e\u003cp\u003eContigs from metagenomic data were assembled with metaSPAdes\u003csup\u003e\u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e48\u003c/span\u003e\u003c/sup\u003e. MAGs were binned using Metawrap, which is a wrapper consensus tools that uses third-party binning tools\u003csup\u003e\u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e49\u003c/span\u003e\u003c/sup\u003e. MetaBAT2, CONCOCT and MaxBin were set in Metawrap and all bins were considered\u003csup\u003e\u003cspan additionalcitationids=\"CR51\" citationid=\"CR50\" class=\"CitationRef\"\u003e50\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e52\u003c/span\u003e\u003c/sup\u003e. Quality for the bins were assessed by CheckM\u003csup\u003e\u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e53\u003c/span\u003e\u003c/sup\u003e. Taxonomic identification was performed with GTDB-Tk\u003csup\u003e54\u003c/sup\u003e, and completeness/contamination thresholds set at \u0026ge;\u0026thinsp;90% and \u0026le;\u0026thinsp;5%, respectively, for high-quality bins. Average nucleotide identity (ANI) was calculated to infer novelty relative to GTDB genomes. A de novo maximum likelihood phylogeny tree that merges query and reference genomes to visualise evolutionary context was build using GTDB toolkit.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec24\" class=\"Section3\"\u003e\u003ch2\u003e4.3.4. Functional and Resistome Annotation\u003c/h2\u003e\u003cp\u003eOpen reading frames were predicted using Prodigal\u003csup\u003e\u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e55\u003c/span\u003e\u003c/sup\u003e, and functional annotation performed against COG, KEGG, Pfam, E.C. and CAZy databases\u003csup\u003e\u003cspan additionalcitationids=\"CR57 CR58 CR59\" citationid=\"CR56\" class=\"CitationRef\"\u003e56\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e60\u003c/span\u003e\u003c/sup\u003e via Prokka and EggNOG-mapper\u003csup\u003e\u003cspan citationid=\"CR61\" class=\"CitationRef\"\u003e61\u003c/span\u003e,\u003cspan citationid=\"CR62\" class=\"CitationRef\"\u003e62\u003c/span\u003e\u003c/sup\u003e. Antimicrobial resistance genes were identified by searching against the CARD database \u003csup\u003e\u003cspan citationid=\"CR63\" class=\"CitationRef\"\u003e63\u003c/span\u003e\u003c/sup\u003e. Open Reading Frame (ORFs) sequences predicted in assembled contigs generated from both metagenomics and metatranscriptomics data were searched against the CARD database to elucidate the distribution of Antibiotic Resistance Genes (ARGs). Positive matches were categorized into three tiers: high-confidence identity (id)\u0026thinsp;\u0026gt;\u0026thinsp;90% and percentage length of reference sequence (plrs)\u0026thinsp;\u0026gt;\u0026thinsp;80%); moderate-confidence (30% \u0026lt; id\u0026thinsp;\u0026le;\u0026thinsp;60% and plrs\u0026thinsp;\u0026gt;\u0026thinsp;50%); and low-confidence (30% \u0026le; id\u0026thinsp;\u0026lt;\u0026thinsp;60% and plrs\u0026thinsp;\u0026ge;\u0026thinsp;30. Identified ARGs were further annotated using resistance gene identifiers and linked to contig taxonomy.\u003c/p\u003e\u003cp\u003eMetagenomic and metatranscriptomic contigs were screened using BLASTx searches against the DBETH database \u003csup\u003e\u003cspan additionalcitationids=\"CR65\" citationid=\"CR64\" class=\"CitationRef\"\u003e64\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR66\" class=\"CitationRef\"\u003e66\u003c/span\u003e\u003c/sup\u003e. Predicted toxins were identified at the species and genus level by applying thresholds of \u0026gt;\u0026thinsp;=\u0026thinsp;70% sequence identity and alignment coverage\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;50% of the database protein length. The resulting toxin-genus/species networks were visualized using Gephi.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec25\" class=\"Section3\"\u003e\u003ch2\u003e4.3.5. Pathogen, Co-infection, and Viral Genome Analysis\u003c/h2\u003e\u003cp\u003eComplete or near-complete SARS-CoV-2 genomes were reconstructed from RNA-seq data using two approaches: Consensus mapping to NC_045512.2 reference and denovo assembly with metaSPADES. Lineage assignment was performed using PANGOLIN\u003csup\u003e\u003cspan citationid=\"CR67\" class=\"CitationRef\"\u003e67\u003c/span\u003e\u003c/sup\u003e. Mutational analyses, including entropy calculations for the spike protein, were conducted using custom scripts leveraging Biopython and the scikit-bio package. Coinfecting viruses and bacteria were detected by Kraken2/Bracken in metagenomic and metatranscriptomic data and called as true co-infections only if confirmed by clinical assessment through manual curation in Varsmetagen. Phylogenetic analyses were performed using IQ-TREE, with support assessed by 1,000 ultrafast bootstraps\u003csup\u003e\u003cspan citationid=\"CR68\" class=\"CitationRef\"\u003e68\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec26\" class=\"Section3\"\u003e\u003ch2\u003e4.3.6. Differential Abundance and Expression\u003c/h2\u003e\u003cp\u003eMicrobiome and metatranscriptome differential analyses were performed from centered log-ratio (CLR)-transformed abundance tables using the DESeq2 R packages\u003csup\u003e\u003cspan citationid=\"CR69\" class=\"CitationRef\"\u003e69\u003c/span\u003e\u003c/sup\u003e, controlling for FDR with Benjamini-Hochberg correction. Functional and taxonomic associations with clinical covariates were tested by Mann-Whitney U or Kruskal-Wallis as appropriate, with further FDR multiple testing corrections.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec27\" class=\"Section3\"\u003e\u003ch2\u003e4.3.7. Host Transcriptome Analysis\u003c/h2\u003e\u003cp\u003eHuman reads from RNA-seq were mapped to the GRCh38 transcriptome with STAR\u003csup\u003e\u003cspan citationid=\"CR70\" class=\"CitationRef\"\u003e70\u003c/span\u003e\u003c/sup\u003e, and gene-level counts generated with featureCounts\u003csup\u003e\u003cspan citationid=\"CR71\" class=\"CitationRef\"\u003e71\u003c/span\u003e\u003c/sup\u003e. Differential expression analysis was performed in DESeq2\u003csup\u003e69\u003c/sup\u003e, and gene set enrichment for up/downregulated genes assessed using Enrichr (Reactome pathways), with FDR\u0026thinsp;\u0026lt;\u0026thinsp;0.05 as the significance threshold.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec28\" class=\"Section3\"\u003e\u003ch2\u003e4.3.8. Data Visualization and manual curation\u003c/h2\u003e\u003cp\u003eQuality control and organism detection tables and default visualizations were generated by Varsmetagen for inspection by a Physician or Laboratory analyst. Analysis and interpretation of findings was conducted following the College of American Pathologist guidelines, with two independent analysis and one review\u003csup\u003e\u003cspan citationid=\"CR72\" class=\"CitationRef\"\u003e72\u003c/span\u003e\u003c/sup\u003e. Custom visualizations (heatmaps, NMDS, PCA plots) were created with ggplot2 and seaborn\u003csup\u003e\u003cspan citationid=\"CR73\" class=\"CitationRef\"\u003e73\u003c/span\u003e,\u003cspan citationid=\"CR74\" class=\"CitationRef\"\u003e74\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv id=\"Sec29\" class=\"Section2\"\u003e\u003ch2\u003e4.4. Statistical analyses\u003c/h2\u003e\u003cp\u003eMetadata variables were analyzed using ANOVA or the Chi-square test, depending on the type and distribution of the variables. Alpha diversity was compared across clinical outcome groups using the Kruskal-Wallis test, while beta diversity was assessed using PERMANOVA with multiple test corrections applied via the Benjamini-Hochberg method. Functional variables were evaluated using the Kruskal-Wallis test with Bonferroni correction. The comparative analysis of ARG between critical and benign outcome groups was conducted using the Mann-Whitney U test. Differential gene expression analysis was performed using DESEQ with Benjamini-Hochberg correction. Differential abundance of microbial transcripts was assessed using either the Kruskal-Wallis or Mann-Whitney test, as appropriate, with Benjamini-Hochberg correction. To address potential inflation of results, the λ factor was calculated, and only tests with 0.95 \u0026lt; λ \u0026lt; 1.05 were considered. Microbial taxonomic markers distinguishing benign and critical groups were identified using the Wilcoxon rank-sum test with Benjamini-Hochberg correction, with effect sizes ranging from \u0026minus;\u0026thinsp;0.194 to 0.267. Effect sizes were calculated to quantify the magnitude of group differences, thereby complementing statistical significance testing and facilitating interpretation of biological relevance. The prevalence test using metatranscriptomic data was conducted with Fisher\u0026rsquo;s test and Benjamini-Hochberg correction. In all analyses, a p-value or adjusted p-value threshold of 0.05 was used to determine significance. Where applicable, the robustness of the tests was evaluated using QQ-plots and other visual resources.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec30\" class=\"Section2\"\u003e\u003ch2\u003e4.5. Data availability\u003c/h2\u003e\u003cp\u003eThe sequencing data underpinning this study has been deposited in the NCBI SRA under Bioproject PRJNA1329314. The SARS-CoV-2 genomes have been submitted to GISAID, with identification details provided in Supplementary Table \u003cspan refid=\"MOESM4\" class=\"InternalRef\"\u003eS4\u003c/span\u003e. Processed tables and structured data have been consolidated in a multiomics data lake, with access potentially available upon request.\u003c/p\u003e\u003c/div\u003e"},{"header":"Declarations","content":"\u003ch2\u003eAcknowledgements\u003c/h2\u003e\n\u003cp\u003eThis study was supported by the Conselho Nacional de Desenvolvimento Cient\u0026iacute;fico e Tecnol\u0026oacute;gico (CNPq) in the call MCTIC/CNPq/FNDCT/MS/SCTIE/Decit N\u0026ordm; 07/2020 process 402669/2020-7 and by the MCTI Corona Omica network. We gratefully acknowledge their financial support, which made this research possible. The authors would like to express their sincere gratitude to the PROADI team for their critical support throughout the course of this study. The contributions from the staff at Laborat\u0026oacute;rio cl\u0026iacute;nico, Hospital Israelita Albert Einstein were instrumental in sample processing and analysis. We are also deeply appreciative of Andr\u0026eacute; Mario Doi, Crist\u0026oacute;v\u0026atilde;o Luis Pitangueira Mangueira, Rubia Anita Ferraz Santana and Maiara Antonio for their technical assistance and unwavering commitment to the project. Appreciation is also extended to the technical IT and data teams of Genesis Genomics for their vital role in developing and maintaining the cloud infrastructure, which enabled seamless data integration and advanced multiomics analysis. Their expertise in infrastructure provisioning and scalable analytics greatly enhanced the quality and scope of this research.\u003c/p\u003e\n\u003ch2 id=\"_Toc208656507\"\u003eAuthor contributions\u003c/h2\u003e\n\u003cp\u003eD.A., F.M.M., T.F.A., M.C.C., and J.R.R.P. conceptualized and designed the study. F.M.M., N.M.Z.T., A.L.B., and A.P.M.S. conducted the wet-lab experiments. D.A., R.R., A.C.S.S., R.H., T.M., and E.D. carried out the bioinformatics analyses. P.H.R.S. and B.M.F. performed the statistical analyses with contributions from W.S. and D.A. The architecture for omics data integration was conceived and assembled by W.S., R.S.R., P.H.R.S., B.M.F., and D.A. F.S., A.S.P., K.G.O, R.S., A.R.O.S. and J.R.R.P. were responsible for articulating patient inclusion and sample collection. A.R.M., F.M.M., J.R.R.P and D.A. carefully discussed results and provided insights to this manuscript. The manuscript was authored by D.A., F.M.M., A.R.M., and J.R.R.P., with input from all co-authors, and all authors have reviewed and approved the final version.\u003c/p\u003e\n\u003ch2 id=\"_Toc208656508\"\u003eCompeting interests\u003c/h2\u003e\n\u003cp\u003eThe authors declare no competing interests.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eYamamoto, S. \u003cem\u003eet al.\u003c/em\u003e The human microbiome and COVID-19: A systematic review. \u003cem\u003ePLoS One\u003c/em\u003e \u003cstrong\u003e16\u003c/strong\u003e, e0253293 (2021).\u003c/li\u003e\n\u003cli\u003eLi, Z. \u003cem\u003eet al.\u003c/em\u003e Impact of SARS-CoV-2 infection on respiratory and gut microbiome stability: a metagenomic investigation in long-term-hospitalized COVID-19 patients. \u003cem\u003eNPJ Biofilms Microbiomes\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, 1\u0026ndash;11 (2024).\u003c/li\u003e\n\u003cli\u003eXie, L. \u003cem\u003eet al.\u003c/em\u003e The clinical outcome of COVID-19 is strongly associated with microbiome dynamics in the upper respiratory tract. \u003cem\u003eJournal of Infection\u003c/em\u003e \u003cstrong\u003e88\u003c/strong\u003e, (2024).\u003c/li\u003e\n\u003cli\u003eNatalini, J. G., Singh, S. \u0026amp; Segal, L. N. The dynamic lung microbiome in health and disease. \u003cem\u003eNat Rev Microbiol\u003c/em\u003e \u003cstrong\u003e21\u003c/strong\u003e, 222 (2022).\u003c/li\u003e\n\u003cli\u003eP\u0026eacute;rez-Cobas, A. E., Ginevra, C., Rusniok, C., Jarraud, S. \u0026amp; Buchrieser, C. The respiratory tract microbiome, the pathogen load, and clinical interventions define severity of bacterial pneumonia. \u003cem\u003eCell Rep Med\u003c/em\u003e \u003cstrong\u003e4\u003c/strong\u003e, 101167 (2023).\u003c/li\u003e\n\u003cli\u003evon Ameln Lovison, O. \u003cem\u003eet al.\u003c/em\u003e Unveiling the role of the upper respiratory tract microbiome in susceptibility and severity to COVID-19. \u003cem\u003eFront Cell Infect Microbiol\u003c/em\u003e \u003cstrong\u003e15\u003c/strong\u003e, 1531084 (2025).\u003c/li\u003e\n\u003cli\u003eMan, W. H., De Steenhuijsen Piters, W. A. A. \u0026amp; Bogaert, D. The microbiota of the respiratory tract: gatekeeper to respiratory health. \u003cem\u003eNat Rev Microbiol\u003c/em\u003e \u003cstrong\u003e15\u003c/strong\u003e, 259 (2017).\u003c/li\u003e\n\u003cli\u003eAzim, A. \u003cem\u003eet al.\u003c/em\u003e Peripheral airways type 2 inflammation, neutrophilia and microbial dysbiosis in severe asthma. \u003cem\u003eAllergy\u003c/em\u003e \u003cstrong\u003e76\u003c/strong\u003e, 2070\u0026ndash;2078 (2021).\u003c/li\u003e\n\u003cli\u003eLiu, C. \u003cem\u003eet al.\u003c/em\u003e Microbial dysbiosis and childhood asthma development: Integrated role of the airway and gut microbiome, environmental exposures, and host metabolic and immune response. \u003cem\u003eFront Immunol\u003c/em\u003e \u003cstrong\u003e13\u003c/strong\u003e, 1028209 (2022).\u003c/li\u003e\n\u003cli\u003eSong, W., Yue, Y. \u0026amp; Zhang, Q. Imbalance of gut microbiota is involved in the development of chronic obstructive pulmonary disease: A review. \u003cem\u003eBiomedicine \u0026amp; Pharmacotherapy\u003c/em\u003e \u003cstrong\u003e165\u003c/strong\u003e, 115150 (2023).\u003c/li\u003e\n\u003cli\u003eZhu, T., Jin, J., Chen, M. \u0026amp; Chen, Y. The impact of infection with COVID-19 on the respiratory microbiome: A narrative review. \u003cem\u003eVirulence\u003c/em\u003e \u003cstrong\u003e13\u003c/strong\u003e, 1076\u0026ndash;1087 (2022).\u003c/li\u003e\n\u003cli\u003eMa, S. \u003cem\u003eet al.\u003c/em\u003e Metagenomic analysis reveals oropharyngeal microbiota alterations in patients with COVID-19. \u003cem\u003eSignal Transduct Target Ther\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, (2021).\u003c/li\u003e\n\u003cli\u003eGaleana-Cadena, D. \u003cem\u003eet al.\u003c/em\u003e Microbiome in the nasopharynx: Insights into the impact of COVID-19 severity. \u003cem\u003eHeliyon\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, e31562 (2024).\u003c/li\u003e\n\u003cli\u003eWilliamson, E. J. \u003cem\u003eet al.\u003c/em\u003e Factors associated with COVID-19-related death using OpenSAFELY. \u003cem\u003eNature 2020 584:7821\u003c/em\u003e \u003cstrong\u003e584\u003c/strong\u003e, 430\u0026ndash;436 (2020).\u003c/li\u003e\n\u003cli\u003eYang, J. \u003cem\u003eet al.\u003c/em\u003e Prevalence of comorbidities and its effects in patients infected with SARS-CoV-2: a systematic review and meta-analysis. \u003cem\u003eInternational Journal of Infectious Diseases\u003c/em\u003e \u003cstrong\u003e94\u003c/strong\u003e, 91\u0026ndash;95 (2020).\u003c/li\u003e\n\u003cli\u003ePopkin, B. M. \u003cem\u003eet al.\u003c/em\u003e Individuals with obesity and COVID-19: A global perspective on the epidemiology and biological relationships. \u003cem\u003eObesity Reviews\u003c/em\u003e \u003cstrong\u003e21\u003c/strong\u003e, e13128 (2020).\u003c/li\u003e\n\u003cli\u003eTegally, H. \u003cem\u003eet al.\u003c/em\u003e Detection of a SARS-CoV-2 variant of concern in South Africa. \u003cem\u003eNature 2021 592:7854\u003c/em\u003e \u003cstrong\u003e592\u003c/strong\u003e, 438\u0026ndash;443 (2021).\u003c/li\u003e\n\u003cli\u003eCallaway, E. Delta coronavirus variant: scientists brace for impact. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e595\u003c/strong\u003e, 17\u0026ndash;18 (2021).\u003c/li\u003e\n\u003cli\u003eCao, Y. \u003cem\u003eet al.\u003c/em\u003e Omicron escapes the majority of existing SARS-CoV-2 neutralizing antibodies. \u003cem\u003eNature 2021 602:7898\u003c/em\u003e \u003cstrong\u003e602\u003c/strong\u003e, 657\u0026ndash;663 (2021).\u003c/li\u003e\n\u003cli\u003eHarvey, W. T. \u003cem\u003eet al.\u003c/em\u003e SARS-CoV-2 variants, spike mutations and immune escape. \u003cem\u003eNature Reviews Microbiology 2021 19:7\u003c/em\u003e \u003cstrong\u003e19\u003c/strong\u003e, 409\u0026ndash;424 (2021).\u003c/li\u003e\n\u003cli\u003eMlcochova, P. \u003cem\u003eet al.\u003c/em\u003e SARS-CoV-2 B.1.617.2 Delta variant replication and immune evasion. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e599\u003c/strong\u003e, 114\u0026ndash;119 (2021).\u003c/li\u003e\n\u003cli\u003eO\u0026rsquo;Toole, \u0026Aacute;. \u003cem\u003eet al.\u003c/em\u003e Tracking the international spread of SARS-CoV-2 lineages B.1.1.7 and B.1.351/501Y-V2 with grinch. \u003cem\u003eWellcome Open Res\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, 121 (2021).\u003c/li\u003e\n\u003cli\u003eMusuuza, J. S. \u003cem\u003eet al.\u003c/em\u003e Prevalence and outcomes of co-infection and superinfection with SARS-CoV-2 and other pathogens: A systematic review and meta-analysis. \u003cem\u003ePLoS One\u003c/em\u003e \u003cstrong\u003e16\u003c/strong\u003e, e0251170 (2021).\u003c/li\u003e\n\u003cli\u003eManna, S., Baindara, P. \u0026amp; Mandal, S. M. Molecular pathogenesis of secondary bacterial infection associated to viral infections including SARS-CoV-2. \u003cem\u003eJ Infect Public Health\u003c/em\u003e \u003cstrong\u003e13\u003c/strong\u003e, 1397\u0026ndash;1404 (2020).\u003c/li\u003e\n\u003cli\u003eMerenstein, C., Bushman, F. D. \u0026amp; Collman, R. G. Alterations in the respiratory tract microbiome in COVID-19: current observations and potential significance. \u003cem\u003eMicrobiome\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, 1\u0026ndash;14 (2022).\u003c/li\u003e\n\u003cli\u003eXie, L. \u003cem\u003eet al.\u003c/em\u003e The clinical outcome of COVID-19 is strongly associated with microbiome dynamics in the upper respiratory tract. \u003cem\u003eJournal of Infection\u003c/em\u003e \u003cstrong\u003e88\u003c/strong\u003e, 106118 (2024).\u003c/li\u003e\n\u003cli\u003eDe Maio, F. \u003cem\u003eet al.\u003c/em\u003e Nasopharyngeal Microbiota Profiling of SARS-CoV-2 Infected Patients. \u003cem\u003eBiol Proced Online\u003c/em\u003e \u003cstrong\u003e22\u003c/strong\u003e, 1\u0026ndash;4 (2020).\u003c/li\u003e\n\u003cli\u003eHoque, M. N. \u003cem\u003eet al.\u003c/em\u003e SARS-CoV-2 infection reduces human nasopharyngeal commensal microbiome with inclusion of pathobionts. \u003cem\u003eSci Rep\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, 1\u0026ndash;17 (2021).\u003c/li\u003e\n\u003cli\u003eNath, S. \u003cem\u003eet al.\u003c/em\u003e Upper respiratory tract microbiome profiles in SARS-CoV-2 Delta and Omicron infected patients exhibit variant specific patterns and robust prediction of disease groups. \u003cem\u003eMicrobiol Spectr\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, (2023).\u003c/li\u003e\n\u003cli\u003eAncona, G. \u003cem\u003eet al.\u003c/em\u003e Gut and airway microbiota dysbiosis and their role in COVID-19 and long-COVID. \u003cem\u003eFront Immunol\u003c/em\u003e \u003cstrong\u003e14\u003c/strong\u003e, 1080043 (2023).\u003c/li\u003e\n\u003cli\u003eBellocchio, L. \u003cem\u003eet al.\u003c/em\u003e COVID-19 on Oral Health: A New Bilateral Connection for the Pandemic. \u003cem\u003eBiomedicines 2024, Vol. 12, Page 60\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, 60 (2023).\u003c/li\u003e\n\u003cli\u003eBowers, R. M. \u003cem\u003eet al.\u003c/em\u003e Minimum information about a single amplified genome (MISAG) and a metagenome-assembled genome (MIMAG) of bacteria and archaea. \u003cem\u003eNat Biotechnol\u003c/em\u003e \u003cstrong\u003e35\u003c/strong\u003e, 725\u0026ndash;731 (2017).\u003c/li\u003e\n\u003cli\u003ePasolli, E. \u003cem\u003eet al.\u003c/em\u003e Extensive Unexplored Human Microbiome Diversity Revealed by Over 150,000 Genomes from Metagenomes Spanning Age, Geography, and Lifestyle. \u003cem\u003eCell\u003c/em\u003e \u003cstrong\u003e176\u003c/strong\u003e, 649-662.e20 (2019).\u003c/li\u003e\n\u003cli\u003eMaras, J. S. \u003cem\u003eet al.\u003c/em\u003e Multi-omics analysis of respiratory specimen characterizes baseline molecular determinants associated with SARS-CoV-2 outcome. \u003cem\u003eiScience\u003c/em\u003e \u003cstrong\u003e24\u003c/strong\u003e, (2021).\u003c/li\u003e\n\u003cli\u003eWu, J. \u003cem\u003eet al.\u003c/em\u003e Dysbiosis of oropharyngeal microbiome and antibiotic resistance in hospitalized COVID-19 patients. \u003cem\u003eJ Med Virol\u003c/em\u003e \u003cstrong\u003e95\u003c/strong\u003e, e28727 (2023).\u003c/li\u003e\n\u003cli\u003eXu, H. \u003cem\u003eet al.\u003c/em\u003e Pulmonary microbial spectrum of Burkholderia multivorans infection identified by metagenomic sequencing. \u003cem\u003eFront Med (Lausanne)\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, 1577363 (2025).\u003c/li\u003e\n\u003cli\u003eSiew, S. W., Musa, S. M., Sabri, N. \u0026lsquo;Azyyati, Farida Asras, M. F. \u0026amp; Ahmad, H. F. The Microbiome and Metabolome Analyses of Pre-Treated Healthcare Wastes During Covid-19 Pandemic Reveal Potent Pathogens, Antibiotics Residues, and Antibiotic Resistance Genes Against Beta-Lactams. \u003cem\u003eSSRN Electronic Journal\u003c/em\u003e (2022) doi:10.2139/SSRN.4240492.\u003c/li\u003e\n\u003cli\u003eGaleeva, J. S. \u003cem\u003eet al.\u003c/em\u003e Microbial communities of the upper respiratory tract in mild and severe COVID-19 patients: a possible link with the disease course. \u003cem\u003eFrontiers in Microbiomes\u003c/em\u003e \u003cstrong\u003e2\u003c/strong\u003e, 1067019 (2023).\u003c/li\u003e\n\u003cli\u003eGreninger, A. L. \u003cem\u003eet al.\u003c/em\u003e Rapid metagenomic identification of viral pathogens in clinical samples by real-time nanopore sequencing analysis. \u003cem\u003eGenome Med\u003c/em\u003e \u003cstrong\u003e7\u003c/strong\u003e, 1\u0026ndash;13 (2015).\u003c/li\u003e\n\u003cli\u003eMartin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. \u003cem\u003eEMBnet J\u003c/em\u003e \u003cstrong\u003e17\u003c/strong\u003e, 10\u0026ndash;12 (2011).\u003c/li\u003e\n\u003cli\u003eLi, H. \u0026amp; Durbin, R. Fast and accurate short read alignment with Burrows\u0026ndash;Wheeler transform. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e25\u003c/strong\u003e, 1754\u0026ndash;1760 (2009).\u003c/li\u003e\n\u003cli\u003eSchmieder, R. \u0026amp; Edwards, R. Quality control and preprocessing of metagenomic datasets. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e27\u003c/strong\u003e, 863\u0026ndash;864 (2011).\u003c/li\u003e\n\u003cli\u003eKopylova, E., No\u0026eacute;, L. \u0026amp; Touzet, H. SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e28\u003c/strong\u003e, 3211\u0026ndash;3217 (2012).\u003c/li\u003e\n\u003cli\u003eCallahan, B. J. \u003cem\u003eet al.\u003c/em\u003e DADA2: High-resolution sample inference from Illumina amplicon data. \u003cem\u003eNat Methods\u003c/em\u003e \u003cstrong\u003e13\u003c/strong\u003e, 581\u0026ndash;583 (2016).\u003c/li\u003e\n\u003cli\u003eEstaki, M. \u003cem\u003eet al.\u003c/em\u003e QIIME 2 Enables Comprehensive End-to-End Analysis of Diverse Microbiome Data and Comparative Studies with Publicly Available Data. \u003cem\u003eCurr Protoc Bioinformatics\u003c/em\u003e \u003cstrong\u003e70\u003c/strong\u003e, e100 (2020).\u003c/li\u003e\n\u003cli\u003eQuast, C. \u003cem\u003eet al.\u003c/em\u003e The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e41\u003c/strong\u003e, D590\u0026ndash;D596 (2013).\u003c/li\u003e\n\u003cli\u003eLu, J. \u003cem\u003eet al.\u003c/em\u003e Metagenome analysis using the Kraken software suite. \u003cem\u003eNat Protoc\u003c/em\u003e \u003cstrong\u003e17\u003c/strong\u003e, 2815\u0026ndash;2839 (2022).\u003c/li\u003e\n\u003cli\u003eNurk, S., Meleshko, D., Korobeynikov, A. \u0026amp; Pevzner, P. A. metaSPAdes: a new versatile metagenomic assembler. \u003cem\u003eGenome Res\u003c/em\u003e \u003cstrong\u003e27\u003c/strong\u003e, 824\u0026ndash;834 (2017).\u003c/li\u003e\n\u003cli\u003eUritskiy, G. V., Diruggiero, J. \u0026amp; Taylor, J. MetaWRAP - A flexible pipeline for genome-resolved metagenomic data analysis. \u003cem\u003eMicrobiome\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, 1\u0026ndash;13 (2018).\u003c/li\u003e\n\u003cli\u003eKang, D. D. \u003cem\u003eet al.\u003c/em\u003e MetaBAT 2: An adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies. \u003cem\u003ePeerJ\u003c/em\u003e \u003cstrong\u003e2019\u003c/strong\u003e, e7359 (2019).\u003c/li\u003e\n\u003cli\u003eAlneberg, J. \u003cem\u003eet al.\u003c/em\u003e Binning metagenomic contigs by coverage and composition. \u003cem\u003eNature Methods 2014 11:11\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, 1144\u0026ndash;1146 (2014).\u003c/li\u003e\n\u003cli\u003eWu, Y. W., Simmons, B. A. \u0026amp; Singer, S. W. MaxBin 2.0: an automated binning algorithm to recover genomes from multiple metagenomic datasets. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e32\u003c/strong\u003e, 605\u0026ndash;607 (2016).\u003c/li\u003e\n\u003cli\u003eParks, D. H., Imelfort, M., Skennerton, C. T., Hugenholtz, P. \u0026amp; Tyson, G. W. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. \u003cem\u003eGenome Res\u003c/em\u003e \u003cstrong\u003e25\u003c/strong\u003e, 1043\u0026ndash;1055 (2015).\u003c/li\u003e\n\u003cli\u003eChaumeil, P. A., Mussig, A. J., Hugenholtz, P. \u0026amp; Parks, D. H. GTDB-Tk v2: memory friendly classification with the genome taxonomy database. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e38\u003c/strong\u003e, 5315\u0026ndash;5316 (2022).\u003c/li\u003e\n\u003cli\u003eHyatt, D. \u003cem\u003eet al.\u003c/em\u003e Prodigal: Prokaryotic gene recognition and translation initiation site identification. \u003cem\u003eBMC Bioinformatics\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, 1\u0026ndash;11 (2010).\u003c/li\u003e\n\u003cli\u003eTatusov, R. L. \u003cem\u003eet al.\u003c/em\u003e The COG database: An updated vesion includes eukaryotes. \u003cem\u003eBMC Bioinformatics\u003c/em\u003e \u003cstrong\u003e4\u003c/strong\u003e, 1\u0026ndash;14 (2003).\u003c/li\u003e\n\u003cli\u003eKanehisa, M. \u0026amp; Goto, S. KEGG: Kyoto Encyclopedia of Genes and Genomes. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e28\u003c/strong\u003e, 27\u0026ndash;30 (2000).\u003c/li\u003e\n\u003cli\u003eBateman, A. \u003cem\u003eet al.\u003c/em\u003e The Pfam protein families database. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e32\u003c/strong\u003e, D138\u0026ndash;D141 (2004).\u003c/li\u003e\n\u003cli\u003eBairoch, A. The ENZYME data bank. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e22\u003c/strong\u003e, 3626\u0026ndash;3627 (1994).\u003c/li\u003e\n\u003cli\u003eCantarel, B. I. \u003cem\u003eet al.\u003c/em\u003e The Carbohydrate-Active EnZymes database (CAZy): an expert resource for Glycogenomics. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e37\u003c/strong\u003e, D233\u0026ndash;D238 (2009).\u003c/li\u003e\n\u003cli\u003eSeemann, T. Prokka: rapid prokaryotic genome annotation. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e30\u003c/strong\u003e, 2068\u0026ndash;2069 (2014).\u003c/li\u003e\n\u003cli\u003eCantalapiedra, C. P., Hern̗andez-Plaza, A., Letunic, I., Bork, P. \u0026amp; Huerta-Cepas, J. eggNOG-mapper v2: Functional Annotation, Orthology Assignments, and Domain Prediction at the Metagenomic Scale. \u003cem\u003eMol Biol Evol\u003c/em\u003e \u003cstrong\u003e38\u003c/strong\u003e, 5825\u0026ndash;5829 (2021).\u003c/li\u003e\n\u003cli\u003eAlcock, B. P. \u003cem\u003eet al.\u003c/em\u003e CARD 2023: expanded curation, support for machine learning, and resistome prediction at the Comprehensive Antibiotic Resistance Database. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e51\u003c/strong\u003e, D690\u0026ndash;D699 (2023).\u003c/li\u003e\n\u003cli\u003eAltschul, S. F., Gish, W., Miller, W., Myers, E. W. \u0026amp; Lipman, D. J. Basic local alignment search tool. \u003cem\u003eJ Mol Biol\u003c/em\u003e \u003cstrong\u003e215\u003c/strong\u003e, 403\u0026ndash;410 (1990).\u003c/li\u003e\n\u003cli\u003eLiu, B., Zheng, D., Zhou, S., Chen, L. \u0026amp; Yang, J. VFDB 2022: a general classification scheme for bacterial virulence factors. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e50\u003c/strong\u003e, D912\u0026ndash;D917 (2022).\u003c/li\u003e\n\u003cli\u003eChakraborty, A., Ghosh, S., Chowdhary, G., Maulik, U. \u0026amp; Chakrabarti, S. DBETH: A Database of Bacterial Exotoxins for Human. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e40\u003c/strong\u003e, D615\u0026ndash;D620 (2012).\u003c/li\u003e\n\u003cli\u003eO\u0026rsquo;Toole, \u0026Aacute;. \u003cem\u003eet al.\u003c/em\u003e Assignment of epidemiological lineages in an emerging pandemic using the pangolin tool. \u003cem\u003eVirus Evol\u003c/em\u003e \u003cstrong\u003e7\u003c/strong\u003e, (2021).\u003c/li\u003e\n\u003cli\u003eMinh, B. Q. \u003cem\u003eet al.\u003c/em\u003e IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era. \u003cem\u003eMol Biol Evol\u003c/em\u003e \u003cstrong\u003e37\u003c/strong\u003e, 1530\u0026ndash;1534 (2020).\u003c/li\u003e\n\u003cli\u003eLove, M. I., Huber, W. \u0026amp; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. \u003cem\u003eGenome Biol\u003c/em\u003e \u003cstrong\u003e15\u003c/strong\u003e, 1\u0026ndash;21 (2014).\u003c/li\u003e\n\u003cli\u003eDobin, A. \u003cem\u003eet al.\u003c/em\u003e STAR: ultrafast universal RNA-seq aligner. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e29\u003c/strong\u003e, 15\u0026ndash;21 (2013).\u003c/li\u003e\n\u003cli\u003eLiao, Y., Smyth, G. K. \u0026amp; Shi, W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e30\u003c/strong\u003e, 923\u0026ndash;930 (2014).\u003c/li\u003e\n\u003cli\u003eHardwick, S. A., Deveson, I. W. \u0026amp; Mercer, T. R. Reference standards for next-generation sequencing. \u003cem\u003eNature Reviews Genetics 2017 18:8\u003c/em\u003e \u003cstrong\u003e18\u003c/strong\u003e, 473\u0026ndash;484 (2017).\u003c/li\u003e\n\u003cli\u003eWickham, H. ggplot2. \u003cem\u003eWiley Interdiscip Rev Comput Stat\u003c/em\u003e \u003cstrong\u003e3\u003c/strong\u003e, 180\u0026ndash;185 (2011).\u003c/li\u003e\n\u003cli\u003eWaskom, M. L. seaborn: statistical data visualization. doi:10.21105/joss.03021.\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"nature-portfolio","isNatureJournal":true,"hasQc":false,"allowDirectSubmit":false,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"","title":"Nature Portfolio","twitterHandle":"","acdcEnabled":false,"dfaEnabled":false,"editorialSystem":"ejp","reportingPortfolio":"","inReviewEnabled":true,"inReviewRevisionsEnabled":false},"keywords":"Respiratory Microbiome, COVID-19, SARS-CoV-2, Multi-Omics, 16S/ITS Amplicon Sequencing, RNA Sequencing, Metagenomics, Disease Severity","lastPublishedDoi":"10.21203/rs.3.rs-8013295/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-8013295/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eThe outcomes of SARS-CoV-2 infection exhibit significant heterogeneity, and the role of the airway microbiota remains insufficiently understood. In a cohort comprising 561 participants, an integrated multi-omics approach was employed. This approach combined 16S/ITS amplicons (n\u0026thinsp;=\u0026thinsp;542), total RNA metatranscriptomics (n\u0026thinsp;=\u0026thinsp;495), and total DNA metagenomics (n\u0026thinsp;=\u0026thinsp;113), with 101 samples analyzed using all three omics. Metatranscriptomes resolved lineage dynamics across major infection waves and identified microbial coinfections; shotgun metagenomes recovered 190 MAGs, including 13 putative novel species, expanding the catalog of respiratory taxa. Cross‑omics testing showed species‑level biomarkers and resistome burden were more linked to clinical outcomes than genus‑level signals, highlighting the value of high‑resolution profiling. Functional analyses revealed outcome‑based shifts in microbial processes and host metatranscriptomic signatures enriched for interferon programs and tissue remodeling pathways. This work delineates mechanisms associated with COVID‑19 severity, uncovers candidate biomarkers, and provides a multi‑omics framework for integrated pathogen\u0026ndash;microbiome\u0026ndash;host surveillance in future respiratory pandemics.\u003c/p\u003e","manuscriptTitle":"Respiratory Microbiome Dynamics in COVID-19: A Comprehensive Multi-Omics Study","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-11-26 09:09:20","doi":"10.21203/rs.3.rs-8013295/v1","editorialEvents":[],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"communications-biology","isNatureJournal":true,"hasQc":false,"allowDirectSubmit":false,"externalIdentity":"commsbio","sideBox":"Learn more about [Communications Biology](http://www.nature.com/commsbio/)","snPcode":"","submissionUrl":"","title":"Communications Biology","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"ejp","reportingPortfolio":"Communications Series","inReviewEnabled":true,"inReviewRevisionsEnabled":false}}],"origin":"","ownerIdentity":"52c677a9-0101-43cf-a733-4ac8e58ce8a4","owner":[],"postedDate":"November 26th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"under-review","subjectAreas":[{"id":57739773,"name":"Biological sciences/Computational biology and bioinformatics/High-throughput screening"},{"id":57739774,"name":"Health sciences/Diseases/Infectious diseases/Viral infection"},{"id":57739775,"name":"Health sciences/Biomarkers/Prognostic markers"},{"id":57739776,"name":"Biological sciences/Microbiology/Microbial communities/Clinical microbiology"},{"id":57739777,"name":"Biological sciences/Microbiology/Microbial communities/Metagenomics"}],"tags":[],"updatedAt":"2025-11-26T09:09:21+00:00","versionOfRecord":[],"versionCreatedAt":"2025-11-26 09:09:20","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-8013295","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-8013295","identity":"rs-8013295","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}