URT swabs were collected during spring 2020 from 114 adult participants enrolled at two separate sites: 1) 86 participants were enrolled in Nashville, Tennessee: 65 with mild-to-moderate COVID-19 and 21 asymptomatic SARS-CoV-2 uninfected controls, and 2) 28 participants were enrolled in Bronx, New York, who were hospitalized due to severe COVID-19. In all participants, SARS-CoV-2 infection was tested by RT-qPCR.

The 65 participants who had mild-to-moderate COVID-19 were enrolled as part of a randomized clinical trial conducted at Vanderbilt University Medical Center (VUMC) to investigate the effect of nasal irrigations on disease course, as previously described (Kimura et al., 2020). These participants were all seen on an ambulatory basis and none were hospitalized. Concurrently, 21 SARS-CoV-2 RT-qPCR negative adults without COVID-19 or other acute respiratory disease symptoms were recruited from the VUMC community (clinicians, students, faculty, and staff). All participants enrolled at VUMC were given swabs and viral preservation media and performed mid-turbinate swabs as directed; only the enrollment swabs were included in this analysis. Informed consent was obtained from all participants, and this study was approved by the VUMC Institutional Review Board (IRB).

Nasopharyngeal swabs were obtained from patients with COVID-19-like symptoms admitted to the Montefiore Medical Center in the Bronx, New York, by hospital staff at admission. Leftover test swab aliquots in viral preservation media from 28 patients who tested positive for SARS-CoV-2 were sent to VUMC for further laboratory processing. All patients or their surrogate decision-makers signed an informed consent at enrollment, which was approved by the IRB at Albert Einstein College of Medicine.

Further details are available in the Online Supplement . To characterize the URT microbiome, all samples were processed at VUMC with the same method. DNA was extracted with the PowerSoil HTP Kit (Qiagen), the target hypervariable region V4 of the 16S ribosomal RNA (rRNA) gene was amplified with previously published primers to construct libraries (Kozich et al., 2013), and libraries were pooled and then sequenced on an Illumina MiSeq with 2 × 250 bp reads (Shilts M et al., 2020). A ZymoBIOMICS mock community controls (Zymo) and 19 negative controls were processed concurrently with the samples. Sequenced reads were run in R version 4.0.3 through dada2 (Callahan et al., 2016) version 1.18.0 to remove low-quality reads, construct amplicon sequence variants (ASVs), and assign taxonomy against the SILVA reference database (Pruesse et al., 2007). Potential contaminants were removed with the “prevalence” method in decontam (Davis et al., 2018) version 1.10.0. Samples with <1,000 reads were removed (N = 7). For ASVs of interest, if the species was not determined with the dada2 workflow, we used the standard nucleotide basic local alignment search tool (BLAST) to search its sequence against the NCBI 16S rRNA sequences (Bacteria and Archaea) database, excluding uncultured/environmental sample sequences, available at https://blast.ncbi.nlm.nih.gov/Blast.cgi. Sequences were deposited to the Sequence Read Archive at NCBI under BioProjects PRJNA726992 and PRJNA726994.

Bacterial load was assessed using universal 16S rRNA primers as previously described (Barman et al., 2008). Prior to analysis, bacterial copy number was log transformed.

As we had data from two separate sites (Tennessee and New York), and URT sampling methods differed between these two sites (mid-turbinate self-swab versus healthcare provider-performed nasopharyngeal swab, respectively), we looked for microbiome trends that were robust to differences in sites and sampling methods. Therefore, we split the samples from the two sites into five prespecified severity groups (uninfected controls and mild-, moderate-, severe-, and very severe-SARS-CoV-2-infected) to examine if changes in the URT microbiome were consistent over a spectrum of disease severities. As described below, severity groups were chosen to be both clinically relevant and to have an approximately even sample size from each site.

The participants from Tennessee who did not have COVID-19 or other respiratory illness symptoms and were SARS-CoV-2 negative by RT-qPCR were designated the “uninfected” control group (N = 20).

