This was a prospective observational study on adult (age ≥ 18 years old) hospitalized patients who tested positive for SARS-CoV-2 on reverse transcription polymerase chain reaction (RT-PCR). Patients were included if they had at least one nasopharyngeal swab and throat swab (NPSTS) sample taken during hospitalization within 3 weeks of hospital admission after informed consent. Patients who were previously vaccinated against SARS-CoV-2, received antibiotics 3 months prior to hospitalization, or had missing data on clinical severity or antimicrobial therapy were excluded. Samples that were inadequate for DNA extraction and 16S rRNA gene amplicon sequencing for microbiota profiling were also excluded. Blood samples for cytokine profiling were taken as early as possible after hospital admission. Mechanically ventilated intensive care unit (ICU) patients without COVID-19 and healthy volunteers working in the same hospital environment were recruited for controls. This study was approved by The Joint Chinese University of Hong Kong – New Territories East Cluster Clinical Research Ethics Committee (2020.076).

COVID-19 severity was classified as asymptomatic, mild/moderate, or severe/critical based on the highest severity level at hospital discharge as previously described (28). Medical records including clinical notes, imaging, laboratory results, and observation charts were manually reviewed to determine the severity of COVID-19 according to the following criteria: Asymptomatic patients had no symptoms despite SARS-CoV-2 infection. Mild/moderate group included patients who had symptoms of fever, cough, myalgia, sore throat, and rigors related to SARS-CoV-2 infection but did not require oxygen therapy. Severe/critically ill patients with COVID-19 included those who had dyspnea, respiratory rate ≥ 30, or required oxygen therapy or mechanical ventilation for SARS-CoV-2 infection due to respiratory failure.

Serial NPSTSs were collected for SARS-CoV-2 viral load quantification and 16S rRNA sequencing analysis. During the early COVID-19 pandemic, all patients with confirmed COVID-19 were hospitalized as part of the Hong Kong public health strategy. Patients were only discharged after testing negative for SARS-CoV-2 on RT-PCR. Patients who were asymptomatic or had mild disease were, sometimes, hospitalized for longer than 2 weeks despite clinical recovery. Therefore, serial NPSTS samples could be collected from all severity groups during hospitalization. The time points used in this study were the first sample within the first week of hospitalization and the second sample between the second and third weeks of hospitalization. Samples were stored in −80°C for 0.1–2.5 years until completion of recruitment for further analysis.

Total RNA was extracted from mixed NPSTSs using the QIAamp Viral RNA Mini Kit (QIAGEN, Hilden, Germany). Primer–probe set targeting the N gene (2019-nCoV_N1-F: 5′-GAC CCC AAA ATC AGC GAA AT-3′; 2019-nCoV_N1-R: 5′-TCT GGT TAC TGC CAG TTG AAT CTG-3′; and 2019-nCoV_N1-P: 5′-FAM-ACC CCG CAT TAC GTT TGG TGG ACC-BHQ1-3′) was used to detect SARS-CoV-2 RNA by real-time RT-PCR as previously described (Lui et al., 2020). The detection limit of real-time RT-PCR was 694 copies/ml, and samples were considered negative if Ct values exceeded 39.9 cycles.

Total DNA was extracted from mixed NPSTS using the QIAamp DNA Mini Kit (QIAGEN, Hilden, Germany) to characterize respiratory microbiota using 16S rRNA gene amplicon sequencing. The molecular process, including the DNA extraction, 16S PCR amplification, and library preparation, were performed at separate locations to avoid contamination. For quality control, negative controls (blank DNA extraction and PCR controls), positive controls (ZymoBIOMICS Microbial Community DNA Standard, catalog no. D6305), and technical replicates (randomly selected DNA samples) were also included. In brief, the 16S rRNA gene hypervariable V3-V4 region (~450 bp) was targeted (341F: 5′-CCT ACG GGN GGC WGC AG-3′; 806R: 5′-GGA CTA CNV GGG TWT CTA AT-3′), with barcodes indexed to each amplicon set for multiplexing sequencing on an Illumina MiSeq for PE300 reads (Chen et al., 2019). QIIME2 (v2022.2) with the latest SILVA ribosomal RNA database (v138 SSU Ref NR 99 dataset) was used to classify amplicon sequence variants (ASVs), with operational taxonomic table showing the proportion of bacterial reads per sample at different taxonomic levels after removing reads assigned to archaea, mitochondria, or chloroplasts.

