Dysbiosis and structural disruption of the respiratory microbiota in COVID-19 patients with severe and fatal outcomes

The COVID-19 outbreak has caused over three million deaths worldwide. Understanding the pathology of the disease and the factors that drive severe and fatal clinical outcomes is of special relevance. Studying the role of the respiratory microbiota in COVID-19 is especially important as the respiratory microbiota is known to interact with the host immune system, contributing to clinical outcomes in chronic and acute respiratory diseases. Here, we characterized the microbiota in the respiratory tract of patients with mild, severe, or fatal COVID-19, and compared it to healthy controls and patients with non-COVID-19-pneumonia. We comparatively studied the microbial composition, diversity, and microbiota structure between the study groups and correlated the results with clinical data. We found differences in the microbial composition for COVID-19 patients, healthy controls, and non-COVID-19 pneumonia controls. In particular, we detected a high number of potentially opportunistic pathogens associated with severe and fatal levels of the disease. Also, we found higher levels of dysbiosis in the respiratory microbiota of patients with COVID-19 compared to the healthy controls. In addition, we detected differences in diversity structure between the microbiota of patients with mild, severe, and fatal COVID-19, as well as the presence of specific bacteria that correlated with clinical variables associated with increased risk of mortality. In summary, our results demonstrate that increased dysbiosis of the respiratory tract microbiota in patients with COVID-19 along with a continuous loss of microbial complexity structure found in mild to fatal COVID-19 cases may potentially alter clinical outcomes in patients. Taken together, our findings identify the respiratory microbiota as a factor potentially associated with the severity of COVID-19.

