Cohort description

Between December 7th, 2022, and January 25th, 2023, we collected BALF and blood samples from a total of 64 individuals who were intubated due to severe pneumonia and tested positive for SARS-CoV-2. The individuals were sampled from among eight medical centers in the Beijing area of China. The median age of the probands was 80.00 (70.00–85.00), and 48 (75.00%) were male. Among the 64 probands, a total of 45 died within 28 d of intubation (non-survivors) and the remaining 19 were discharged after recovery (survivors). Samples were also collected from an additional cohort of severe Hospital Acquired Pneumonia/Community-Acquired Pneumonia (HAP/CAP) patients from November to December 2019 (before the COVID-19 pandemic), consisting of 20 survivors and seven non-survivors; there were no significant differences in the age or gender distributions compared to the COVID-19 patient cohort (Supplementary Table. 1).

Clinical parameters associated with 28-d mortality

We first investigated anthropometric characteristics, clinical parameters, and cytokine profiles associated with clinical outcomes in invasive ventilated COVID-19 patients. We found no significant differences between survivors and non-survivors with respect to age distribution (Mann–Whitney U-test, one-tailed p = 0.49), age, or blood parameters (Table 1); however, the survivors had significantly lower levels of SOFA (11.00, 9.50–14.00 in non-survivors compared to 8.00, 7.00–11.00 in survivors; p < 0.01), lower APACHE II scores (26.00, 24.00–27.00 in non-survivors and 19.00, 17.00–21.00 in survivors; p < 0.01), and higher levels of PaO2/FiO2 (132.86, 86.65–231.00 in non-survivors and 188.00, 145.00–344.00 in survivors; p < 0.05) at the time of intubation (Table 1). Cytokine levels in the BALF and plasma indicated that the BALF of non-survivors had significantly higher levels of interleukin (IL)-17 (4.39, 3.53–6.42 in non-survivors and 3.53, 2.73–5.47 in survivors; one-tailed p < 0.05), IL-6 (254.58, 21.09–782.93 in non-survivors and 21.16, 10.75–305.34 in survivors; one-tailed p < 0.05), IL-10 (4.66, 2.25–9.83 in non-survivors and 2.17, 1.71–4.13 in survivors; one-tailed p < 0.01), and TNF-α (22.14, 11.83–58.82 in non-survivors and 9.24, 4.68–27.76 in survivors; one-tailed p < 0.01), and significantly lower levels of IL-5 (1.81, 1.20–3.76 in non-survivors and 2.86, 1.84–4.07 in survivors; one-tailed p < 0.05) (Table 1). We also investigated the percentage of co-morbidities in non-survivors and survivors, but found no significant differences with respect to diabetes, hypertension, or respiratory system diseases.

Table 1 Clinical Characteristics of the Study PopulationSARS-CoV-2 genome analysis

Using long reads generated from the Oxford Nanopore Technology (ONT) platform and short reads from Illumina sequencing (Fig. 1a), we assembled and examined the genomes of the SARS-CoV-2 strains infecting our cohort. A total of 27 complete or near-complete and five partial genomes were assembled from the 64 patients (Supplementary Table. 2). Phylogenetic analysis indicated that each genome belonged to one of the two SARS-CoV-2 lineages that dominated this wave of Omicron infections, BA.5 and BF.7 (Fig. 1b, c). We examined the distribution of the BA.5 and BF.7 lineages between survivors and non-survivors, and found no significant difference in the survival rate between those infected with the two lineages; two survivors and 10 non-survivors were infected with BA.5, whereas six survivors and 14 non-survivors were infected with BF.7 (p = 0.41, Spearman’s rank sum test) (Fig. 1d). Single nucleotide polymorphism (SNP) analysis revealed a total of 319 and 243 SNPs in the BA.5 and BF.7 lineages, respectively. A recent genomic analysis of Beijing Omicron strains incorporated 682 genomes assembled primarily from the general population of COVID-19 patients in Beijing during the same time period as our cohort collection; 131 of the SNPs we identified here were also reported in that analysis (Supplementary Table. 3). There was no significant enrichment of individual SNPs in our cohort of invasive ventilated COVID-19 patients compared to the percentages among the 682 genomes (p > 0.05, Fisher’s exact test for all SNPs), nor were there significant differences in particular SNPs between the survivor and non-survivor groups (p > 0.05, Fisher’s exact test for all SNPs). Thus, it was not likely that particular SARS-CoV-2 Omicron lineages or individual SNPs were responsible for a higher occurrence of invasively ventilated COVID-19 cases or increased mortality.