Data distribution was assessed using Shapiro–Wilk test, and descriptive statistics such as mean and standard error as well as median and interquartile range (IQR) were used to summarize data. A phylogenetic tree was generated by inserting the representative reads into the SILVA 128 reference tree using the SATe-enabled phylogenetic placement method. Alpha diversity was assessed using Shannon, Simpson, richness, and evenness, and Wilcoxon rank sum test was used for assessing the pairwise difference between the defined groups. Beta diversity was assessed using unweighted and weighted GuniFrac and Bray–Curtis distance. Principal coordinate analysis was used to assess beta diversity between different groups using permutational multivariate analysis of variance (PERMANOVA) with 9,999 permutations using the adonis2 in the Vegan R package (v2.6-4). In the effect size analysis using a single multivariable model, antibiotic-controlled association between metadata variables, including intubation, ICU, severity, peak CRP, hospitalized time, antivirus, peak viral load, Charlson comorbidity index, age, and gender, was tested by adding antibiotics into the model formula. An exploratory analysis on discriminative bacterial taxa between patients who required ICU admission and those that did not was estimated using linear discriminant analysis (LDA) effect size (LEfSe) with the default setting, with further comparisons of the relative abundances using nonparametric Mann–Whitney–Wilcoxon rank sum test and Tukey’s honest significant difference post-hoc test (Segata et al., 2011). Logistic regression and receiver operating characteristic (ROC) curve with the calculation of area under the ROC curve (AUC) were used to evaluate the potential markers identified for prediction of need for ICU admission. Delong’s test was used to assess the differences in AUCs. All other data visualization was performed using the ggplot package in R. A two-sided p-value < 0.05 or a false discovery rate–adjusted p-value (p_adj or q) < 0.1 was used as the threshold for statistical significance.

A 3 ml of EDTA blood sample was taken from recruited patients with COVID-19 and immediately cooled and transported to the laboratory for processing. Plasma was separated by centrifugation (2,000g for 10 min) at 4°C and stored in 300 µl of aliquots at − 70°C until analysis. Milliplex human cytokine multiplex assay using the Bio-plex 200 System (Bio-Rad Laboratories, Inc. CA, USA) was used to determine levels of 32 cytokines including sCD40L, EGF, Eotaxin, FGF-2, Flt-3L, Fractalkine, GRO-α, IFN-α2, IFN-γ, IL-1β, IL-1RA, IL-3, IL-5, IL-6, IL-7, IL-8, IL-10, IL-12 p40, IL-12 p70, IL-13, IL-15, IL-18, IP-10, MCP-1, MCP-3, MDC, MIG, MIP-1β, TGF-α, TNF-α, TNF-β, and VEGF (Ling et al., 2021). Associations of cytokine factors with clinical variants and bacterial genera were analyzed by Fit a Negative Binomial Generalized Linear Model using the glm.nb in the MASS R package and a Spearman’s rank-order correlation test using the cor.test in the Stats R package, respectively.

A total of 314 NPSTS samples from 197 subjects that generated high-quality 16S sequence reads were analyzed ( Figure S1 ). COVID-19 group consisted of 171 adult hospitalized patients (asymptomatic = 14, mild/moderate = 106, and severe/critical = 51) ( Table 1 ). Meanwhile, 11 mechanically ventilated ICU adult hospitalized patients without COVID-19 and 15 adult healthy volunteers who worked in the same hospital as healthcare workers or departmental staff were included as controls. The baseline characteristics of the recruited subjects are shown in Table 1 and Table S1 . Overall, patients with severe/critical COVID-19 were older and were more likely to have comorbidities. Peak median viral load was similar across different COVID-19 severity ( Table 1 , p = 0.086). Antibiotic use was highest in patients with severe/critical COVID-19 and non–COVID-19 ICU patients. Use of specific COVID-19 treatments varied across the spectrum of COVID-19 severity.