The COVID-19 outbreak has caused over three million deaths worldwide. Understanding the pathology of the disease and the factors that drive severe and fatal clinical outcomes is of special relevance. Studying the role of the respiratory microbiota in COVID-19 is especially important as the respiratory microbiota is known to interact with the host immune system, contributing to clinical outcomes in chronic and acute respiratory diseases. Here, we characterized the microbiota in the respiratory tract of patients with mild, severe, or fatal COVID-19, and compared it to healthy controls and patients with non-COVID-19-pneumonia. We comparatively studied the microbial composition, diversity, and microbiota structure between the study groups and correlated the results with clinical data. We found differences in the microbial composition for COVID-19 patients, healthy controls, and non-COVID-19 pneumonia controls. In particular, we detected a high number of potentially opportunistic pathogens associated with severe and fatal levels of the disease. Also, we found higher levels of dysbiosis in the respiratory microbiota of patients with COVID-19 compared to the healthy controls. In addition, we detected differences in diversity structure between the microbiota of patients with mild, severe, and fatal COVID-19, as well as the presence of specific bacteria that correlated with clinical variables associated with increased risk of mortality. In summary, our results demonstrate that increased dysbiosis of the respiratory tract microbiota in patients with COVID-19 along with a continuous loss of microbial complexity structure found in mild to fatal COVID-19 cases may potentially alter clinical outcomes in patients. Taken together, our findings identify the respiratory microbiota as a factor potentially associated with the severity of COVID-19.
The Coronavirus Disease 2019 (COVID-19) outbreak, declared a pandemic by the World Health Organization on March 11, 2020, is caused by the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2). As of May 2021, SARS-CoV-2 has infected more than 150 million people and caused over three million deaths worldwide 1 . COVID-19 shows a broad spectrum of clinical manifestations ranging from asymptomatic infection and mild respiratory symptoms to severe pneumonia and death 2,3 and has been linked to demographic factors and comorbidities 4,5 . To date, aberrant immune responses against SARS-CoV-2 antigens have been shown to be critically involved in severe clinical outcomes and other secondary inflammatory conditions that remain after initial COVID-19 infection 3,6 .
Studying the role of the human microbiota in COVID-19 is particularly relevant, as the respiratory microbiota is known to interact with the host immune system, contributing to clinical outcomes in chronic and acute respiratory diseases 7 . The respiratory microbiota plays a central role in shaping pulmonary immunity by enhancing innate and adaptive immune responses. This suggests that host immunity is regulated by interactions with bacterial communities in the respiratory tract.
Some studies suggest that interactions between microorganisms and the host immune system are speciesspecific, denoting that even small variations in the diversity and composition of the microbiota could have important consequences on host health [7][8][9] . In the case of COVID-19, severe to fatal clinical outcomes are often associated with the presence of comorbidities known to exhibit an altered (dysbiotic) microbiota 10 (e.g., diabetes type II, obesity, age, and heart disease). Furthermore, in a wide range of microbiome-associate diseases (MADs), dysbiosis is a common feature that may impact disease progression 11,12 . Nonetheless, few studies that characterize the respiratory microbiota in COVID-19 and the presence of dysbiosis are available to date [13][14][15][16][17][18][19][20] .
To gain insight into the association between respiratory microbiota and COVID-19 severity, we characterized the microbiota in the respiratory tract of patients with mild, severe, or fatal COVID-19, and compared it with healthy controls and patients with non-COVID-19-pneumonia. We performed comprehensive analyses to evaluate the respiratory microbiota as a risk factor and its implications for patients with COVID-19. We comparatively studied the microbial composition, diversity, and structure of the microbiota between the study groups and correlated the results with clinical data. These analyses let us detect firstly, differences in bacterial abundance between groups, secondly, higher levels of dysbiosis in the respiratory microbiota of patients with COVID-19, thirdly, differences in diversity structure between the microbiota of patients with mild, severe, and fatal COVID-19, and lastly, the presence of specific bacteria that correlated with clinical variables associated with increased risk of mortality. In summary, our results demonstrate that increased dysbiosis of the respiratory tract microbiota in patients with COVID-19, along with a continuous loss of microbial complexity structure found in mild to fatal COVID-19 cases, may potentially alter clinical outcomes in patients. Taken together, our findings identify the respiratory microbiota as a factor potentially associated with the severity of COVID-19.
Demographics, health-related characteristics, and symptomatology are described in Table 1. Overall, 52 patients were male (54.7%) with a median age of 45 years old (IQR: 21). Regarding health conditions, 58.2% of the participants had at least one comorbidity, with DM2 (17%), hypertension (17%), smoking (17%), and obesity (35%) being the most represented in the cohort. The median days of symptom onset was seven, and 52.6% of the individuals received antibiotic treatment prior to hospitalization. In addition, we found important associations between some health/demographic characteristics and severity. For example, patients with fatal COVID-19 were predominantly male (73.6%, Wilcoxon rank-sum test p = 0.01), significantly older (median = 58, Wilcoxon ranksum test p = 6.57e−07), with higher BMI (median = 30.4, Wilcoxon rank-sum test p = 0.05), and most of them received prior antibiotic treatment (78.9%, Wilcoxon rank-sum test p = 0.002) compared to patients with severe and mild COVID-19. Also, a higher number of days after symptoms onset was found in the non-COVID-19 pneumonia group (median = 10, Wilcoxon rank-sum test p = 0.01).