Fig. 1

Analysis of SARS-CoV-2 genomes in the BALF samples. a Coverage distribution of reads mapped to reference SARS-CoV-2 genome, in which we found 8 complete genomes and 43 > 50% completeness. COVID-19 patients (n = 63). b Overview of SNP variant types in the assembled genome, in total we found 34 insertions, 18 deletions, and 583 substitutions. COVID-19 patients (n = 51). c Phylogenetic relationships of the 51 genomes assembled in this study vs. 682 SARS-CoV-2 genomes sampled in Beijing general population from the same period and deposited in GISAID. 20 genomes belong to BA.5 and 24 to BF.7 lineages. COVID-19 patients (n = 44). d The correlation between SARS-CoV-2 genomes of BF.7 and BA.5 lineages and 28-day mortality (survivors or non-survivors) based on Spearman’s rank test, no significant differences were found between the 28-day mortality rate in the two strains. COVID-19 patients (n = 32)

Host transcriptomic signatures in survivors and non-survivors of COVID-19

Human transcripts accounted for the majority of the BALF metatranscriptomic data, averaging 56.47% of the reads. These transcripts were used to profile host gene expression in the lung tissue of invasive ventilated COVID-19 patients and to identify signatures specific to COVID-19 patients through comparison to data from severe HAP/CAP patients before the pandemic. First, in comparing gene expression in the invasive ventilated COVID-19 patients compared to the severe HAP/CAP patients, we found a total of 2,324 significantly up-regulated genes (|fold change [FC]| > 1.5, p < 0.05) and 33,668 down-regulated genes (|FC | > 1.5, p < 0.05) (Supplementary Fig. S1a). Pathway enrichment analysis of the differentially expressed genes suggested that the intrinsic apoptotic signaling and response to virus pathways were significantly elevated in COVID-19 compared to HAP/CAP patients, reflecting a distinct host response to SARS-CoV-2 infections (Supplementary Fig. S1b).

Next, we compared gene expression between invasive ventilated COVID-19 survivors and non-survivors. In survivors, there were 489 significantly up-regulated genes, with the most strongly up-regulated genes including EGFR, PTPN23, and PANX1. There were 187 significantly up-regulated genes in non-survivors, the most strongly up-regulated of which included CXCL10, CCL20, and CCL8 (Fig. 2a). Pathway enrichment analysis indicated that genes associated with viral processes were significantly up-regulated in survivors. Neutrophil chemotaxis and neutrophil migration could be observed to be significantly enriched in non-survivors (Fig. 2b). Further, we compared the intersection of external differential genes (compared to HAP/CAP) with internal differential genes (COVID-19 survivors compared to non-survivors), and we identified 268 and 97 significantly up-regulated genes in survivors and non-survivors in COVID-19 patients, respectively (Fig. 2c). These genes demonstrate both up- and down-regulated genes in COVID-19 survivors, but exclude the portion shared by HAP/CAP that is specific to COVID-19 patients. Among the overlapping differential KEGG pathways in COVID-19 patients vs HAP/CAP patients, and COVID-19 Survivors vs. Non-survivors. Cytokine-cytokine receptor interaction, Chemokine signaling pathway, and IL − 17 signaling pathway were significantly enriched in non-survivors (Fig. 2d), agreeing with previous reports on COVID-19 patients.12 Further examination of the intersection of internal differential genes in COVID-19 patients and internal differential genes in HAP/CAP patients, eight of which were shared with survivors and non-survivors of invasive ventilated COVID-19: IGHA1, CACNA1C-AS3 upregulated in survivors of both cohorts, and BPIFA1, RIT1, CCDC146, CFAP300, ENSG00000290744, WTAP upregulated in non-survivors of both cohorts (Supplementary Fig. S1c). The gene BPIFA1 encodes a product with antimicrobial activity13 and was significantly down-regulated in survivors of COVID-19 patients with bacterial co-infections, suggesting attenuated host responses as a result of co-infection. IGHA1, which is associated with an antibody immune response, was significantly up-regulated in survivors, suggesting an increased immune response in COVID-19 patients. These genes were significantly associated with clinical outcomes in both COVID-19 and HAP/CAP patients. In addition, COVID-19 non-survivors showed significant enrichment of genes associated with neutrophil chemotaxis and lymphocyte migration compared to the HAP/CAP non-survivors, indicating that these factors may specifically explain clinical outcomes in those infected with COVID-19 (Fig. 1d).

Fig. 2

