Host-Microbiome Associations in Saliva Predict COVID-19 Severity

Abstract

Established evidence indicates that oral microbiota plays a crucial role in modulating host immune responses to viral infection. Following Severe Acute Respiratory Syndrome Coronavirus 2 - SARS-CoV-2 - there are coordinated microbiome and inflammatory responses within the mucosal and systemic compartments that are unknown. The specific roles that the oral microbiota and inflammatory cytokines play in the pathogenesis of COVID-19 are yet to be explored. We evaluated the relationships between the salivary microbiome and host parameters in different groups of COVID-19 severity based on their Oxygen requirement. Saliva and blood samples (n = 80) were collected from COVID-19 and from non-infected individuals. We characterized the oral microbiomes using 16S ribosomal RNA gene sequencing and evaluated saliva and serum cytokines using Luminex multiplex analysis. Alpha diversity of the salivary microbial community was negatively associated with COVID-19 severity. Integrated cytokine evaluations of saliva and serum showed that the oral host response was distinct from the systemic response. The hierarchical classification of COVID-19 status and respiratory severity using multiple modalities separately (i.e., microbiome, salivary cytokines, and systemic cytokines) and simultaneously (i.e., multi-modal perturbation analyses) revealed that the microbiome perturbation analysis was the most informative for predicting COVID-19 status and severity, followed by the multi-modal. Our findings suggest that oral microbiome and salivary cytokines may be predictive of COVID-19 status and severity, whereas atypical local mucosal immune suppression and systemic hyperinflammation provide new cues to understand the pathogenesis in immunologically naïve populations.

Keywords: COVID-19; inflammatory cytokines; machine learning; network analysis; saliva microbiome.

Figure 1.. Analysis of salivary microbiome between the 3 groups, control (36), mild/moderate (blue) and severe (yellow).

A: Principal coordinates analysis (PCoA) of Jaccard dissimilarity demonstrating the three clusters formed by the microbial constituents. Ellipsoids show the deviation of spread of each group (p<0.05, ADONIS). B. Violin plots of observed features diversity between the three groups (p<0.05, Wilcoxon rank-sum test). C.Violin plots of Shannon diversity between the three groups (p<0.05, Wilcoxon rank-sum test). D. Principal coordinates analysis (PCoA) of Euclidean distances of PhILR distances, demonstrating the three clusters between the conditions.

Figure 2.. Salivary and blood cytokines.

Heatmaps of statistically significant (Kruskal-Wallis) cytokines in control, mild/moderate and severe measured in A. Blood, and B. Saliva. Correlogram of C. Controls. D. All COVID patients. E. Mild/moderate category. F. Severe category, between blood (rows) and salivary (columns) cytokines. Areas coloured based on the Spearman correlation Rho and any correlations not statistically significant (P>0.05) were removed and replaced with a white square.

Figure 3.. Hierarchical feature selection for identifying biomarkers predictive of COVID-19 severity.

A) Hierarchical classification of COVID severity and the sub-models with relative explained variance of each group derived from 1st principal component of PCA. Modality-specific biomarker overlap between sub-models of B) abundance-based paradigms and C) inferred-interaction paradigms. In the inferred interaction paradigm, features are edges in the SSPN which can be decomposed into nodes which is not the case for the abundance-based paradigms in B. Clairvoyance feature selection results for D) sub-model y1 and E) sub-model y2 color coordinated by modality. Each marker on the scatter plot represents a unique hyperparameter/feature set combination that yields a specific accuracy using 10-fold cross-validation repeated with 10 different random states.

Figure 4.. Aggregate network representations of fitted sub-models showcasing predictive capacity – Aggregate network for sub-model A) y1 and B) y2.

The edge weights can be interpreted as predictive capacity for COVID severity. For y1, positive values indicate that an increase in perturbation results in an increased likelihood that a sample is classified as y2 (i.e., infected) relative to uninfected. For y2, positive values indicate that an increase in perturbation results in an increased likelihood that a sample is classified as severe relative to mild/moderate. Node size is proportional to weighted degree as a measurement of network connectivity and, by extension, biomarker importance. C) ASVs in the SSPN as proportions of presence in the various conditions, demonstrating that the majority of the ASVs are due to differences in abundances of common taxa across all conditions, and not rare species.