We selected the earliest collected NPSTS samples from each subject (hospitalized COVID-19 = 171, non–COVID-19 ICU = 11, and healthy = 15) to profile the upper airway microbial communities. The median time interval between time of hospitalization and sampling time of first sample was 3 days. A total of 2,649,726 high-quality 16S rRNA V3-V4 reads were generated, ranging between 1,026 and 118,093 reads per sample (13,450 ± 11,799). As shown in Figure 1 and Figure S2 , Firmicutes was the most abundant microbial phyla in the surveyed samples (mean relative abundance ± SD of 34.70 ± 1.20%), followed by Bacteroidota (30.06 ± 1.30%), Proteobacteria (17.86 ± 1.32%), and 11 other phyla. Proteobacteria were lower in both hospitalized COVID-19 (17.70 ± 1.44%, Mann–Whitney U-test, p < 0.01) and non–COVID-19 ICU (6.52 ± 2.10%, p <0.001) patients compared with healthy volunteers (27.95 ± 3.62%). In contrast, Bacteroidetes was higher in both in both hospitalized COVID-19 (30.47 ± 1.30%, p < 0.001) and non–COVID-19 ICU (44.09 ± 8.78%, p < 0.01) patients compared with healthy volunteers (30.06 ± 1.30%).

Alpha diversity at the ASV level as measured by Shannon, Simpson, and evenness was progressively reduced as severity of COVID-19 increased ( Figure 2A , Figure S3A ). Similarly, distinct clustering of the microbial communities from patients with COVID-19, particularly those with severe/critical severity, from healthy controls was observed in weighted GUniFrac, unweighted GUniFrac, and Bray–Curtis distance at the ASV level (p < 0.001). Antibiotic treatment (yes vs. no) was consistently the most significant covariate that had the largest effect on the overall structure of the upper airway microbiota in all beta diversity metrics (unweighted GUniFrac: R^{2 =} 0.0442, p < 0.001; weighted GUniFrac: R^{2 =} 0.0307, p < 0.001; and Bray–Curtis: R^{2 =} 0.0216, p < 0.001) ( Figure 2A , Figure S3B , Table S2 ). Similarly, the microbial community was also significantly affected by intubation (yes vs. no) and ICU admission (yes vs. no). The effect of severity (severe/critical vs. non-severe/critical), peak CRP (≤ 7 vs. > 7), hospitalized time (≤ 1 vs. ≥ 2 weeks), and antivirus treatment (yes vs. no) on structure of microbial community was not detected in all measures of beta diversity ( Figure 2A , Figure S3C ). Notably, 86.7% (OR = 40.3; 95% CI, 8.3–392.5; p < 0.001) and 62.5% (OR = 14.5; 95% CI, 5.5–40.6; p < 0.001) of hospitalized patients with COVID-19 who were intubated and admitted to the ICU, respectively, were treated with antibiotics prior to NPSTS sample collection in this study.

Because antimicrobial use was the most significant factor that affected upper airway microbiota dysbiosis, we further assessed the difference in microbiota without the effect of antimicrobials. Removal of NPSTS samples taken after antimicrobial use resulted in 157 samples from 137 antibiotic-naïve hospitalized patients with COVID-19 (14 asymptomatic, 98 mild/moderate, and 25 severe/critical cases), five non–COVID-19 ICU patients, and 15 healthy volunteers. The median time interval between time of hospitalization and sampling time of first sample was 3 days. All measures of alpha diversity were no longer significantly different across COVID-19 severity groups when confounding by antimicrobial was removed ( Figure 2B , Figure S4A ). After exclusion of the samples taken after antimicrobials, microbiota from patients with severe/critical COVID-19 still showed dispersive distribution by beta diversity analysis when compared with other groups, but the differences were reduced ( Figure 2B , Figure S4B ). After removal of the samples taken after antimicrobial use, the confounding effect of intubation and ICU admission on upper respiratory microbiota in COVID-19 was no longer consistently found ( Figure S4C ).

