Characterization of the Upper Respiratory Bacterial Microbiome in Critically Ill COVID-19 Patients

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 K2 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.


Introduction
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,

Statistical and Correlation Network Analysis
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.

Genome Reconstruction and Functional Annotation of Moraxella Lincolnii
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].

Patient Data
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.

Composition and Alteration of URT Microbiota Taxa in COVID-19 Patients
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.

Distinct URT Microbiota Diversity in COVID-19 Patients
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)

Bacterial Microbiota Associated with Clinical Outcome in COVID-19 Patients
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.
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.

Genomic Feature and Functional Potential of Moraxella Lincolnii
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.

Discussion
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.

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

Informed Consent Statement:
Informed consent was obtained from the study participants.

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