Respiratory microbiota composition differs between severity levels of COVID-19 patients and controls.
Of the 95 analyzed samples belonging to the upper respiratory tract, we identified a total of 4514 Amplicon Sequence Variants (ASVs). Regarding the analysis of the relative abundance at the phylum level (Fig. 1A), Firmicutes, Bacteroidetes, and Proteobacteria were the most dominant phyla among our severity groups and controls. In general, these phyla are present in all group samples but there were changes in relative abundance associated with disease severity. Overall, we found Firmicutes, Actinobacteria, and TM7 to be increased in patients with COVID-19, while Bacteroidetes and Proteobacteria were found to be decreased.
Analysis of relative abundance at the genus level revealed genera that differed significantly between patients with COVID-19 patients and controls (Fig. 1B). Overall, we found Veillonella, Staphylococcus, Corynebacterium, Neisseria, Actinobacillus, and Selenomonas enriched in patients with COVID-19 but reduced in the healthy controls. Conversely, we found Haemophilus and Alloiococcus enriched in the healthy controls but reduced in patients with COVID-19. In addition, there were differences in the abundance of some genera between the severity levels for COVID-19. For example, Streptococcus and Staphylococcus showed increasing abundance from mild to fatal COVID-19. By contrast, Haemophilus and Actinomyces showed the opposite pattern, where the highest abundance was associated with mild COVID-19 and the lowest with fatal COVID-19. In addition, we found Corynebacterium highly abundant only in severe COVID-19, whereas Actinobacillus was found highly abundant only in fatal and mild COVID-19.
Additionally, we also characterized the lower respiratory tract microbiota between patients with severe and fatal COVID-19, finding a different composition from that found in the upper respiratory tract ( Supplementary  Fig. S1). We found that fatal patients showed a higher abundance of Proteobacteria and Bacteroidetes compared Table 1. Demographic data of the cohort. Only upper respiratory samples (OPS and NPS) are included in the information. BMI Body Mass Index, DM2 Diabetes Mellitus Type 2. *Respiratory diseases: either asthma, COPD, or ILD. P values denote statistical significant differences given by Wilcoxon rank-sum test (*< 0.05, **< 0.005, ***< 0.0005).  www.nature.com/scientificreports/ Beta diversity. Beta diversity analyses showed differences in the microbiota composition between patients with different severity levels of COVID-19 and controls (Fig. 3A,B, Supplementary Fig. S2). In particular, the PCoA analysis ( Fig. 3A) showed differences between severity levels and control groups which are supported by the PERMANOVA result (F = 2.7, p = 0.007). We observe changes in the direction of the ellipses belonging to mild COVID-19 and healthy controls. While the ellipses for all severe, fatal COVID-19 and non-COVID-19 pneumonia occurred in the same direction, the ellipse for mild COVID-19 was almost orthogonal. Furthermore,  www.nature.com/scientificreports/ when analyzing the results of the Ružička metric we detected that the intra-treatment similarities of the healthy controls were significantly higher than the similarities between the diseased samples, meaning that the COVID-19 associated microbiota (regardless of severity level) exhibited significantly higher levels of dysbiosis compared to the healthy controls ( Supplementary Fig. S2). The LefSe analysis allowed us to identify the differentially abundant taxa associated with the groups that were compared (Fig. 3B). We observed that all the severe COVID-19 groups and the two control groups showed differentially abundant taxa or biomarkers. In particular, for mild COVID-19, we found Prevotella melaninogenica and P. pallens, Veillonella parvula, Neisseria subflava, Fusobacterium, and Actinomyces to be highly abundant. In the case of severe COVID-19, we found Megasphaera and CW040 as the most prevalent. In the case of fatal COVID-19, Rothia dentocariosa, Streptococcus infantis, and Veillonella dispar were the most significant. In addition, the highest number of differentially abundant taxa was found to be associated with healthy controls (e.g., Streptococcus, Flavobacterium, and Oribacterium, and f_Veillonellaceae). Finally, for the non-COVID-19-pneumonia group, we found Corynebacterium, Prevotella nigrescens, Capnocytophaga, and Enterobacteriaceae to be the most abundant.
The Lefse analysis allowed us to detect enriched or depleted bacteria in the different risk factor groups for the variables analyzed (Fig. 4B). We found depleted Neisseria subflava in the high-risk samples for troponin and APACHE. On the other hand, Veillonella dispar was found interestingly depleted in the low-risk samples for APACHE, BUN, myoglobin, and urea. However, we also found some bacterial groups to be consistently enriched in the high-risk samples for several clinical variables. For example, Corynebacterium was found to be enriched in the high-risk samples for lymphocytes count and urea, while Actinomyces was enriched for BUN and urea. Additionally, four ASV´s of the genus Prevotella (Prevotella melaninogenica; Prevotella; [Prevotella]; [Prevotella]_s) were found significantly enriched in the high-risk samples for myoglobin, BUN, troponin, and lymphocyte count.
The structure of the respiratory microbiota is different between patients with different severity levels of COVID-19. Network analysis for the microbiota associated with the severity levels for COVID-19 revealed differences at the structural level (Fig. 5A,B). Graphical representation of the networks showed a different arrangement for each and a continuum of loss of complexity across COVID-19 severity groups (from mild to fatal) (Fig. 3A). In particular, the microbiota network associated with mild COVID-19 was the largest and most connected (nodes = 148, edges = 4758) compared to severe COVID-19 (nodes = 84, edges = 688) and fatal www.nature.com/scientificreports/ COVID-19 (nodes = 74, edges = 75). Interestingly, in patients with fatal outcomes, the respiratory microbiota network was highly disaggregated and poorly connected with multiple isolated nodes (nodes = 74, edges = 75). The metric calculation of the networks illustrates that the topology associated with each COVID-19 severity level is different (Fig. 3B, Supplementary Table S1). For example, the mild COVID-19 network showed the highest values for the average number of neighbors, density, and clustering. By contrast, the severe disease network was characterized by higher centralization and heterogeneity, whereas the fatal disease network showed the highest values for diameter, characteristic path length, and connected components.