Gene Differential Expression Analysis in BALF of COVID-19 Patients. a Volcano plot comparing gene expression between the survival and non-survival groups of COVID-19 patients. 489 genes were up-regulated and 187 genes were down-regulated in the survival group (abs(FC) > 1.5, p-value < 0.05), including red-highlighted genes that were functionally characterized or implicated in subsequent analyses. b Gene Ontology (GO) enrichment analysis (Biological Process) of the differentially expressed genes identified in (a), where the representative genes involved in each GO term are indicated in parentheses. c Volcano plot displaying genes with differential expression between COVID-19 patients and HAP/CAP patients, which are also differentially expressed in the survival group of COVID-19 patients. These genes are specific to COVID-19 survival patients and include 268 up-regulated and 97 down-regulated genes. d KEGG enrichment analysis of the differentially expressed genes identified in (c). HAP/CAP patients (n = 27), COVID-19 patients (n = 63), non-survivors of the invasive ventilated COVID-19 patients (n = 45), survivors of the invasive ventilated COVID-19 patients (n = 18)

The lower respiratory tract microbiome and COVID-19 outcome

We next analyzed the lung microbial community composition of each patient by mapping metatranscriptomic data to reference genomes of viral, bacterial, and fungal species. Comparative to negative controls, BALF samples of our cohorts exhibited significant differences in microbial compositions in each comparison, in addition to significantly higher concentrations of nucleic acid after extraction and library construction (Supplementary Fig. S2; Supplementary Table. 4). Of the fragments that could not be mapped to the human transcriptome (i.e., microbial reads), sequences from the virome accounted for an average of 34.82% of the total microbial reads; these were primarily from SARS-CoV-2 (averaging 75.97% of the viral reads) and Human betaherpesvirus (averaging 8.89% of the viral reads) (Fig. 3a). Bacterial reads accounted for 40.52% of all microbial reads, with the most abundant phyla being Proteobacteria, Firmicutes, Actinobacteriota, and Bacteroidota (averaging 39.49%, 38.31%, 12.60%, and 7.64% of the bacterial reads, respectively). At the genus level, the most abundant taxa included Acinetobacter, Pseudomonas, and Klebsiella (averaging 12.14%, 9.48%, and 1.76% of the bacterial reads, respectively), all of which belong to the phylum Proteobacteria and are pathogens that have previously been reported in co-infections with COVID-19 (Fig. 3b). Importantly, fungal reads were found in all of the invasive ventilated COVID-19 individuals and accounted for an average of 24.65% of the microbial reads, with most of these mapping to Candida species (an average of 51.11% of fungal reads) (Fig. 3c).

Fig. 3

Alternation of microbiome and virulence factors in COVID-19 and HAP/CAP patients. a Composition and relative abundance of the top 10 viruses in invasive ventilated COVID-19 patients. b Relative abundance of the top ten bacterial phylum levels in invasive ventilated COVID-19 patients. c Relative abundance of fungi belonging to the genus Candida in COVID-19 patients. d Box plot of differential virus between invasive ventilated COVID-19 and HAP/CAP patients, significance was derived from Wilcoxon tests. e Box plot of differential bacteria between invasive ventilated COVID-19 and HAP/CAP patients, significance was derived from Wilcoxon tests. f Box plot of significantly different fungi between invasive ventilated COVID-19 and HAP/CAP patients. g Box plot of significantly different bacterium between survivors and non-survivors of the invasive ventilated COVID-19 patients. h Box plot of significantly different fungi between survivors and non-survivors of the invasive ventilated COVID-19 patients. i Relative abundance of the top ten virulence factors in invasive ventilated COVID-19 patients. j Heatmap of significantly different virulence factors between survivors and non-survivors of invasive ventilated COVID-19 patients. k Heatmap of differential virulence factors between invasive ventilated COVID-19 and HAP/CAP patients. HAP/CAP patients (n = 27), COVID-19 patients (n = 63), non-survivors of the invasive ventilated COVID-19 patients (n = 45), survivors of the invasive ventilated COVID-19 patients (n = 18). One-tailed wilcoxon rank-sum test was used for all significance statistics, *p-value < 0.05, **p-value < 0.01, ***p-value < 0.001, ****p-value < 0.0001