All SARS-CoV-2-infected participants from Tennessee, who had mild-to-moderate illness and were not hospitalized, were given a symptom score questionnaire. At enrollment, participants were asked to rate their symptoms over the last 24 h on an ordinal scale, with 0 indicating the symptom was not present and 7 indicating the symptom was the most severe. Symptoms included were cough, eye redness, nasal congestion, headache, sore throat, sputum production, fatigue, coughing blood, shortness of breath, nausea/vomiting, diarrhea, muscle or joint pain, chills, loss of smell or taste, loss of the ability to think clearly, inability to sleep well, and inability to breathe easily. Summed symptom scores ranged from 0 to 78, with a median (interquartile range) of 30 (15.5–44.5); two participants did not fill out the symptom score questionnaire and so were excluded ( Figure E1 ). Those with a summed symptom score under the median were labeled the “mild” COVID-19 group (N = 27) and those with a summed symptom score at or over the median were labeled the “moderate” COVID-19 group (N = 28).

All SARS-CoV-2-infected participants from New York were hospitalized. Within this group, patients were divided into the “severe” group, who were not admitted to the ICU (N = 15), and the “very severe” group, who were admitted to the ICU (N = 13).

A summary of participant samples included/excluded can be found in Figure E1 . Descriptive statistics were used to characterize the five disease severity groups for the participants who were included in the analysis. A Kruskal–Wallis, Wilcoxon rank-sum, or Pearson chi-squared test, as appropriate, was used to test for differences in variables between the groups.

Microbiome data processing, as described below, was performed in R version 3.5.1 using the wrapper MGSAT (Tovchigrechko, 2020). ASV abundances within each sample were normalized using simple proportions. ASVs with average absolute counts <10 and those with an average relative abundance <0.0005 were aggregated into category “other” which was used when calculating relative abundances but otherwise disregarded. Ninety-five ASVs remained after this filtering. The R package vegan (Oksanen et al., 2014) version 2.5-2 was used to calculate richness and alpha- and beta-diversity at the ASV level. Hill numbers N0, N1, and N2 were used to assess, respectively, richness, the exponential Shannon index, and the inverted Simpson alpha-diversity index (Hill, 1973). Pairwise differences in microbial community composition between all samples were assessed with the Bray–Curtis dissimilarity index, computed on simple proportions, and the PERMANOVA test as implemented in adonis2 (Anderson, 2001) was used to test for significant differences between overall microbial composition and severity groups; age and sex were added to the model. A Tukey HSD post-hoc test was run to examine the significance of pairwise group comparisons. The betadisper function in vegan was used to test for differences in variance between the groups. Variance in Bray–Curtis dissimilarities between each severity group was examined further with the Kruskal–Wallis test and ordinal logistic regression, as described below.

In addition to bacterial load, we preselected the following variables extracted from the bacterial 16S rRNA microbiome data to apply the Kruskal–Wallis test for significant differences between the COVID-19 severity groups: richness, Shannon alpha-diversity, Simpson alpha-diversity, pairwise Bray–Curtis dissimilarities within each severity group, and the 95 ASVs that remained after filtering. The Benjamini–Hochberg correction was applied to adjust for multiple comparisons. Eta-squared (effect size) and its 95% confidence intervals (CI) were found over 1,000 bootstrap replications with the kruskal_effsize function in the R package rstatix version 0.6.0 (Tomczak and Tomczak, 2014; Kassambara, 2020).

To further explore how the microbiome changed as COVID-19 severity increased, we next performed ordinal logistic regression with the lrm function in the R package rms (version 6.2-0) (Harrell, 2021), setting the COVID-19 severity groups as the dependent variable and the microbiome components as the independent variables. Interquartile odds ratios (OR) and their 95% CIs were found with rms::summary.rms and P-values for each independent variable were calculated with the rms::anova function.

As we were interested only in microbiome parameters that consistently increased/decreased as disease severity increased, to minimize overfitting the models, and because some of the microbiome variables may have high collinearity due to how they are calculated (e.g., the richness, Shannon, and Simpson alpha-diversity indices assess a similar phenomenon; the ASV relative abundances could be correlated due to normalization), we applied a strict selection criteria as to which microbiome variables would be added to the ordinal logistic regression models. In addition to bacterial load, microbiome variables were only added to the model if 1) the adjusted Kruskal–Wallis test result P <0.1 and 2) its median either consistently increased or decreased as COVID-19 severity increased.

