key: cord-0681988-0rnm1mbl
authors: Bai, Xiangning; Narayanan, Aswathy; Skagerberg, Magdalena; Ceña-Diez, Rafael; Giske, Christian G.; Strålin, Kristoffer; Sönnerborg, Anders
title: Characterization of the Upper Respiratory Bacterial Microbiome in Critically Ill COVID-19 Patients
date: 2022-04-23
journal: Biomedicines
DOI: 10.3390/biomedicines10050982
sha: 04d16ddb06a813cd2c46cf3c4e0d5181c722030e
doc_id: 681988
cord_uid: 0rnm1mbl

The upper respiratory tract (URT) microbiome can contribute to the acquisition and severity of respiratory viral infections. The described associations between URT microbiota and severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection are limited at microbiota genus level and by the lack of functional interpretation. Our study, therefore, characterized the URT bacterial microbiome at species level and their encoded pathways in patients with COVID-19 and correlated these to clinical outcomes. Whole metagenome sequencing was performed on nasopharyngeal samples from hospitalized patients with critical COVID-19 (n = 37) and SARS-CoV-2-negative individuals (n = 20). Decreased bacterial diversity, a reduction in commensal bacteria, and high abundance of pathogenic bacteria were observed in patients compared to negative controls. Several bacterial species and metabolic pathways were associated with better respiratory status and lower inflammation. Strong correlations were found between species biomarkers and metabolic pathways associated with better clinical outcome, especially Moraxella lincolnii and pathways of vitamin K(2) biosynthesis. Our study demonstrates correlations between the URT microbiome and COVID-19 patient outcomes; further studies are warranted to validate these findings and to explore the causal roles of the identified microbiome biomarkers in COVID-19 pathogenesis.

The upper respiratory tract (URT) is the primary portal of entry for the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [1] , which has caused the coronavirus disease 2019 (COVID- 19) pandemic. Infection with SARS-CoV-2 may cause epithelial barrier dysfunction enhancing inflammatory responses and dysbiosis in the respiratory tract, which may worsen the pathogenic processes [2] . It has been evidenced that the URT microbiota may influence the susceptibility and severity of respiratory viral infections [3] . Co-infections of SARS-CoV-2 with other respiratory viruses and bacteria are well described in COVID-19 patients [4] [5] [6] [7] [8] [9] . Studies have claimed that bacterial co-infections are more frequently encountered in COVID-19 compared to other viral infections [10, 11] . However, this has not been confirmed and any increase in bacterial co-infections may rather be due to, e.g., the length of hospital stay and ventilation time. For other viral infections, 

