Nasopharyngeal Microbiome Community Composition and Structure Is Associated with Severity of COVID-19 Disease and Breathing Treatment

Abstract

Viral infections are known to modulate the upper respiratory tract microbiome, but few studies have addressed differences in the nasopharyngeal microbiome following SARS-CoV-2 infection. Using nasopharyngeal swab medical waste samples from 79 confirmed SARS-CoV-2 positive and 20 SARS-CoV-2 negative patients, we assessed microbiome composition with metagenomic sequencing. COVID-19 status and breathing assistive device use was associated with differences in beta diversity, principal component analyses, community composition and abundance of several species. Serratia more frequently appeared in COVID-19 patient samples compared to negative patient samples, and Serratia, Streptococcus, Enterobacter, Veillonella, Prevotella, and Rothia appeared more frequently in samples of those who used breathing assistive devices. Smoking and age were associated with differences in alpha diversity. Cross-sectional differences in the microbiome were apparent with SARS-CoV-2 infection, but longitudinal studies are needed to understand the dynamics of viral and breathing treatment modulation of microbes.

1. Introduction

Although much research has been rapidly produced regarding the severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) pandemic, only a handful of studies have examined associations of viral infection with the nasopharyngeal (NP) microbiome [1,2,3,4,5]. It is hypothesized that SARS-CoV-2 infection may compromise immunity and increase bacterial attachment, growth, and possibly lead to co-infections or “super infections” [5,6]. Bacterial infection can further exacerbate inflammatory responses causing tissue damage and further dissemination of both the virus and bacteria [6].

The microbiota is known to promote viral infections through several mechanisms including enhancement of virion stability [7,8], contributing to viral replication [9], and suppressing antiviral immunity [10]. Alternatively, commensal microbiota can suppress viral infections directly by binding to viruses and thereby suppressing replication [11], destabilizing morphology [12], and/or inhibiting infectivity [13], or via other indirect mechanisms [14]. A recent review covers the interactions of microbiota and viruses extensively [15]. Importantly, interactions between the commensal microbiota and viral pathogens have been linked to clinical outcomes [8].

Studies of the NP microbiome in SARS-CoV-2 positive and negative samples have found a range of observations; from no differences in those with mild COVID-19 to species-specific alterations that differed in reports from different cities and patient populations [2,3,4,16,17]. In one study of co-infections in the respiratory tract (n = 100 positive) that found limited viral genetic variation among patients, the authors suggested that the presence of particular species (e.g., Clostridium botulinum, Bacillus cereus and Halomonas spp.) was potentially a more important factor on clinical outcomes [5]. While it is difficult to draw strong conclusions from this study due to a lack of control group, the observation underscores the need to understand whether SARS-CoV-2 infection enhances the respiratory tract for opportunistic coinfections or if an impoverished microbiome increases susceptibility to severe disease.

The present study uses shallow shotgun sequencing of a cross-sectional set of medical waste samples from Louisiana, USA that were confirmed PCR positive for SARS-CoV-2 across a range of severity (asymptomatic to hospitalization and death). We hypothesize that COVID-19 status will be associated with a shift in the NP microbiome, especially in hospitalized patients who have been previously observed to undergo rapid and dramatic shifts in the gut microbiome upon admission to the hospital [18].

2. Materials and Methods

2.1. Samples

A total of 99 medical waste samples from patients tested for SARS-CoV-2 infection at Ochsner Health (79 SARS-CoV-2+ and 20 SARS-CoV-2−; detected by RT-PCR testing on the Abbott m2000 device, Abbott Laboratories) were included in the analysis. The samples were nasopharyngeal swabs in viral transport media (VTM) collected under IRB-approved protocol 2019.334 along with metadata in 3–8 April 2020 in Louisiana, USA. The utility of this type of medical waste sample for microbiome studies has been previously established for Influenza [19] and SARS-CoV-2 [2,20]. Demographics associated with the samples are presented in