Discussion
Due to the massive efforts of the global research community, many investigations on the role of the microbiota in COVID-19 have appeared. Despite the clear focus on gut microbiota 6,21,22 there have been some studies that have characterized the respiratory microbiota and its impact on COVID-19. However, most of them focus on comparing patients with COVID-19 to those without COVID-19 [13][14][15][16][17][18][19][20] , and do not include disease severity levels thus hindering information that could potentially help to understand disease progression. In turn, to date, few studies have characterized the respiratory microbiota in COVID-19 among severity levels 23,24 . Here, we analyzed the respiratory microbiota from an eco-evolutionary perspective, by testing microbial composition and community structure of a large cohort of patients with different levels of severity (mild, severe, and fatal) and linked the results to clinical variables to gain insight into the mechanisms by which the microbiota may affect host response against the disease.
As in recent studies investigating the etiology of COVID-19 3,5 , we found that demographic and health-related factors showed strong associations with severity. Male sex, high values for BMI, age over 50 years, and previous antibiotic treatment were significantly associated with patients with fatal COVID-19 (Table 1), which potentially favor the development of a fatal state of the disease.
Similarly, in previous work exploring the respiratory microbiota associated with COVID-19 24-26 , we found significantly lower microbial diversity in the microbiota of COVID-19 patients than in the healthy controls (Fig. 1A). This result is relevant since, in general, a higher diversity correlates with a better response of microbial systems A B www.nature.com/scientificreports/ against perturbations (e.g., disease). A more diverse microbiota may persist after disease (e.g., by functional redundancy) or may recover to a previous state (e.g., resilience) 27 , having direct consequences on host health 11 . Furthermore, we found differences in the abundance of some bacteria between our study groups (Fig. 1A). In particular, as in other respiratory diseases 28 , we observed a higher ratio of Firmicutes/Bacteroidetes in patients with COVID-19. Firmicutes was detected highly abundant while Bacteroidetes was especially decreased in the microbiota associated with COVID-19 patients. This is of particular interest since, in murine models, it has been shown that Bacteroidetes can down-regulate the ACE2 expression 29 . Although the correlation was observed in the gut microbiota, the low abundance of members of that phylum in the severe and fatal patients in this study opens the question of whether this process may also take place in the respiratory tract.
Regarding the composition at the genus level, most of the differences we found were in potentially pathogenic bacteria (Fig. 1B). A gradual increase of Streptococcus was identified from mild to fatal COVID-19 cases. Although Streptococcus is usually a commensal member of the respiratory microbiota, it can become pathogenic in the face of environmental disturbances. In higher abundance, this genus has been associated with viral acute respiratory infections 30,31 . In addition, genera such as Veillonella, Staphylococcus, and Actinomyces also showed high abundances at different COVID-19 severity levels. Specifically, Veillonella and Actinomyces have been found as opportunistic pathogens in COVID-19 32 . Additionally, Staphylococcus is one of the most common causative agents of secondary infections in respiratory diseases such as influenza 33 .
Regarding the characterization of the microbiota in the lower respiratory tract, we detected opposite patterns in the abundances of some phyla and genera between severe and fatal COVID-19 ( Supplementary Fig. S1), suggesting that the phenomena that we observe in the upper tract, due to particular environmental conditions could be different in the lungs. For instance, respiratory diseases generally entail inflammatory processes that increase mucus production which in turn, favor the presence of biofilms and anaerobic bacteria in the lungs 34 . This could be the reason behind the depletion of aerobic bacteria such as Staphylococcus and the increase of anaerobic groups such as Streptococcus, Abiotrophia, and Mycoplasma in the samples of the lower respiratory tract.
In terms of beta diversity, we found some differences between the analyzed groups. For example, we observed in the PCoA analysis that the samples belonging to severe and fatal COVID-19, as well as to the non-COVID-19-pneumonia group, are more similar in terms of microbial composition, than those from healthy controls and mild COVID-19 patients (Fig. 2A). Moreover, the PCoA analysis displayed a change in the tendency of variation between some of the analyzed groups. This observation implies that the features that define the microbial variation within each group are distinct, at least between mild COVID-19/healthy controls and severe forms of the disease, including non-COVID-19 pneumonia. Furthermore, our dysbiosis analysis let us detect that COVID-19 microbiota showed higher levels of dysbiosis than from the healthy controls ( Supplementary Fig. S2). Several MADs exhibit this behavior in which microbiota instability (dysbiosis) is present not by the dominance of one or a few bacteria but because of a higher heterogeneity/stochasticity of microbial groups 12 . Dysbiosis involves a disruption in the bidirectional interactions between the host immune system and the microbial communities, which    www.nature.com/scientificreports/ can alter the functions provided by these communities and reshape the entire host-microbiota interaction 11,35 . Microbiota stability has been shown to be a hallmark of host health and homeostasis 9,11,31 . For example, some reports suggest that there is a homeostatic mechanism that maintains the lung epithelium in a state of interferon primacy with antiviral activity against other respiratory infections such as influenza. This particular antiviral response is stimulated by specific pathogen-associated molecular patterns (PAMPs) that are induced by microbial communities 36 . Furthermore, in other viral diseases it has been shown that, due to different mechanisms (e.g. alteration of the epithelium and increased adhesion of respiratory pathogens), microbial dysbiosis can promote viral infection, potentially including SAR-CoV-2 24 . A common question when studying MADs is whether dysbiosis favors the disease or is caused by it. In the case of COVID-19, the clinical outcome is highly correlated with comorbidities such as hypertension, diabetes, and obesity 10 , which are often associated with dysbiosis in the gut microbiota 32 . This observation, together with the highly distributed antibiotic consumption in COVID-19 patients (53.6% in our cohort, regardless of severity), warrants a reflection on the possibility that most of the patients could be dysbiotic at the time of illness. This particular scenario in which a previous dysbiotic microbiota caused by comorbidities and/or antibiotic consumption may affect susceptibility and outcome of COVID-19 has been suggested previously 20,37 . In other respiratory diseases such as COPD and asthma, it has been shown that dysbiosis in the respiratory microbiota can lead to a dysregulated immune response, increasing inflammatory processes 9,11,38 . Considering that aberrant immune responses are determinant in the progression of COVID-19, a previous dysbiotic respiratory microbiota could be affecting disease progression.
LefSe analysis (Fig. 3B) showed a differential abundance of microbial groups. For example, we found that most of the groups associated with the healthy controls belong to the so-called "normal" respiratory microbiota (e.g., Streptococcus, Oribacterium, and Veillonellaceae) 39 . By contrast, when we looked at the results of the microbiota associated with the COVID-19 and non-COVID-19 pneumonia groups, other potentially pathogenic microbial groups appeared. Specifically, in patients with non-COVID-19 pneumonia, we found bacteria associated with nosocomial infections such as Corynebacterium 40,41 . For mild COVID-19, we found some microbial groups associated with disease or bacteremia like Prevotella melaninogenica, V. parvula and Neisseria subflava 10,41 . For the case of severe COVID-19, we found Megasphaera which has been associated with the risk of ventilator-associated pneumonia (VAP) in other studies characterizing the microbiota of COVID-19 16 . Additionally, we found Rothia dentocariosa very abundant in deceased patients. This bacteria has been found as a causative agent of secondary www.nature.com/scientificreports/ pneumonia in H1N1 infection 33 and more recently has been associated with disease progression in previous studies characterizing the respiratory microbiota of COVID-19, proposing it as a biomarker for the disease 26,42 . From a clinical standpoint, it makes sense that a higher mortality predictor such as APACHE score correlates with poor survival in patients with COVID-19. Other clinical factors such as BUN or urea have also been used as markers of severity in respiratory diseases such as community-acquired pneumonia 43 . Acknowledging that there are multiple pathophysiological considerations still unexplained in SARS-CoV-2 infection, and the multisystemic involvement that has been observed in COVID-19 44 , biochemical markers of organ dysfunction such as lymphopenia, elevated myoglobin, and serum troponin levels, such as those found in this study, can help predict mortality in these patients 45 . Although the association of these factors with specific microbial groups in the respiratory tract has not been previously reported, the findings of this work pave the way for further study of the relationship between respiratory microbiota and clinical outcomes. The identification of pathogenic bacteria such as Actinomyces, Prevotella, and Corynebacterium in association with two or more clinical factors further supports the current line of research attempting to correlate the gut-lung axis with pulmonary disease 46 .
Recent studies, as well as this work, suggest that particularly anaerobic bacteria inhabiting the respiratory tract may be involved in the pathogenesis of COVID-19 and the host immune system. For example, in patients with high levels of blood urea, we found Corynebacterium, which has been previously linked to hyperuricemia 47 , and has also been found as an opportunistic pathogen in patients with immunosuppression or severe disease 48 . In addition, Prevotella has also been found to be increased in studies with patients with severe disease and has been linked to cardiac injury and high-risk mortality 45,49 . In this work, we found this specific genus associated with four clinical variables predictive of mortality in patients with COVID-19 (Fig. 4A,B). This finding is of particular interest considering previous evidence that Prevotella potentiates a Th17-mediated response through IL-8, CCL20, and IL-6 secretion 49,50 ; both the Th17 response and its cytokines are currently associated with the host immune response to SARS-CoV-2 51 .
Finally, the co-occurrence arrangement of ecological networks lets us identify structural patterns that reflect variations in the biological properties of COVID-19-associated microbial communities. For example, we found that all networks are distinguishable in terms of topological metrics such as density, clustering, and heterogeneity. It is worth mentioning that such metrics are potentially related to the stability of the systems as well as to other ecological properties such as resilience and redundancy 52 . In particular, we found a striking pattern of reduction in structural complexity from mild to fatal COVID-19. The loss of complexity manifests itself in a reduction in the number of nodes, edges (connections), density, and clustering, and moves from a highly connected and dense network (mild COVID-19) to a highly disaggregated and disconnected network (fatal COVID-19) (Fig. 5A,B, Supplementary Table S1).
These structural changes can lead to the generation of hypotheses about the consequences at the microbial community level. For example, changes in structural patterns could potentially be reflected in alterations in the ecological relationships between microorganisms. A common feature in MADs is that commensal/neutral bacteria can become pathogenic in the face of disease 32 . That is the case of bacteria such as Prevotella, Veillonella, Streptococcus, Actinomyces, or Megasphaera, which have been found as opportunistic pathogens in other COVID-19 microbiota characterization studies 16,[30][31][32] and were also found in this study [severe and fatal associated microbiota (Fig. 3B)]. The shift from neutral to deleterious interactions in specific bacteria could be the result of a loss of interactions that maintain the function and stability of microbial systems. This, in turn, could lead to exacerbated growth of potentially pathogenic microbial groups, but also the depletion of beneficial bacteria, altering the entire environment and possibly compromising the functions provided by the microbiota to the host.