Age, sex, race, presence of comorbidities, and current smoking were a priori selected to be added to the model as independent variables due to their associations with COVID-19 severity. The race/ethnicity categories were simplified to Black, Hispanic, White, or Other. Comorbidities were reduced to presence of any comorbidity [comorbidities included obesity (defined as a body mass index > 30), diabetes, hypertension, heart disease, or lung disease]. Due to the small number of participants on antibiotics (N = 2) or inhaled steroids (N = 1), we did not add these variables to any models.

Two different main ordinal logistic regression models were built, containing as the independent variables the following: model 1—age, sex, presence of any comorbidities, current smoking, race, bacterial load, bacterial richness, and relative abundance of Corynebacterium_unclassified.ASV0002 and model 2—within-group pairwise Bray–Curtis dissimilarities, which had to be tested separately from all the other variables due to fundamental differences in data structure. The following data were missing for some participants: race/ethnicity (N = 9), comorbidities (N = 3), current smoking (N = 1), and bacterial load (N = 6). We first ran ordinal logistic regression for model 1 (designated model 1A) with only the complete cases (N = 87). Next, we imputed the missing data over 20 imputations with Hmisc ( Harrell and , with contributions from Charles Dupont and many others, 2020) (version 4.4.2) function aregImpute and ran ordinal logistic regression (designated model 1B) using Hmisc::fit.mult.impute with the imputed data so that all cases were included (N = 103). Further details are available in the Online Supplement .

Figures were generated with the R package ggplot2 (Wickham, 2009) version 3.3.3.

A flowchart showing participant samples that were retained for analysis can be found in Figure E1 . Demographic and clinical characteristics of the 103 study participants whose samples were included for analysis can be found in Table 1 . The hospitalized participants (severe or very severe COVID-19) tended to be older, were more likely to be Black or Hispanic, and had a comorbidity rate higher than the non-hospitalized participants. Among hospitalized patients, those who were admitted to the ICU had a higher mortality rate, were more likely to be placed on a ventilator, and had longer hospital stays than those who were not admitted to the ICU. As the incidence of similar symptoms was captured between both cohorts, tables reporting symptom status by COVID-19 severity ( Table E1 ) and age quartile ( Table E2 ) are available in the Online Supplement .

Abundant ASVs were similar among the five severity groups ( Figure 1A ). The most abundant ASVs were Staphylococcus_unclassified.ASV0001 (23.6%), Corynebacterium_unclassified.ASV0002 (13.6%), Corynebacterium_unclassified.ASV0003 (7.0%), Dolosigranulum_pigrum.ASV0006 (4.1%), and Corynebacterium_unclassified.ASV0004 (4.5%). However, when all ASVs were included, bacterial community composition (PERMANOVA P < 0.001, R² = 0.08) and dispersion (betadisper test P = 0.005) differed among the severity groups ( Figure 1B ). Along principal coordinate analysis (PCoA) axis 1 (explaining 7.6% of the variance), the severity groups consistently moved along the same gradient ( Figures 1B and E2 ). A Tukey HSD post-hoc test was run to examine group centroid pairwise comparisons: uninfected control participants compared with the very severe COVID-19 patients were significant [adjusted P = 0.005, difference (95% CI) = 0.08 (0.02–0.14)], while mild compared to very severe (adjusted P = 0.054, difference (95% CI) = 0.06 (−0.001 to 0.12)) and uninfected controls compared with severe [adjusted P = 0.06 difference (95% CI) = 0.06 (−0.002 to 0.12))] approached significance.