Because upper airway microbiota is significantly confounded by use of antimicrobial use, longitudinal change in microbiota during hospitalization was assessed after exclusion of samples taken after antimicrobial therapy (102 samples). Assessment of these samples showed that upper airway microbiota alpha and beta diversity did not change overtime during 2 weeks of hospitalization in antibiotic-naïve hospitalized patients with COVID-19 and healthy individuals ( Figure 3 , Figure S5 ).

Differentially abundant bacterial genera associated with demographic and clinical variants were characterized using LDA by LEfSe in the surveyed NPSTS samples from patients with COVID-19 in an exploratory analysis (n = 171) (LDA score > 2, p < 0.05) ( Figure 4 , Table S3 ). Interestingly, 15 and 10 bacterial genera showed consistent increase or decrease in the relative abundance in patients with antimicrobial use, intubation, and/or ICU admission. Among these, the enrichment of four bacterial genera (Enterococcus, Limosilactobacillus, Sneathia, and Pseudomonas) and the depression of eight bacterial genera (Alloprevotella, Prevotella, Butyrivibrio, Hespellia, Lachnoanaerobaculum, Oribacterium, Solobacterium, and Centipeda) were also significantly associated with the severity of patients with COVID-19. Because 63% (32 of 51) of patients with severe/critical COVID-19 in this study were admitted to the ICU, we further selected early NPSTS samples prior to antimicrobial use within the first week of hospitalization (n = 116) to identify bacterial markers that may be used to predict need for ICU admission ( Figure 5 ). A total 13 discriminative bacterial genera were observed using LDA, with Enterobacter, Mageeibacillus, Fannyhessea, Scardovia, Howardella, and Bulleidia being higher but with Alloprevotella, Campylobacter, Leptotrichia, Centipeda, Hespellia, Catonella, and Acinetobacter being lower in patients who required ICU admission. The differential abundances of these 13 bacterial genera in the upper airway tract within the first week of hospitalization predicted the need for ICU admission with a combined AUC of 0.95 (95% CI of 0.91–0.99). Comparatively, prediction based on clinical demographics and comorbidity (age, gender, and Charlson comorbidity index) only achieved a combined AUC of 0.71 (95% CI of 0.57–0.85). Combination of clinical and microbiota (AUC of 0.96 (95% CI 0.93–1.00) did not improve on the predictive performance compared with the use of microbiota alone ( Figure S6 , p = 0.486).

Plasma was collected from 90 patients with COVID-19 to measure the levels of 32 cytokines. Median (IQR) time to cytokine profiling was 4 (2–8) days after hospital admission. Using the Negative Binomial Generalized Linear Model, changes in numerous cytokine factors were significantly associated with clinical factors ( Figure 6 ). For example, increased IL-5, IL-6, IL-10, and MIG but deceased MIP-1β were observed in patients with severe/critical COVID-19 or those admitted to ICU, with satisfactory AUCs for predictive performance ( Table S4 ). Those patients with higher IL-6 also had a higher chance of intubation and antibiotics and/or antivirus use. Older patients or those with higher Charlson comorbidity had higher levels of IL-6, IL-10, IP-10, MCP-1, and MIG, but lower levels of Fractalkine, IFN-χ, IL-12 p70, IL-13, and TGF-α. To understand the potential of upper airway microbiota dysbiosis on plasma cytokine levels, Spearman correlations between bacterial genera and cytokines were explored ( Figure S7 , Table S5 ). Next, we focused on the 13 bacterial genera that were associated with a need for ICU admission and found that only Acinetobacter, Hespellia, and Campylobacter were associated with MIG, IL-18, Fractalkine, and IL-1β after removal of samples taken after antimicrobial therapy ( Figure 7 ).