Compared to the severe HAP/CAP patients from 2019, our cohort of invasive ventilated COVID-19 patients had significantly different microbial signatures. First, invasive ventilated COVID-19 patients showed significant enrichment of SARS-CoV-2 (99.74%, 0–100 in invasive ventilated COVID-19 patients compared to 0 in HAP/CAP patients; p = 8.5e-14) Rhinovirus A (0%, 0–25.09 compared to 0; p = 0.0011) and Alphaherpes virus (0%, 0–5.82 compare to 0; p = 0.00079) among the viral reads (Fig. 3d). In the bacteriome, invasive ventilated COVID-19 patients were significantly enriched in Halanaerobiaeota at the phylum level (0.006%, 0–2.70 compared to 0%, 0–0.0051 in HAP/CAP patients; p = 3.7e-07) and Pseudomonas at the genus level (4.85%, 0.0026–42.12 compared to 0.46%, 0–56.55 in HAP/CAP patients; p = 0.04) (Fig. 3e). Several fungal species were also significantly enriched in invasive ventilated COVID-19 patients, such as Candida glabrata (0%, 0–86.95 compared to 0%, 0–64.66 in HAP/CAP patients; p = 5.7e-03) and Candida parapsilosis (0.0049%, 0–95.07 compared to 0%, 0–0.44 in HAP/CAP patients; p = 6.6e-07) (Fig. 3f).

In contrast, we found no significant differences in the percentages of SARS-CoV-2 reads out of all viral reads between invasive ventilated COVID-19 survivors and non-survivors (99.75%, 0.05–99.96 and 99.52%, 0–100, respectively; p = 0.48) (Supplementary Fig. S3), or other viruses. However, among bacterial reads, there were significant differences in the percentages of Alloprevotella (0.0006%, 0–1.41 in survivors and 0.023%, 0–32.86 in non-survivors; p < 0.05), Caulobacter (0.0002%, 0–1.14 in survivors and 0.0046%, 0–5.12 in non-survivors; p = 0.026), Escherichia-Shigella (0.000018%, 0–1.17 in survivors and 0.0016%, 0–5.91 in non-survivors; p = 0.023) and Ralstonia (0.0011%, 0–0.37 in survivors and 0.064%, 0–1.58 in non-survivors; p = 0.024) (Fig. 3g). Among non-survivors, the fungal reads showed significant enrichment of Aspergillus sydowii (0.029%, 0–0.56 in survivors and 0.048%, 0–5.88 in non-survivors; p = 0.05) and Penicillium rubens Wisconsin (0.00019%, 0–3.51 in survivors and 0.042%, 0–23.09 in non-survivors; p = 0.04). These enriched species indicated likely synergistic effects between bacterial, fungal, and viral pathogens that contributed to worse outcomes (Fig. 3h). Sputum culture results for a number of patients in the COVID-19 cohort were indeed positive for Escherichia or Aspergillus (Supplementary table. 5). Among the KEGG pathways enriched in COVID-19 Survivors compared to HAP/CAP Survivors (Supplementary Fig. S1d) and COVID-19 Non-survivors compared to HAP/CAP Non-survivor (Supplementary Fig. S1e), Coronavirus disease -- COVID-19 pathways and Pathogenic Escherichia coli infection and Salmonella infection pathways were both enriched, indicating the importance of co-infection of bacteria, especially Escherichia/ Salmonella with SARS-CoV-2 is signature of the COVID-19 cohort non-survivors and play important role in shaping the clinical outcome. We also found significant differences between severe HAP/CAP survivors and non-survivors in the abundance of Cardiobacterium (0%, 0–3.39 in survivors and 0.00028%, 0–0.0032 in non-survivors; p < 0.05), Ligilactobacillus (0.0036%, 0–0.66 in survivors and 0.035%, 0.004–1.41 in non-survivors; p = 0.011), and Lactobacillus (0.0013%, 0–0.15 in survivors and 0.025%, 0.00094–3.59 in non-survivors; p = 0.0083) among the bacterial reads (Supplementary Fig. S4); in the mycobiome, Aspergillus aculeatus (0%, 0–1.94 in survivors and 0.023%, 0–0.88 in non-survivors; p = 0.0014) and Aspergillus flavus (0%, 0–0.6 in survivors and 0%, 0–11.08 in non-survivors; p = 0.0091) were significantly enriched in non-survivors (Supplementary Fig. S5).