Table 1. Summary characteristics of SARS-CoV-2 positive and negative patients whose nasopharyngeal microbiome was analyzed.

	SARS-CoV-2 Positive		SARS-CoV-2 Negative
	n	%	n	%
Total	79	100%	20	100
Inpatient	14	18%	1	5%
Any Breathing Assist Device	12	15%	0	0%
Disposition (dead)	5	6%	2	10%
Severity 1
low	2	3%	N/A
mild	43	54%
moderate	20	25%
high	2	3%
severe	12	15%
Proton Pump Inhibitor	14	18%	5	25%
Smoker (yes)	8	10%	3	15%
Inhaler	21	27%	9	45%
Antibiotic	15	19%	2	10%
Average BMI (±SD)	32.6	±8.2	28.8	±5.9
Average Age (±SD)	51.6	±17.0	50.3	±14.0
Race
Black	44	56%	11	55%
White	35	44%	9	45%
Sex
Male	27	34%	7	35%
Female	52	66%	13	65%

¹ Low, no or very few symptoms; mild, symptoms that were treated at home; moderate, symptoms were severe enough to require an emergency visit but not admission; high, symptoms required admission; severe, symptoms required breathing assistive devices in the hospital. BMI, body mass index.

and are representative of the demographic breakdown of testing at the time.

2.2. Clinical Attributes

Patient records were manually searched for information on clinical attributes. Chart data were collected only for comorbidities and concomitant medications that were active up to the date of testing. For example, if someone was tested on 8 April after being admitted 4 April, they would be classified as admitted, any breathing treatment prior to 8 April would be recorded and any medications would be included. Disposition (alive/dead) was recorded 3 months after the test date. Patients were classified as having breathing treatment if they had treatment with oxygen mask, non-rebreather mask, nasal cannula, BiPAP, CPAP, or ventilator. Smoker status, BMI, age, and race were all variables recorded at the last visit prior to testing. If any of these variables were unknown, the sample was excluded from that specific analysis. Severity of COVID-19 was recorded based on initial reports from the test (e.g., “patient reports coughing, shortness of breath”) and the greatest severity recorded 3 months after testing (e.g., patient initially went to the emergency department, treated symptoms at home and was later admitted would be classified as “severe”).

A total of 13 attributes were analyzed. Five categorial attributes (COVID-19 status, inpatient/outpatient, breathing assistance, disposition, and severity) were chosen that were hypothesized to have a direct association between COVID-19+ disease and microbiome; four categorial attributes (proton-pump-inhibitor, smoking status, antibiotic usage, and inhaler usage) and two continuous attributes (BMI, age) were chosen that were hypothesized to have an indirect association between COVID-19+ disease and microbiome; and two categorical attributes that were hypothesized to have no association between COVID-19+ disease and microbiome (race, gender). Data for each patient is presented in Supplemental Tables S1 and S2.

2.3. Metagenomic Sequencing and Classification

NP swabs in VTM were shipped overnight to the BioInfoExperts, LLC, BSL-2 laboratory. DNA was extracted using the Zymo Quick-DNA Kit (Zymo, Irvine, CA USA) and eluted in 50 μL DNase/RNase-Free water as per manufacturer’s instructions and shipped to CosmosID. DNA samples were quantified via Qubit4 (Thermofisher, Waltham, MA USA) and normalized for library preparation via Nextera XT DNA Library Preparation Kit (Illumina, San Diego, CA USA) following manufacturer’s protocol. DNA libraries were then quantified via Qubit4 (Thermofisher, Waltham, MA USA) and pooled in equimolar ratio for sequencing via Illumina NovaSeq 6000 (Illumina, San Diego, CA USA).