Conclusions
Overall, this work provides insight into the role of the respiratory microbiota in COVID-19 disease. Although experimental validation is needed, our data suggest that environmental and host-related factors could be affecting the respiratory microbiota prior to SARS-CoV-2 infection, potentially compromising the immunological response of the host against disease and promoting secondary bacterial infections. We hypothesize that high levels of dysbiosis and poor microbial structural complexity in the respiratory microbiota of COVID-19 patients could be the result of antibiotic intake and comorbidities. Although further validation of our results is needed, an earlier dysbiotic state as that found in our study may have consequences at the host and microbial community level. On the one hand, increased dysbiosis in diseased patients could be modifying the PAMPs that stimulate a homeostatic antiviral response, allowing better conditions for SAR-CoV-2 replication. Additionally, the loss of structural complexity may lead to the emergence of opportunistic pathogens that, through ecological competition, may cause depletion of beneficial bacteria and promote secondary bacterial infections that worsen the clinical outcome. In summary, the findings of this work contribute to understanding the pathology of COVID-19 by identifying the respiratory microbiota as a potential factor affecting disease outcomes. It is worth mentioning that the main limitation of our study is the number of healthy subjects enrolled, which were all the eligible individuals in our research center at that time. Further research involving a larger number of healthy subjects is needed to confirm these findings. In addition, studies looking for the specific mechanisms by which dysbiotic respiratory tract microbiota compromise immune responses against virus infections are needed. Finally, further investigations using lower respiratory samples are needed to disentangle the behaviour of microbial communities in the lungs of patients with COVID-19.