We also examined the toxicity factors that may have contributed to clinical outcomes. Among invasive ventilated COVID-19 patients, the most common toxicity factors included VF0273 (flagella), which was present in 92% of patients, VF0084 (Xcp secretion system), which was present in 65% of patients, and VF0467 (acinetobactin), which was present in 60% of patients (Fig. 3i). Among invasive ventilated COVID-19 patients compared to severe HAP/CAP patients in 2019, there was significant enrichment (one-sided Wilcoxon test, p < 0.05) of factors involved in motility (VF0237 [flagella], VF0430 [flagella], and VF0473 [polar flagella]); immune modulation (VF0367 [LPS] and VF0309 [PDIM]); adherence (VF0525 [PfbA], VF0418 [Scm], VF0145 [CBPs], and VF0354 [EfaA]); nutritional/metabolic factors (VF0151 [PsaA]), and the effector delivery system (VF0084 [Xcp secretion system]) (Fig. 3k). Toxicity factors that were significantly enriched in non-survivors of COVID-19 cohort included VF0504 (AdeFGH efflux pump) (0 in survivors and 0%, 0–2.85 in non-survivors; p = 0.04), VF0944 (HSI-3) (0 in survivors and 0%, 0–4.27 in non-survivors; p = 0.04), VF0469 (phospholipase D) (0 in survivors and 0%, 0–2.91 in non-survivors; p = 0.03), VF0472 (PNAG) (0 in survivors and 0%, 0–11.09 in non-survivors; p = 0.03), and VF0571 (RcsAB) (0 in survivors and 0%, 0–10.38 in non-survivors; p = 0.04) (Fig. 3j). We found also significant positive correlation between cytokine levels in the BALF and blood with the common toxicity factors. In the BALF, the strongest positive correlation was between IL-1β and VF0504 (AdeFGH efflux pump) and VF0470 (Phospholipase C). In the blood, all the cytokine levels except IL-6 were positively correlated with six toxicity factors including VF0368 (BvrR-BvrS), VF0003 (Capsule), VF0560 (Capsule), VF1138 (Curli fibers), VF0228 (Enterobactin) and VF0521 (ESX-3) (Supplementary Fig. S6). There were also significant differences between severe HAP/CAP survivors and non-survivors, however, none overlapped with the differential virulence factors we found in the COVID-19 cohort (Supplementary Fig. S7).

Correlational analysis of lung microbiome, host transcriptome, and immune responses

Last, we probed potential interactions between the microbiome and invasively ventilated COVID-19 host responses by analyzing correlational networks between and among omics datasets (Fig. 4). The virome had the largest number of significant correlations with host gene expression (76 pairs with |r | > 0.5) and with cytokine expression (12 pairs with |r | > 0.5). The most highly correlated pairs were Chelonid alphaherpesvirus 5 and RPS15 (r = 0.66), Chelonid alphaherpesvirus 5 and SRRM2 (r = 0.72). Chelonid alphaherpesvirus 5 and IL-4 were significantly positively correlated (r = 0.60). Chelonid alphaherpesvirus 5 and Human betaherpesvirus 6 A were showed positive correlations with IL-2 (both r = 0.59). The abundance levels of specific bacteria were positively correlated with the differential expression of 14 genes and three cytokines. Notably, the strongest negative correlation (r = -0.70) was between the genus Delftia and PSMC3IP expression levels, and the strongest positive correlation was between the genus Exiguobacterium and RPS15 expression (r = 0.66). Three fungal taxa were correlated with differential host gene expression, including correlations between Sordaria macrosporak-hell and OR5BD1P, Y. lipolytica CLIB89W29 and MUC5B (r = -0.62 and r = 0.603, respectively) (Fig. 4a). Positive correlations between bacteria and cytokines included Exiguobacterium and IL-4 (r = 0.64), Catellicoccus and IL-2 (r = 0.61), and an uncultured Catellicoccus strain and IL-2 (r = 0.60) (Fig. 4b). Potential synergistic effects were also observed within the microbiome, such as a significant correlation between flagella (VF0273) and 23 bacterial genera; the strongest correlation was with Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium (r = 0.68). Notably, there were also strong positive correlations between flagella (VF0273) and the fungal strains Y. lipolytica CLIB122 and Y. lipolytica CLIB89W29 (r = 0.65 and r = 0.66, respectively) (Supplementary Fig. S8).

Fig. 4

Multi-omics correlation analysis of COVID-19 patients. a Heatmaps display the correlation analysis of viral, fungal, bacterial, and virulence factor, with differentially expressed genes in COVID-19 patients. Spearman correlation coefficients with adjusted p-values < 0.05 are marked in colored cells. b Heatmaps depicting correlation analysis between differentially expressed genes of viral and bacterial with cytokine levels in invasive ventilated COVID-19 patients. Positions with adjusted p-values < 0.05 are shown in color code, and the Spearman correlation coefficient numbers are within each grid. Only absolute correlation coefficients greater than 0.5 are displayed. COVID-19 patients (n = 63)