The genomic DNA (gDNA) of the NP samples were extracted by the standardized International Human Microbiota Standards (IHMS) Protocol Q (http://www.microbiomestandards.org, accessed on 26 January 2021) [19] with some modifications. Sequencing libraries were prepared with the Nextera DNA Flex kit (Illumina, CA, USA) following the manufacturer's instructions. Libraries were normalized with Qubit assay, and then sequenced on NovaSeq6000 (NovaSeq Control Software 1.7.0/RTA v3.4.4) with a 151nt (Read1)-10nt(Index1)-10nt(Index2)-151nt(Read2) setup using 'NovaSeqXp' workflow in 'S4 mode flowcell. The Bcl to FastQ conversion was performed using bcl2fastq_v2.20.0.422 from the CASAVA software suite. The quality scale used is Sanger/phred33/Illumina 1.8+.

The raw sequencing data were pre-processed using our in-house bioinformatics pipeline as described previously [20] . Briefly, the adapter and low-quality reads (a quality score of less than Q30) were removed using Trim galore (v0.6.4) (https://www.bioinformatics. babraham.ac.uk/projects/trim_galore/, accessed on 30 March 2021). After the quality trimming, Bowtie2 (v2.3.5.1) [21] was used in combination with SAMtools (v1.19) [22] and BEDtools (v2.29.2) [23] to identify and remove human DNA sequences. The non-human reads were then used for downstream analysis. The bacterial taxonomic assignment and abundance estimation was conducted with MetaPhlAn 3.0 [24] using default parameters. Functional profiling was performed using the HMP Unified Metabolic Analysis Network 3 (HUMAnN 3.0), which quantifies gene families and microbial pathways in microbial community from metagenomic sequencing data [25] .

Alpha diversity of bacterial communities was assessed with microbial richness (number of detected taxa), Shannon and Simpson diversity indices using R function esti-mate_richness. Differences in alpha diversity between groups were assessed by testing the significance of these indexes using Wilcoxon rank sum test. Beta diversity was measured by Bray-Curtis and weighted UniFrac distances using R package Phyloseq (v1.30.0) [26] . Samples were clustered according to bacterial composition using non-metric multidimensional scaling (NMDS) approach with Bray-Curtis distance in Phyloseq (v1.30.0). Permutational multivariate analysis of variance (PERMANOVA) was performed to test the differences in bacterial composition between groups using vegan package (Adonis function) [27] using a Bray-Curtis dissimilarity method. Given the small sample size, different methods were used to determine and verify specific differences in bacterial taxa and metabolic pathways between groups. In addition to Wilcoxon rank sum test, LEfSe algorithm [28] was used to identify specific bacterial taxa and metabolic pathways as taxonomic and functional biomarkers. Kruskal-Wallis test was used to process the dataset with LEfSe alpha values set at 0.05. The threshold used to consider a discriminative feature for the logarithmic linear discriminant analysis (LDA) score was set at >2.

Correlation analyses were performed using the Spearman's rank correlation coefficient rho (library "psych", function "corr.test"). Correlation network of bacterial species, pathways, and clinical parameters was generated based on the Spearman's correlation coefficient. The input variables were species biomarkers, metabolic pathways, and clinical markers reflecting COVID-19 outcome. The integration network was constructed using R package bnlearn [29] , and only edges of correlation significance test were plotted. Visualization of the network was performed using Cytoscape (v3.6.1) [30] . Benjamini-Hochberg correction was used to adjust p-values in the case of multiple testing. Due to small sample size and exploratory purpose of this study, factors with adjusted p-value below 0.1 were considered statistically significant; whenever no significant association was identified after correction, results for unadjusted analysis were given, where raw p-value below 0.05 was considered significant.

To understand the genetic characteristics and functional potential of species biomarker (i.e., Moraxella lincolnii) that was found to be associated with clinical outcome, we performed genome reconstruction through metagenome assembly and functional prediction. Briefly, pre-processed reads were mapped to the reference genome of Moraxella lincolnii strain CCUG 9405 (GCA_002014765.1) using Bowtie2 (v2.3.5.1). The mapped reads belonging to Moraxella lincolnii were extracted using SAMtools (v1.19), then assembled using MEGAHIT v1.1.3 [31] and kmer lengths starting from 21 to 141. For further confirmation, the assemblies were mapped to the reference Moraxella lincolnii genome using BLASTn. The draft genome of Moraxella lincolnii was further annotated using RAST (Rapid Annotation using Subsystem Technology) Server (https://rast.nmpdr.org/rast.cgi, accessed on 13 July 2021), which provides high-quality gene calling and functional annotation including a mapping of genes to subsystems and metabolic reconstruction [32] .

NP samples from 37 critically ill COVID-19 patients were collected at a median (range) of 7 (2-35) days after symptom onset. The median (range) age of the patients was 61 years . The majority (n = 30/37, 81.1%) were males. The median (range) body mass index (BMI) was 29.96 (19.27-55) . Out of the 37 patients, 25 (67.6%) had other comorbidities including hypertension, diabetes, chronic lung disease, ischemic heart disease, heart failure, systemic inflammatory disease, transplanted, dementia, neurologic disease, and malignancy. Symptoms at admission to hospital included fever, cough, shortness of breath, chest pain, gastrointestinal problems, and loss of taste or smell. Six patients received antibiotics prior to sample collection. Twenty-nine patients were discharged, and eight patients were deceased.

The median (range) Ct value of E gene and SARS-CoV-2 specific RdRp/ORF1/N2 gene was 24.8 (13.8-39.2) and 23.9 (13.8-37.8), respectively. The median Ct value of RdRp/ORF1/N2 gene was used to define higher (Ct value ≤ 23.9, n = 17) and lower viral load (Ct value > 23.9, n = 20). Patients were further divided into different groups based on clinical variables reflecting COVID-19 outcome, i.e., SpO2 (91-100%, n = 6; 81-90%, n = 12; ≤80%, n = 12; supplemental oxygen support, n = 7), PaO2/FiO2 ratio (≥300 mm Hg, n = 6; 200-299 mm Hg, n = 17; <200 mm Hg, n = 14), respiratory SOFA score (0-1, n = 6; 2, n = 17; 3-4, n = 14). The characteristics of COVID-19 patients including levels of inflammatory markers CRP, lymphocytes, D-dimer, ferritin, and IL-6 of COVID-19 patients are shown in Table 1 and Supplementary Table S1 .

In total, 260 bacterial species belonging to 95 genera, 59 families, and 33 orders were identified ( Figure 1A and Supplementary Table S2 ). The most abundant bacterial genera were Cutibacterium, Corynebacterium, and Staphylococcus, among which Corynebacterium was significantly enriched in COVID-19 patients compared to controls (p = 0.0478, Wilcoxon rank sum test). Among other genera with relative abundance above 0.2%, four were significantly decreased in COVID-19 patients (p < 0.028, Wilcoxon rank sum test) ( Figure 1B) . The most abundant bacterial species were Cutibacterium acnes, Corynebacterium accolens, Corynebacterium pseudodiphtheriticum, and Staphylococcus aureus. Among 25 species with relative abundance of >0.2%, three were significantly decreased in patients (p < 0.018, Wilcoxon rank sum test) ( Figure 1C ). Using LEfSe, eight out of 95 genera, Ralstonia, Lactobacillus, Atopobium, Dialister, Porphyromonas, Slackia, Neisseria, and Rothia were significantly decreased in COVID-19 patients, and Corynebacterium was significantly enriched in patients (LDA score = 5.052, adjusted p = 0.059) ( Figure 1D ), which was consistent with the Wilcoxon rank sum test. At species level, five species were significantly decreased in COVID-19 patients (LDA score > 3.3, adjusted p < 0.056) ( Figure 1E ). Notably, we observed a high abundance of respiratory bacteria that commonly cause pneumonia in critically ill COVID-19 patients, e.g., Staphylococcus aureus, Haemophilus influenzae, and Moraxella catarrhalis ( Figure 1F ). Significantly differentially abundant genera and species between COVID-19 patients and negative controls are labelled with A-N as indicated, more details are shown in Figure 1D ,E. The size of each node represents their relative abundance. (B,C) Bar plots of main bacterial taxa at genus and species levels (average relative abundance > 0.2%) between patients and controls. * Statistically significant difference (Benjamini-Hochberg adjusted p < 0.06). (D,E) Taxonomic biomarkers at genus (D) and species level (E) identified by linear discriminative analysis (LDA) effect size (LEfSe) analysis between patients (in red) and controls (in green). LDA scores (log 10) for the enriched taxa in controls are represented on the positive scale, while LDA-negative scores indicate enriched taxa in patients. The LEfSe alpha value was set at 0.05, and the threshold used to consider a discriminative feature for the LDA score was set at >2. (F) Heat map of abundant bacterial species (average abundance > 0.2%) among individuals between patients and controls. The relative abundance of bacterial species is represented by color gradient as indicated. The species were ordered by decreasing relative abundance.

A significant reduction in alpha diversity of bacterial microbiota at genus level was found in NP samples from COVID-19 patients compared to those from SARS-CoV-2 negative controls, as measured by the microbial richness (adjusted p = 0.045, Wilcoxon rank sum test), Shannon and Simpson diversity indices (adjusted p = 0.034, Wilcoxon rank sum test) (Figure 2A ). Similarly, a decrease in microbial richness, Shannon and Simpson indices of alpha diversity at species level was observed in COVID-19 patients (adjusted p < 0.1, Wilcoxon rank sum test) ( Figure 2B ). No significant difference in beta diversity of bacterial microbiota was found between samples from patients and controls at genus or species level, as assessed with Bray-Curtis and weighted UniFrac dissimilarities (data not shown). NMDS based on Bray-Curtis distance showed no significant separation between patients and controls at genus or species levels (p > 0.05, PERMANOVA); however, controls were more diversely distributed than COVID-19 patients ( Figure 2C ,D).

To assess the potential effect of antibiotics on differences in the bacterial microbiota between groups, we excluded the six COVID-19 patients who had received antibiotics within three months prior to sampling. Similarly, we observed that the alpha diversity of bacterial microbiota at species level in COVID-19 patients was marginally significantly decreased compared to controls, as assessed with Shannon and Simpson diversity indices (p = 0.056 and 0.054, respectively, Wilcoxon rank sum test) (Supplementary Figure S1A ), while no difference in beta diversity was found between patients and controls (Supplementary Figure S1B ,C). The differentially abundant bacterial species between 31 patients who had not been given antibiotics and controls were consistent with those observed between all patients and controls with two exceptions (Supplementary Figure S1D) . As information about the use of antibiotics was unavailable in SARS-CoV-2-negative individuals, we did not adjust this factor in subsequent comparative analyses performed between patients and controls.

Given that the bacterial microbiota diversity was decreased in COVID-19 patients, we were interested to determine if samples with higher viral load (Ct value ≤ 23.9, n = 17) showed decreased microbiota diversity compared to those with lower viral load (Ct value > 23.9, n = 20), but no significant difference was observed (data not shown). No significant correlation was observed between the abundance of bacterial species and Ct values of the SARS-CoV-2 specific gene. shown). NMDS based on Bray-Curtis distance showed no significant separation between patients and controls at genus or species levels (p > 0.05, PERMANOVA); however, controls were more diversely distributed than COVID-19 patients ( Figure 2C,D) . To assess the potential effect of antibiotics on differences in the bacterial microbiota between groups, we excluded the six COVID-19 patients who had received antibiotics within three months prior to sampling. Similarly, we observed that the alpha diversity of bacterial microbiota at species level in COVID-19 patients was marginally significantly decreased compared to controls, as assessed with Shannon and Simpson diversity indices (p = 0.056 and 0.054, respectively, Wilcoxon rank sum test) (Supplementary Figure S1A) , while no difference in beta diversity was found between patients and controls (Supplementary Figure S1B,C) . The differentially abundant bacterial species between 31 patients who had not been given antibiotics and controls were consistent with those observed between all patients and controls with two exceptions (Supplementary Figure S1D) . As information about the use of antibiotics was unavailable in SARS-CoV-2-negative individuals, we did not adjust this factor in subsequent comparative analyses performed between patients and controls.

Given that the bacterial microbiota diversity was decreased in COVID-19 patients, we were interested to determine if samples with higher viral load (Ct value ≤ 23.9, n = 17) 

Patients were categorized into groups based on respiratory and inflammatory status reflecting the severity of COVID-19 outcome, as mentioned above. The association between bacterial microbiota and the clinical parameters was analysed by LefSe and Spearman's correlation analyses. LEfSe showed that Moraxella lincolnii and Propionibacterium namnetense were enriched in patients who had the highest PaO2/FiO2 ratio (≥300 mm Hg) or lowest respiratory SOFA score 0-1 (Moraxella lincolnii: LDA score = 4.307, adjusted p = 0.079; Propionibacterium namnetense: LDA score = 4.255, adjusted p = 0.0089, respectively) ( Figure 3A) . Moraxella lincolnii and Propionibacterium namnetense were also enriched in patients with the highest SpO2 91-100% (LDA score = 4.234, p = 0.014; LDA score = 4.077, adjusted p = 0.026, respectively) ( Figure 3B ). Spearman's correlation corroborated the association between microbiota and clinical markers obtained by LefSe. Several bacterial species correlated to certain clinical markers ( Figure 3C ), among which, Moraxella lincolnii and Propionibacterium namnetense correlated positively to PaO2/FiO2 ratio, while inversely to inflammation marker CRP and D-dimer, respectively (p < 0.05). These data imply that Moraxella lincolnii and Propionibacterium namnetense may be associated with improved clinical outcome, i.e., better respiratory status and lower inflammation level. respectively) ( Figure 3B ). Spearman's correlation corroborated the association between m crobiota and clinical markers obtained by LefSe. Several bacterial species correlated to ce tain clinical markers ( Figure 3C ), among which, Moraxella lincolnii and Propionibacteriu namnetense correlated positively to PaO2/FiO2 ratio, while inversely to inflammation mark CRP and D-dimer, respectively (p < 0.05). These data imply that Moraxella lincolnii and Pr pionibacterium namnetense may be associated with improved clinical outcome, i.e., better re piratory status and lower inflammation level. (SpO2) level or oxygen support. The biomarkers were identified by linear discriminative analysis (LDA) effect size (LEfSe) analysis. LDA scores (log 10) for the enriched species in a given group are represented with colors as shown. The LEfSe alpha value was set at 0.05, the threshold used to consider a discriminative feature for the LDA score was set at >2. (C) Correlation between bacterial species and markers of respiratory and inflammatory status. Spearman's correlation rho values are represented by color gradient as indicated (red is for positive, green is for negative correlation). Only correlations with p < 0.05 are shown on the plots.

Functional analysis was performed to understand the potential role of the bacterial microbiota in COVID-19 outcome. Intriguingly, several pathways correlated positively to PaO2/FiO2 ratio, while inversely to at least one inflammatory marker, i.e., CRP, Ddimer, or ferritin in COVID-19 patients (p < 0.05, Spearman's correlation) ( Figure 4A ). These included: (i) superpathways of menaquinol-7, menaquinol-11, menaquinol-12, and menaquinol-13 biosynthesis; (ii) mono-trans, poly-cis decaprenyl phosphate biosynthesis pathway; (iii) superpathway of tetrahydrofolate biosynthesis; (iv) gondoate biosynthesis (anaerobic); (v) flavin biosynthesis I; (vi) biotin biosynthesis II; (vii) L-histidine degradation II; (viii) polyisoprenoid biosynthesis. Interestingly, most of these pathways belong to the same superclass 'vitamin biosynthesis'; in particular, the four superpathways of menaquinone biosynthesis are known as vitamin K 2 biosynthesis. These results might indicate a potential role of vitamin K 2 in better clinical outcome, i.e., better respiratory status and lower inflammation in COVID- 19. way; (iii) superpathway of tetrahydrofolate biosynthesis; (iv) gondoate biosynthesis (anaerobic); (v) flavin biosynthesis I; (vi) biotin biosynthesis II; (vii) L-histidine degradation II; (viii) polyisoprenoid biosynthesis. Interestingly, most of these pathways belong to the same superclass 'vitamin biosynthesis'; in particular, the four superpathways of menaquinone biosynthesis are known as vitamin K2 biosynthesis. These results might indicate a potential role of vitamin K2 in better clinical outcome, i.e., better respiratory status and lower inflammation in COVID-19. Additionally, we observed that the PaO2/FiO2 ratio correlated positively to pathways involved in 1,4-dihydroxy-2-naphthoate biosynthesis I, formaldehyde oxidation I, etc. CRP level correlated inversely to pathways involved in pyrimidine deoxyribonucleotides de novo biosynthesis, heme biosynthesis I (aerobic), etc., while they correlated positively to pyruvate fermentation to propanoate I. The serum ferritin level correlated inversely to pathways involved in heme biosynthesis I (aerobic), 6-hydroxymethyl-dihydropterin diphosphate biosynthesis I, L-lysine biosynthesis I, etc. D-dimer correlated inversely to several functional pathways (p < 0.05, Spearman's correlation) ( Figure 4A ).

To assess the correlations between bacterial microbiota and metabolic pathways contributing to respiratory and inflammatory status, Spearman's correlation analysis was performed between bacterial species and functional pathways in COVID-19 patients. Interestingly, we observed very strong correlations between species and pathway biomarkers that were associated with clinical markers. Moraxella lincolnii correlated strongly to superpathways of menaquinol-7, menaquinol-11, menaquinol-12, and menaquinol-13 biosynthesis (rho = 1, adjusted p < 2.2 × 10 −16 ); mono-trans, poly-cis decaprenyl phosphate biosynthesis (rho = 0.998, adjusted p = 2.33 × 10 −43 ); formaldehyde oxidation I (rho = 0.805, adjusted p = 2.50 × 10 −7 ); 1,4-dihydroxy-2-naphthoate biosynthesis I (rho = 0.781, adjusted p = 1.23 × 10 −6 ), etc. ( Figure 4B ). Most of these pathways, as mentioned above, belong to the same superclass 'vitamin biosynthesis'. Propionibacterium namnetense correlated positively to the superpathway of pyrimidine deoxyribonucleotides de novo biosynthesis (rho = 0.657, adjusted p < 0.0001), and L-histidine degradation II (rho = 0.478, adjusted p = 0.0425) ( Figure 4B ).

Correlation network analysis was performed to unravel any interaction among two species biomarkers, metabolic pathways, and clinical markers reflecting COVID-19 outcome. Results demonstrated species-to-pathways, species-to-clinical markers, and pathway-toclinical markers interconnections, e.g., Moraxella lincolnii to superpathways of menaquinol biosynthesis, Moraxella lincolnii to CRP and the PaO2/FiO2 ratio, and superpathways of menaquinol biosynthesis to CRP and the PaO2/FiO2 ratio ( Figure 4C ). These data highlight the associations between microbiota and metabolic pathways contributing to clinical outcome in COVID-19 patients.

Given that Moraxella lincolnii was found to be associated with respiratory and inflammatory status, and that it also correlated very strongly to the metabolic pathways associated with better respiratory and inflammatory status, particularly vitamin K 2 biosynthesis, we were interested to identify genetic evidence that could support this finding. We therefore performed genome reconstruction and functional annotation of Moraxella lincolnii from one sample with a high abundance of this species. Genomic characteristics of M. lincolnii are summarized in Supplementary Table S3 . Notably, a gene that encodes a key enzyme in menaquinone (vitamin K 2 ) biosynthesis was identified in the M. lincolnii genome, i.e., ubiE (bifunctional demethylmenaquinone methyltransferase/2-methoxy-6-polyprenyl-1,4benzoquinol methylase). Additionally, M. lincolnii possess genes involved in several crucial metabolic pathways that are known to act in an integrated manner to maintain the balance and organism homeostasis, including genes involved in lipid metabolism, amino acid metabolism, glycolysis, pentose phosphate pathway, etc. (Supplementary Table S3b and Supplementary Figure S2) . These results support our hypothesis that M. lincolnii contributes to vitamin K 2 biosynthesis and other metabolic pathways, thereby, possibly being associated with better clinical outcome in COVID-19 patients. Further in vitro studies are in need to validate the effect of M. lincolnii in vitamin K 2 biosynthesis and other biological functions that may play a beneficial role in COVID-19.

In this study, using shotgun whole metagenome sequencing, we characterized the upper respiratory microbiome profile in critically ill COVID-19 patients and correlated the findings to viral load, and respiratory and inflammatory status. No association was found between SARS-CoV-2 loads and bacterial microbiota diversity or differentially abundant bacterial taxa, which is in line with a previous report using metagenome sequencing [14] . The viral load in the URT is highest early in the disease course [33, 34] , and since our sampling occurred a median of 7 days after symptom onset, the peak of viral replication had most probably passed. Intriguingly, in the COVID-19 patients, several bacterial species were associated with markers of both respiratory and inflammatory status; in particular, Moraxella lincolnii and Propionibacterium namnetense were correlated to better respiratory status and low inflammation. M. lincolnii, a poorly characterized bacterium isolated from the human respiratory tract [35] , was found in two of the COVID-19 patients who showed high viral load but with very low levels of inflammatory markers and normal respiratory conditions, implying that M. lincolnii might increase the host immunity against SARS-CoV-2. We acknowledged that this hypothesis remains to be confirmed since only two patients harbored M. lincolnii. It is noteworthy that a recent study has indeed suggested a protective role of M. lincolnii in respiratory health status [36] . Moreover, M. lincolnii has recently been shown to exhibit strong inhibitory activity against nasal S. aureus, mediated by proteins fitting the profile of antimicrobial peptides (AMPs) [37] . AMPs are known to exhibit antiviral and immunomodulatory properties [38] , which possibly could contribute to improved patient outcomes in viral respiratory infections. Research on the recently identified P. namnetense [39] is also very rare, our data appeal for further in-depth studies to investigate the potential role of these two species in COVID-19 outcomes.

To date, the most probable hypothesis on how the respiratory microbiome could influence viral respiratory infections relies on the immunological properties of microbes inhabiting the respiratory tract [3] . At functional level, we found that several metabolic pathways correlated to both better respiratory status and to lower inflammation level. Interestingly, most of these pathways are involved in vitamin biosynthesis, particularly four superpathways of menaquinone (vitamin K 2 ) biosynthesis. Several other pathways that correlated to better respiratory status or/and lower inflammation level belonged to pathway superclass 'Cofactor, Carrier, and Vitamin Biosynthesis', such as 1,4-dihydroxy-2naphthoate biosynthesis I, which is the naphthalenic intermediate in the biosynthesis of vitamin K 2 [40] . Our results may, thus, indicate a potential role of vitamin K 2 in improving clinical outcome in COVID-19. This is supported by the reported correlations between vitamin K deficiency and severe COVID-19 outcome [41, 42] . Vitamin K is associated with an impaired production of inflammatory cytokines and plays an important role in immunomodulation [43, 44] . In addition to attenuating the excessive production of proinflammatory cytokines [45] , vitamin K may protect the integrity of the alveolar-capillary membrane [46] , thereby, possibly improving respiratory status in COVID-19 patients as we observed.

Another remarkable finding was the strong correlation between the two species biomarkers (M. lincolnii, P. namnetense) and metabolic pathways associated with better respiratory status and lower inflammation level. In particular, M. lincolnii correlated strongly to pathways involved in vitamin K 2 biosynthesis. Vitamin K naturally occurs in two biologically active forms, K 1 and K 2 ; of these, vitamin K 2 is predominantly of bacterial origin [47, 48] . In order to confirm this and gain insights into the functions of M. lincolnii, we reconstructed the draft genome of M. lincolnii from metagenome data and performed functional annotation. Strikingly, we found that M. lincolnii carried the gene encoding 'bifunctional demethylmenaquinone methyltransferase/2-methoxy-6-polyprenyl-1,4-benzoquinol methylase', an enzyme catalyzing the last step in menaquinone (vitamin K 2 ) biosynthesis. M. lincolnii possess additional genes involved in several metabolic pathways that have been suggested to be associated with COVID-19, e.g., lipid and amino acid metabolism, heme biosynthesis, glycolysis, pentose phosphate pathway, etc. [49, 50] . Our data indicate a great need for in vitro research to validate the effects of M. lincolnii in vitamin K 2 biosynthesis and other metabolic processes that may play beneficial roles in COVID-19 outcome.

Our study demonstrated a reduced bacterial microbiota diversity in the critical COVID-19 patients compared to the SARS-CoV-2-negative individuals. Some respiratory pathogens, particularly pneumonia-causing bacteria, e.g., Staphylococcus aureus and Haemophilus influenzae, were abundant in the critically ill COVID-19 patients, highlighting the possibility that co-infections with such pathogens may contribute to severe clinical outcome. In contrast, a significant reduction in respiratory commensals, e.g., Neisseria mucosa, Ralstonia pickettii, was found in the COVID-19 patients. Such reductions in healthy commensals might contribute to the susceptibility to and severity of SARS-CoV-2 infection, as suggested in other viral infections [3] , although the causal relationship between URT microbiome alteration and SARS-CoV-2 infection warrants further investigation. It should be noted that most of our critically ill COVID-19 patients had well-known risk factors for disease severity, such as age, gender, and comorbidities [51, 52] , which may be confounders contributing to URT microbiome changes in the patients. Moreover, although the controls were sampled due to suspicions of SARS-CoV-2 infection, we could not rule out potential confounders that may explain the differences in the URT microbiome observed.

This study has limitations. The major flaws were the small sample size and the lack of detailed information of SARS-CoV-2-negative controls; thus, the microbiome changes in COVID-19 patients and microbiome biomarkers identified for clinical outcomes remain to be validated with further studies. Second, only critically ill COVID-19 patients were included, and a single NP sample per patient collected at hospital admission was analysed, our findings may not apply to patients with asymptomatic, mild to moderate COVID-19. Third, although whole metagenome sequencing has obvious advantages, particularly its functional profiling capacity, it suffers from host-derived DNA contamination, which may obscure microbial signatures in low-biomass and highly host-contaminated NP samples. 16S rRNA sequencing should be considered in combination with whole metagenome sequencing to obtain a comprehensive landscape of URT microbial communities and functionality. In spite of these limitations, our study reveals important information for the interpretation of the role of the URT microbiome in SARS-CoV-2 infection. It is noteworthy that the correlations observed in this study do not illustrate a direct causal link between the URT microbiome and SARS-CoV-2 infection as described widely in other microbiome studies [53] [54] [55] ; further in-depth studies are warranted to explore the causal roles of the URT microbiome in the COVID-19 pathogenesis.

In conclusion, our study characterized the URT microbiome in correlation to COVID-19 outcomes. Several bacterial species and metabolic pathways were associated with respiratory and inflammation status in COVID-19 patients. Strong associations were found between two species biomarkers and several pathways that were associated with better clinical outcome; in particular, Moraxella lincolnii and pathways involved in vitamin K 2 biosynthesis. To our knowledge, this is the first study to depict the URT microbiome associated with respiratory status in critical COVID-19 patients. In addition, our study demonstrates a distinct URT microbiome profile in patients with critical COVID-19 compared to non-COVID-19 individuals. These findings aggregately render evidence of the URT microbiome as a possible contributor to COVID-19 outcome. Future in-depth studies are warranted with a larger sample size, serial samples from each patient, samples from other geographic areas, combination of different sequencing techniques, and in vitro assays to elucidate the causal roles of URT microbiome changes in SARS-CoV-2 infection, disease progression, and patient outcomes. This could possibly aid the identification of microbial targets for potential interventions and treatments of COVID-19.

Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/biomedicines10050982/s1, Table S1 . Clinical and laboratory metadata of COVID-19 patients (.xlsx). Table S2 . Bacterial species identified in this study and their relative abundance in each sample (.xlsx). Table S3 . Genomic characteristics and functional annotation of Moraxella lincolnii (.xlsx). Figure S1 . Differences in bacterial microbiota diversity and differentially abundant species between 31 COVID-19 patients who did not take antibiotics prior to sampling and 20 SARS-CoV-2-negative controls. (A) Alpha diversity between COVID-19 patients and SARS-CoV-2-negative controls at bacterial species level assessed by richness, Shannon and Simpson diversity indices. (B,C) Bray-Curtis dissimilarity and non-metric multidimensional scaling (NMDS) based on Bray-Curtis distance of bacterial species composition between patients and controls. (D) Bacterial species biomarkers identified by linear discriminative analysis (LDA) effect size (LEfSe) analysis between 31 COVID-19 patients (in red) and 20 controls (in green). LDA scores (log 10) for the enriched species in controls are represented on the positive scale, while LDA-negative scores indicate enriched species in patients. The LEfSe alpha values was set at 0.05, the threshold used to consider a discriminative feature for the LDA score was set at >2. Figure S2 . Subsystem distribution of Moraxella lincolnii identified in this study. The functional annotation and subsystem categorization was conducted using RAST server. Funding: This study was supported by the Swedish Research Council (2020-02129, 2017-05848, 2016-01675) and ALF Med 20 FoUI-953887. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

The study was approved by the Swedish ethical review authority (Dnr 2020-01558, approval date: 22 April 2020).

Informed consent was obtained from the study participants.

The raw shotgun metagenome sequencing data have been submitted to the NCBI Sequence Read Archive under the accession number PRJNA781460.

A Novel Coronavirus from Patients with Pneumonia in China

The role for the metagenome in the pathogenesis of COVID-19

Relationship between nasopharyngeal microbiota and patient's susceptibility to viral infection

Co-infections of SARS-CoV-2 with multiple common respiratory pathogens in infected patients

Coinfection and Other Clinical Characteristics of COVID-19 in Children

Bacterial and viral co-infections in patients with severe SARS-CoV-2 pneumonia admitted to a French ICU

Bacterial co-infections with SARS-CoV-2

CoV-2-related pneumonia cases in pneumonia picture in Russia in March-May 2020: Secondary bacterial pneumonia and viral co-infections

Metagenomic Analysis Reveals Clinical SARS-CoV-2 Infection and Bacterial or Viral Superinfection and Colonization

Epidemiology and clinical features of COVID-19: A review of current literature

Co-infections in people with COVID-19: A systematic review and meta-analysis

The role of the local microbial ecosystem in respiratory health and disease

Nasopharyngeal Microbiota Profiling of SARS-CoV-2 Infected Patients

Metagenomic Next-Generation Sequencing of Nasopharyngeal Specimens

Metatranscriptomic Characterization of COVID-19 Identified A Host Transcriptional Classifier Associated With Immune Signaling

Nasopharyngeal Microbiome Signature in COVID-19 Positive Patients: Can We Definitively Get a Role to Fusobacterium periodonticum? Front. Cell Infect

Temporal association between human upper respiratory and gut bacterial microbiomes during the course of COVID-19 in adults

Signatures of COVID-19 Severity and Immune Response in the Respiratory Tract Microbiome

Towards standards for human fecal sample processing in metagenomic studies

Fast gapped-read alignment with Bowtie 2

Genome Project Data Processing, S. The Sequence Alignment/Map format and SAMtools

BEDTools: A flexible suite of utilities for comparing genomic features

Integrating taxonomic, functional, and strain-level profiling of diverse microbial communities with bioBakery 3

Species-level functional profiling of metagenomes and metatranscriptomes

phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data

Community Ecology Package. 2020. Available online

Metagenomic biomarker discovery and explanation

Learning Bayesian Networks with the bnlearn R Package

Cytoscape: A software environment for integrated models of biomolecular interaction networks

An ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph

The RAST Server: Rapid annotations using subsystems technology

COVID-19: Recommended sampling sites at different stages of the disease

Temporal profiles of viral load in posterior oropharyngeal saliva samples and serum antibody responses during infection by SARS-CoV-2: An observational cohort study

Moraxella lincolnii sp. nov., isolated from the human respiratory tract, and reevaluation of the taxonomic position of Moraxella osloensis

Differential nasopharyngeal microbiota composition in children according to respiratory health status

Identification of Nasal Gammaproteobacteria with Potent Activity against Staphylococcus aureus: Novel Insights into the

Human Antimicrobial Peptides as Therapeutics for Viral Infections

Propionibacterium namnetense sp. nov., isolated from a human bone infection

Biosynthesis of bacterial menaquinones: The membrane-associated 1,4-dihydroxy-2-naphthoate octaprenyltransferase of Escherichia coli

Reduced vitamin K status as a potentially modifiable risk factor of severe COVID-19

The Association of Low Vitamin K Status with Mortality in a Cohort of 138 Hospitalized Patients with COVID-19

Vitamin K suppresses the lipopolysaccharide-induced expression of inflammatory cytokines in cultured macrophage-like cells via the inhibition of the activation of nuclear factor kappaB through the repression of IKKalpha/beta phosphorylation

Inhibition of TNF-alpha, IL-1alpha, and IL-1beta by Pretreatment of Human Monocyte-Derived Macrophages with Menaquinone-7 and Cell Activation with TLR Agonists In Vitro

The COVID-19 Cytokine Storm; What We Know So Far. Front

Potential Beneficial Effects of Vitamin K in SARS-CoV-2 Induced Vascular Disease?

Vitamin K: Food composition and dietary intakes

The importance of menaquinones in human nutrition

Metabolic Alterations in SARS-CoV-2 Infection and Its Implication in Kidney Dysfunction

Immune and Metabolic Signatures of COVID-19 Revealed by Transcriptomics Data Reuse. Front

Risk factors for severe and critically ill COVID-19 patients: A review

Risk factors for Covid-19 severity and fatality: A structured literature review

Targeting the Gut Microbiota in Coronavirus Disease 2019: Hype or Hope?

Gut microbiome alterations and gut barrier dysfunction are associated with host immune homeostasis in COVID-19 patients

The human microbiome and COVID-19: A systematic review

We acknowledge support from the National Genomics Infrastructure in Genomics Production Stockholm funded by Science for Life Laboratory, the Knut and Alice Wallenberg Foundation and the Swedish Research Council, and SNIC/Uppsala Multidisciplinary Center for Advanced Computational Science for assistance with massively parallel sequencing and access to the UPPMAX computational infrastructure. We thank the core facility at Novum, BEA, Bioinformatics and Expression Analysis, which is supported by the board of research at the Karolinska Institute and the research committee at the Karolinska University Hospital. The computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at (SNIC CENTRE) partially funded by the Swedish Research Council through grant agreement no. 2018-05973. We also would like to thank Robert van Domselaar from the Division of Infectious Diseases, Department of Medicine Huddinge, Karolinska Institute for intellectual input and laboratory assistance.

The authors declare no conflict of interest.