Classification was performed using CosmosID’s analytical pipeline (Cosmos ID, Rockville, MD, USA). Specifically, k-mers for each read were compared to their curated microbial genomics database (Genbook) which contains nearly 160,000 phylogenetically organized genomes and gene sequences. For each sample, an Abundance Score (functionally equivalent to an operational taxonomic unit (OTU) count) was assigned for all OTUs identified at a strain level. OTUs were agglomerated to a genus level in R package phyloseq.

2.4. Alpha Diversity

Three alpha diversity metrics (Shannon, Inverse Simpson, and number of observed OTUs) were calculated using the R package phyloseq. The Kruskal–Wallis Rank Sum (KWRS) test was used to determine whether alpha diversity per sample was different between patients for categorical attributes. The Spearman rank correlation test and a generalized linear regression model were used to determine differences in diversity for continuous attributes. All significance testing was performed in R. The procedure was adapted from that described in The Riffomonas Reproducible Research Tutorial Series (code www.riffomonas.org/minimalR, accessed on 2 January 2021).

2.5. Beta Diversity

Beta diversity testing used 11 categorical attributes. A distance matrix was prepared using centered log-ratio transformed filtered relative abundances. The permutational multivariate analysis of variance (PERMANOVA) which partitions the sums of squares for within- and between-cluster distances was performed on the transformed values with 999 permutations to infer the null distribution. Differences in dispersion, which can confound the PERMANOVA analysis, were also assessed. All beta diversity testing was performed using R package vegan. The procedure was adapted from that described in the Introduction to the Statistical Analysis of Microbiome Data in R (code at https://github.com/Nick243, accessed on 2 January 2021).

2.6. Community Composition

We examined the frequency with which particular genera were present or absent in patients for a subset of attributes: two that were hypothesized to be related to COVID-19 disease: (COVID-19 status, breathing assistance) and one that was hypothesized to be indirectly related (antibiotic usage). For each of the attributes, the total number of patients as well as the percentage of patients for whom each genus was present or absent was visualized. The Fisher’s Exact test was used to determine significantly different ratios using an alpha value of 0.05 and a one-tailed hypothesis.

2.7. Differential Abundance

Three attributes (COVID-19 status, breathing assistance, and antibiotic use) were investigated for differences in relative abundance. A subset (n = 11) of the 21 genera was tested that was either a known opportunistic pathogen, had previously been identified as associated with COVID-19 disease, and/or were identified in the community composition analysis as being present at a significantly different frequency between groups (subset: Actimomyces, Enterobacter, Finegoldia, Neisseria, Peptoniphilus, Prevotella, Rothia, Serratia, Staphylococcus, Streptococcus, and Veillonella). For statistical testing, genera were included only if present in at least 3 or more patients in each group for the three attributes. Significance was assessed using the Kruskal–Wallis Rank Sum (KWRS) test in R. The procedure was adapted from that described by the Riffomonas “Minimal R” protocol (http://www.riffomonas.org/, accessed on 2 January 2021).

3. Results

3.1. Classification

A total of 202 unique genera was observed across all samples. After removing the index organism (Imtechella halorens), the relative abundance was calculated for each genus. Genera that were observed in <5% of samples and/or were found at a combined relative abundance of <0.5% were excluded, resulting in a total of 27 unique genera (Acinetobacter, Actinomyces, Aggregatibacter, Bradyrhizobium, Burkholderia, Burkholderiaceae [unclassified], Corynebacterium, Cutibacterium, Delftia, Dolosigranulum, Enterobacter, Finegoldia, Kocuria, Micrococcus, Mycobacterium, Neisseria, Peptoniphilus, Pilimelia, Porphyromonas, Prevotella, Propionibacteriaceae [unclassified], Propionibacterium, Rothia, Serratia, Staphylococcus, Streptococcus, and Veillonella). The abundance of these 27 genera was used for alpha and beta diversity testing.

3.2. Alpha Diversity

For 11 categorical attributes, diversity as measured using all three metrics (Shannon, Inverse Simpson, and Observed OTUs) was similar among patients (p > 0.05, Supplemental Table S3). One attribute (smoking) was significantly associated with the Inverse Simpson measure of diversity (p = 0.047), although not with the Shannon metric. Age was also significantly associated with Shannon and Inverse Simpson measures of diversity (p = 0.016 and p = 0.007, respectively). For age, both metrics were significant under a generalized linear model regression, although with low R² values (Inverse Simpson: p = 0.002, R² = 0.083; Shannon: p = 0.010, R² = 0.057).

3.3. Beta Diversity

A significant difference in beta diversity among patients was observed for four attributes hypothesized to have a direct impact on the microbiome (

Table 2. Beta diversity in the microbiome associated with patient attribute.

Characteristic	Patient Subgroup	Adonis p-Value	Dispersion p-Value
Breathing assistance	COVID-19+	0.003	0.104
Disposition	COVID-19+	0.010	0.637
Inpatient	COVID-19+	0.004	0.283
Severity	COVID-19+	0.109
Antibiotic	All	0.213
COVID-19 status	All	0.012	0.834
Inhaler	All	0.072
Proton pump inhibitor	All	0.953
Smoker	All	0.070
Race	All	0.001	0.001
Sex	All	0.138

Bold indicates significant difference.

): breathing assistance (p = 0.001); COVID-19 status (p = 0.012); disposition (p = 0.014); and inpatient status (p = 0.040). A significant difference in beta diversity among patients was found for a fifth attribute (race; p = 0.001), which was expected to have a limited impact on microbiome. For the five attributes with significant differences in beta diversity, the dispersion was plotted, and for the four attributes tested for COVID-19+ patients, no significant differences were found (p > 0.05) suggesting a true underlying difference. However, for race, the dispersion between groups was significant (p = 0.001), suggesting that the difference in beta diversity was driven by differential dispersion rather than biological difference.

To investigate the four attributes showing a significant difference in beta diversity and hypothesized to play a direct role on the microbiome, we prepared plots for the first two principal components (Figure 1). As expected, none of the attributes showed complete separation into groups, but observable difference was observed among some patients who required breathing assistance (Figure 1a) and hospital admission (Figure 1b), but fewer for COVID-19 status (Figure 1c).

3.4. Community Composition

Overall, a similar frequency of occurrence of a genus was observed in patients who tested positive as those who tested negative for COVID-19 (Figure 2a; Supplemental Table S4). Serratia was detected significantly more frequently in positive vs. negative patients (p = 0.02). One genus (Kocuria) was detected only in COVID-19+ patients (n = 12); however, this association was detectable just above the alpha value for significance (p = 0.055). Two genera (Finegoldia, Peptoniphilus) were entirely absent in COVID-19 patients who required breathing assistance (n = 12), while present in a subset of COVID-19+ patients who did not require breathing assistance (Finegoldia: n = 24; Peptoniphilus: n = 19), e.g., Finegoldia: p = 0.009; Peptoniphilus: p = 0.027; Figure 2b; Supplemental Table S4. Two genera were absent in patients who did not receive antibiotics (Kocuria and Porphyromonas) but present in a subset of patients on antibiotics (Figure 2c).

3.5. Differential Abundance

For COVID-19 status, the abundance of two genera (Serratia and Veillonella) was significantly different in COVID-19+ vs. COVID-19− patients (p = 0.028 and p = 0.038, respectively; Supplemental Table S5). For breathing assistance, the abundances of six genera (Enterobacter, Prevotella, Veillonella, Streptococcus, Rothia, Serratia) were significantly different in patients who required breathing assistance (p = 0.007–0.024; Supplemental Table S5). For antibiotic usage, the abundance of two genera (Finegoldia, Serratia) were significantly different (p = 0.009 and p = 0.012, respectively).

Boxplots of the relative abundance for each patient showed a clear outlier for Veillonella for COVID-19 status, two clear outliers for Streptococcus and breathing assistance, and two clear outliers for Veillonella and antibiotic usage. These outliers were removed, and the statistical test was performed again, and the p-values remained significant.

We plotted relative abundances as boxplots for each patient for the seven genera showing significantly different abundances (Figure 3). For COVID-19 status, the median relative abundance of Serratia for COVID-19+ patients was clearly well above those for COVID-19− patients, while individual relative abundances of Veillonella were similar between groups (Figure 3a). For breathing assistance, the median relative abundance of all six genera was higher in patients who required breathing assistance compared to those who did not. However, for three of the genera (Serratia, Streptococcus, and Enterobacter), the median differences could be related to the smaller sample size than difference in distribution. On the other hand, for the three genera (Rothia, Prevotella, and Veillonella), higher abundance in patients requiring breathing assistance appeared to be consistent for the majority of patients, suggestive of an underlying biological difference. For antibiotic use, while the median relative abundance of Serratia was higher in individuals taking antibiotics, individual relative abundance of Veillonella was similar between groups. It should be noted that the relative abundance of Finegoldia was clearly higher in the antibiotic usage group.

Results are available via an interactive interface at https://public.tableau.com/app/profile/rebecca.rose5228/viz/BioInfoExpertsOchsner_COVID-19_Microbiome_Study/TaxonomyHierarchy, accessed on 7 February 2021.

4. Discussion

Comparing the nasopharyngeal microbiomes of 79 SARS-CoV-2 positive samples, representing a range of COVID-19 severity, and 20 negative samples, revealed differences in diversity, frequency, and abundance of several species with respect to COVID-19 status. While alpha diversity did not differ according to COVID-19 status, it was different with respect to age and smoking status. It should be noted that both age and smoking status have previously been shown to influence the gut, skin, and upper respiratory tract microbiome [21,22].

Beta diversity of the microbiome was significantly different for those who had COVID-19, used breathing assistive devices, had an inpatient hospital stay, and, sadly, for those who died. These factors often showed high overlap; breathing assistive devices were used only for inpatients suffering COVID-19, although not all who had COVID-19 needed breathing assistance or required an inpatient hospital stay. One patient who died was COVID-19 negative, but that the death was linked to other categories (see Supplemental Table S1 for specific patient details).

It is recognized that race is a social construct and not a biological factor, but race is shown repeatedly to be associated with health outcomes. Surprisingly, we observed racial differences in beta-diversity of the NP space, but with the caveat that the dispersion was significantly different between racial groups, suggesting that the differences in diversity can be driven by the wide variance rather than a biological factor. Research on racial differences with respect to the microbiome is a topic of high interest currently. Investigators have reported that oral microbial richness and composition were different among African- and European-Americans in a study of more than 1500 individuals [23]. Importantly, this difference was noted for self-reported race as well as percent of those with genetically determined African ancestry. The Healthy Life in an Urban Setting (HELIUS) study included 2084 participants and demonstrated ethnic differences in the diversity and abundance of species comprising the gut microbiome, despite the fact that all subjects were residents of Amsterdam, Netherlands [24]. A larger sample size will be required to determine whether actual racial differences exist, but this preliminary finding emphasizes that race should be considered in microbiome studies.

COVID-19 positivity was associated with increased representation of Serratia. It is important to note that Serratia marcescens is a recognized causative agent of human diseases, including pneumonia. Other investigators have noted the paradoxical risk of S. marcescens infection in intensive care units during the current pandemic [25]. Despite taking increased sanitation measures and great care in PPE donning and doffing, highly resistant S. marcescens infections occurred in five ICU patients in Paris, France, likely because of decreased hand hygiene related, ironically, to increased sterile glove use. [25] Kocuria spp. were detected only in the COVID-19 positive patients and, although not at a significant level, it is curious that none of the COVID-19 negative samples contained this genus since it is ubiquitous in its occurrence on human skin and mucus membranes. Kocuria species have been identified as a causative agent in hospital acquired infections [26] but, in this case, Kocuria was present only in a subset of individuals taking antibiotics, a majority of whom (16/17) were not admitted to hospital. This association was not determined to be significant, hence a larger study is needed to confirm the connection, if any.

Samples from patients who required breathing assistance showed a marked difference in the PCoA compared to those who did not and revealed an association with increased abundance of Serratia, Streptococcus, Enterobacter, Veillonella, Prevotella and Rothia. On the other hand, while a subset of COVID-19+ patients had Finegoldia, an opportunistic human pathogen, and Peptoniphilus, typically associated with the gut and vaginal microbiota, these genera were absent in COVID-19+ patients requiring breathing assistance. Breathing assistive devices were employed only for inpatients, although not all inpatients required breathing devices. Other investigators have noted ventilator-associated lung dysbiosis in COVID-19 patients with significantly higher rates of co-infection versus non-COVID-19 patients [27]. In patients with severe COVID-19 (n = 38 positive), that invasive mechanical ventilation greatly increased the risk of secondary infection, particularly bacterial infections in the respiratory system, bloodstream, and urinary tract [28]. The interplay between SARS-CoV-2 infection, breathing treatment, hospital stay, and respiratory microbiomes is complex, and a detailed timeline and/or longitudinal sampling of patients is required to identify causal relationships. Nonetheless, our results indicate that these factors should be considered when to evaluating microbiome structure and dynamics over time.

Reduction in Fusobacterium periodonticum associated with SARS-CoV-2 was reported in a study of Italian patients (n = 18 positive, 12 negative) [4], but was not observed in our study, nor did we observe significant difference in alpha diversity and COVID-19 status. Propionibacteriaceae, or Corynebacteriam accolens, as was reported for a study of 50 patients (n = 40 positive, 10 negative) in Baltimore [3]. Another study found no differences in the nasopharyngeal microbiota of patients who had mild COVID-19 (n = 9 positive, 10 negative) [2]. Analysis of patient samples in Nashville (n = 38 mild to moderate positive, 21 negative) revealed the complexity of interaction between viruses and bacteria in the upper respiratory tract, as differences in the airway microbiome were reported to be dependent on the SARS-CoV-2 viral load [16]. Similarly, in the lower airway, samples obtained through bronchial lavage (n = 142 positives) were enriched with oral bacteria in COVID-19 patients who had worse clinical outcomes, but mortality was better predicted by viral load and host immune response [17]. While microbial differences were present, SARS-CoV-2 and host factors were the most important in disease severity in both studies. Finally, Clostridium botulinum, Bacillus cereus and Halomonas spp. were found in almost all samples in an Indian study focused on more severe COVID-19 [5], which was clearly incompatible with our findings. This study did not include COVID-19 negative controls, which may have helped to distinguish between species commonly found in the NP space for the region in which they were collected versus those that are enhanced or lost during SARS-CoV-2 infection. In summary, viewed across the many studies to date, including the one reported here, differences in patient populations, geography, and climate, or even quantity or genomic diversity of the SARS-CoV-2 virus itself all play some role in the NP microbiome.

Reliance on retrospective medical chart review is recognized to be a limitation of the study reported here. At the time of data collection, direct patient contact or interviewing was difficult or impossible due to limitations of both PPE supply and staffing. It is well known that antibiotics dramatically change the gut microbiome [29]. Thus, future research will need to investigate specific medication used, with respect to the NP microbiome. Additionally, antibiotic use can cause changes in the gut biome for a few weeks or up to six months [30] and antibiotics taken by patients for up to six months earlier may also have had a lasting effect on the NP microbiome, which this study was not designed or powered to detect.

5. Conclusions

While COVID-19 infection was found to be associated with changes in the microbial communities comprising the NP microbiome, breathing treatments had a more varied and profound impact on the type and abundance of microbial taxa comprising the NP microbiome. These findings suggest the NP microbiome should be further evaluated longitudinally over the course of a hospital stay.