The full Kruskal–Wallis results testing for significant differences of the preselected microbiome parameters (i.e., alpha-diversity indices, beta-diversity, bacterial load, and all ASVs passing the abundance cutoff) among the COVID-19 severity groups are available in Table E3 . Only bacterial richness, the within-group Bray–Curtis dissimilarities, and 28 ASVs had an adjusted P <0.1. Of these 30 microbiome parameters passing the adjusted P-value cutoff, only bacterial richness, the within-group Bray–Curtis dissimilarities, and Corynebacterium_unclassified.ASV0002 had medians that consistently increased or decreased as COVID-19 severity increased. These three microbiome parameters, along with bacterial load, were examined further with ordinal logistic regression.

Our results were similar regardless of whether we used only complete cases (model 1A) or the imputed dataset (model 1B); therefore, only the results with the complete cases (model 1A) are presented throughout the rest of the manuscript. The full results are available in the Online Supplement and Tables E4, E5 and Figures E3 – 5 .

Due to their association with COVID-19 severity, age, sex, presence of comorbidities, and current smoking were added to ordinal logistic regression model 1. While increased age [P = 0.65, OR (95% CI) 1.2 (0.55–2.59)], female sex [P = 0.32, OR (95% CI) 1.51 (0.67–3.38)], and presence of comorbidities [P = 0.26, OR (95% CI) 1.9 (0.77–4.69)] were associated with increased disease severity, and current smoking was associated with a reduced risk [P = 0.65, OR (95% CI) 0.61 (0.07–5.07)], only race was significantly associated with COVID-19 severity [P = 0.02, respective ORs (95% CIs) for Black : White, Hispanic : White, and Other : White were 5.31 (1.24–22.73), 10.8 (2.29–50.83), and 3.19 (0.48–21.04)]. Ordinal logistic regression results for each of the tested microbiome parameters are described in detail below.

Six samples failed the bacterial load assay and so were excluded from analysis. Bacterial load was similar in the uninfected controls and mild and moderate COVID-19 participants ( Figure 2 ). Among only those infected with SARS-CoV-2, bacterial load increased as COVID-19 severity increased, although differences between groups were not statistically significant [Kruskal–Wallis adjusted P = 0.23, eta-squared (95% CI) = 0.04 (−0.02 to 0.24)] ( Figure 2 and Table E3 ). However, when we performed ordinal logistic regression, bacterial load was significantly associated with disease severity [P = 0.04, OR (95% CI) = 2.09 (1.04–4.21)].

Bacterial richness and alpha-diversity generally increased as disease severity increased, although alpha-diversity dropped in patients with very severe COVID-19 ( Figure 3 ). This change in richness was near significant with the Kruskal–Wallis test [adjusted P = 0.050, eta-squared (95% CI) = 0.09 (0.01–0.27)], while neither the Shannon nor Simpson indices were significant ( Table E3 ). Richness was further examined with ordinal logistic regression, and it was near significantly associated with COVID-19 severity [P = 0.08, OR (95% CI) = 1.51 (0.96–2.37)].

After calculating pairwise Bray–Curtis dissimilarities between all samples and examining only the within-group pairs, we observed that dissimilarity between samples increased as COVID-19 severity increased ( Figure 3 ). Within-group dissimilarities were significantly different with a Kruskal–Wallis test [adjusted P < 0.001, eta-squared (95% CI) = 0.16 (0.13–0.21)] and were significantly associated with COVID-19 severity in our ordinal logistic regression model [P < 0.001, OR (95% CI) = 2.69 (2.28–3.17)]. When we looked at all Bray–Curtis dissimilarities, within-group dissimilarities in the very severe COVID-19 group were among the highest, similar to the dissimilarities between different COVID-19 severity groups ( Figure E6 ).

Corynebacterium_unclassified.ASV0002 relative abundance was found to be significantly different between severity groups with Kruskal–Wallis test [adjusted P = 0.04, eta-squared (95% CI) = 0.10 (0.03–0.26)], decreasing as disease severity increased ( Figure 4 ). However, in our ordinal logistic regression model, Corynebacterium_unclassified.ASV0002 relative abundance was not significantly associated with disease severity [P = 0.11 and OR (95% CI) = 0.69 (0.42–1.09)]. The full-length sequence matched with 100% identity to both Corynebacterium accolens and Corynebacterium macginleyi.