Methods
Ethics statement. The Science, Biosecurity, and Bioethics Committee of the Instituto Nacional de Enfermedades Respiratorias revised and approved the protocol (B-1020). Informed consent was obtained from all subjects or their legal guardians (next of kin). Additionally, the Institution requested informed consent for

Composition analyses.
To determine if the respiratory microbiota of patients with different severity levels of COVID-19 and controls differed in microbial composition, we constructed stacked barplots of relative abundanceat phylum and genus levels. The relative abundance of each taxon was calculated with phyloseq and plotted with ggplot2 R packages. The colors were randomly assigned using randomcoloR R package.
Alpha diversity. We calculated the Shannon-Wiener diversity index with the "microbiome" R package. To detect potential differences between groups we conducted a Wilcoxon rank-sum test in the "vegan" R package.
Beta diversity. We carried out a Principal Coordinates Analysis (PCoA) with weighted UniFrac distance at ASV level in the "phyloseq" R package. Potential differences in beta diversity were addressed with a Permutational Analysis of Variance (PERMANOVA) with 999 permutations performed with the "vegan" R package. Additionally, we tested for dispersion and stochasticity as a proxy of dysbiosis in microbial communities 68 . For this purpose, we calculated the Ružička similarity metric in the "CommEcol" R package and performed a Wilcoxon rank-sum test to detect potential statistical differences between healthy controls and diseased groups in their intra-treatment sample similarities. Dysbiosis was assumed when the similarities between the healthy microbiota samples were significantly higher than the similarities between the diseased microbiota samples 12 .
Finally, to detect differentially abundant taxa associated with severity levels and controls, we performed a Linear Discriminant Analysis (LDA) with effect size (LefSe) at the ASV level using the web-based tool MicrobiomeAnalyst 69 . Only taxa with an LDA score higher than 1.5 and a p-value < 0.01 obtained from the Kruskal-Wallis test were used. www.nature.com/scientificreports/ Clinical data analyses. To analyze the clinical data associated with our patient's cohort, we transform each clinical variable into a binomial category according to its data distribution. We used cut-points based on the 25 and 75 percentiles for each variable. For example, for a given variable, we classified all samples with values above or equal to the 75 percentile as "1", and all samples with values under the 75 percentile as "2". Subsequently, for all clinical variables (N = 65) we constructed Kaplan-Meier survival curves in SPSS Statistics (version 21) (Chicago, Illinois, USA) by using the hospitalization days as time variable, the mortality status (either deceased or alive) as a dependent variable, and the specific clinical qualitative variables as exposure variable. Only those curves statistically (p < 0.05) and biologically meaningful were retained for the subsequent analyses.
In addition, to determine if there were differentially abundant bacteria associated with the several risk factors for the clinical variables obtained from the Kaplan-Meier curves, we performed a second LefSe analysis. From this result, only taxa with an LDA score higher than 1.5 and a p-value < 0.01 according to the Kruskal-Wallis test were used.
Network structure inference. We inferred the network structure for the microbiota associated with patients with different severity levels of COVID-19. Network calculation was performed in the software CoNet (v1.1.1 beta) 52 by using read counts summarized at the ASV level. One network was constructed for each severity level (samples; mild COVID-19: 37, severe COVID-19: 27, and fatal COVID- 19: 19). Only co-occurrences statistically supported by the three tested methods (Pearson, Spearman, and Kendall) with a correlation > 0.85 and a p-value < 0.01 were established as edges in the graphs. Also, we applied a multi-test correction using the Benjamini-Hochberg procedure. Network visualization was performed in Cytoscape (v. 3.8.2) 70 .
To further characterize the structure we computed metrics of the topology of each network using the Net-workAnalyzer plug-in in Cytoscape 70 and visualized them with a spider chart constructed in fmsb R package.