Integrated transcriptomic and immune analysis reveals distinct mitochondrial and immune signature in COVID-19 ARDS requiring invasive mechanical ventilation

Abstract

Introduction:

Severe COVID-19–related acute respiratory distress syndrome (ARDS) often progresses to respiratory failure requiring invasive mechanical ventilation (IMV). Although several biomarkers (e.g., CRP, suPAR, Ang-2) have been explored to predict COVID-19 disease severity, validated early-presentation biomarkers that specifically prognosticate the need for IMV and mechanistic pathways in ARDS remain limited.

Research question:

What transcriptomic and immune signatures distinguish IMV from non-IMV in patients with COVID-19 ARDS, and how do these molecular changes contribute to disease severity?

Methods:

Patients (n = 36) with SARS-CoV-2 infection requiring oxygen support were enrolled in this prospective study. Blood samples were collected at baseline, day 4, and day 8 from patients on IMV (n = 11) and non-IMV (n = 25). Transcriptomic profiles of peripheral blood mononuclear cells were assessed using RNA-sequencing, followed by gene ontology and KEGG pathway enrichment. Plasma cytokines and chemokines were measured using multiplex immunoassays, peripheral immune cell endotypes were characterized by mass cytometry, and neutrophil extracellular trap (NET) formation was evaluated by kinetic live-cell imaging.

Results:

Clinically, IMV patients had significantly longer hospital (43.0 vs. 19.5 days, p = 0.001) and ICU stays (39.6 vs. 9.1 days, p < 0.01) and higher 28-day mortality (45.5% vs. 12.0%; aOR = 8.78, p = 0.04). Transcriptomic profiling identified 36 DEGs at day 4 and 21 at day 8 in patients on IMV, characterized by downregulation of mitochondrial gene MTARC2 and anticoagulant regulator TFPI, and upregulation of stress-response genes (RNF165, KCNMA1, and AC245014.3) (p < 0.05). Plasma IL-6, IL-8, and IL-10 were elevated in IMV group, with IL-10 showing the strongest predictive value (AUC 0.74, p = 0.0002). Immune profiling revealed higher percentage of CD57⁻ NK cells at baseline in non-IMV group. NET formation did not differ between groups.

Conclusion:

Patients requiring IMV demonstrated distinct transcriptomic and immune profiles associated with severe COVID-19. These findings identify candidate biomarkers and pathways associated with disease severity; however, the results are exploratory and require validation in larger cohorts and functional studies.

1 Introduction

Severe acute respiratory syndrome from SARS-CoV-2, a strain of RNA coronavirus that caused the COVID-19 pandemic, endangered the lives of millions globally (1). Among hospitalized patients, acute respiratory distress syndrome (ARDS) emerged as the most fatal complication, characterized by severe impairment of gas exchange due to intense pulmonary inflammation, resulting in alveolar–capillary leakage, surfactant depletion and dysfunction, and ultimately pulmonary fibrosis (2). This process is largely driven by a catastrophic hyperinflammatory response, termed the cytokine storm. The resulting inflammation is associated with profound oxidative stress and triggers expression of inducible nitric oxide synthase (iNOS), leading to excessive production of reactive oxygen species (ROS) and reactive nitrogen species (RNS). Together, these processes impair mitochondrial function and exacerbate intracellular hypoxia (3). Emerging evidence points to distinct mitochondrial alterations in COVID-19 ARDS, but comprehensive transcriptomic profiling to define disease progression to respiratory failure has not yet been performed (4).

Neutrophil extracellular traps (NETs) have also been implicated in the pathogenesis of severe COVID-19 through their proposed roles in inflammation, endothelial injury, immunothrombosis, and ARDS progression (5). Excessive NET formation has been associated with disease severity and adverse outcomes in COVID-19, supporting investigation of NETosis as a potential contributor to respiratory failure requiring invasive mechanical ventilation (IMV).

Clinically, approximately one in five patients with COVID-19 develops respiratory distress requiring supplemental oxygen or IMV (6), and among those admitted to the ICU, nearly one in three die (7). The ability to identify reliable prognostic biomarkers is therefore critical for recognizing patients at increased risk of respiratory decompensation and progression to IMV.

We hypothesized that mitochondrial dysfunction represents a defining feature of severe COVID-19 ARDS, and that cytokine release and immune exhaustion are interconnected consequences of impaired mitochondrial activity. To test this hypothesis, we evaluated transcriptomic signatures of peripheral blood mononuclear cells (PBMCs), plasma cytokine levels, and immune cell phenotypes in patients with COVID-19, comparing those who required IMV with those who did not.

2 Materials and methods

This prospective observational study was conducted between September 10, 2020, and March 1, 2022, at Houston Methodist Hospital. The study was approved by the Institutional Review Board, and informed consent was obtained from all participants or their legal representatives. Adults (≥18 years) with PCR-confirmed SARS-CoV-2 infection requiring supplemental oxygen were enrolled and stratified according to the WHO clinical progression scale into non-mechanically ventilated (scores 5–6) and mechanically ventilated (scores 7–9) groups.

Fresh whole blood was collected in EDTA tubes at baseline, Day 4, and Day 8, with a maximum total collection volume of 180 mL. PBMCs were isolated within 24 h of collection and stored at −150 °C for subsequent RNA sequencing and mass cytometry (8). Plasma samples were analyzed for cytokine and chemokine levels and used to stimulate healthy donor neutrophils for quantification and visualization of NETosis. All patient data and biospecimens were de-identified, and researchers performing the analyses were blinded to patient group assignments. Survivors were contacted by telephone to assess 28-day mortality outcome.

Sample Size and Power Calculation: Initially, we calculated that a minimum of 20 patients per group (IMV vs. non-IMV) would provide 80% power to detect a 2-fold or greater difference between groups at a significance level of 0.05. During the study period, we enrolled a total of 36 patients, of whom 11 (30%) required IMV and 25 did not. The overall number of patients requiring IMV was lower than initially projected, likely reflecting declining intubation rates with the waning of the COVID-19 Delta variant, widespread vaccination, and advances in supportive care. Nevertheless, because 30% of our cohort (11/36) required IMV, we were able to proceed with comparative analyses between the two groups.

Statistical Methods: We reported means and standard deviations (SD) for continuous variables, and proportions for categorical variables. Histograms and the Shapiro–Wilk test were used to assess the normality of continuous variables and outcomes. Student’s t-test was used to assess mean differences between continuous variables. The chi-squared test (χ²) or Fisher’s exact test, as appropriate, was used to evaluate differences between categorical variables. The primary variable of interest was the IMV status: patients who were on IMV (1) vs. those who were not (0). The primary outcomes included 28-day mortality (1 = Died, 0 = Alive), 30-day readmission (1 = Readmitted, 0 = Not readmitted), hospital length of stay, and ICU length of stay. For binary outcomes, both univariate and multivariate generalized logistic regression models were used. For continuous outcomes, such as hospital and ICU length of stay, generalized linear regression models were employed. All models were adjusted for demographics (age, gender, and BMI), comorbidities, Charlson Comorbidity Index (CCI), and medications to reduce bias. A p-value <0.05 was considered statistically significant. We reported adjusted odds ratios (aOR) with 95% confidence intervals (CI), as well as coefficient estimates (Δβ). All statistical analyses were performed using R version 4.5.1.

2.1 RNAseq and gene expression analysis

PBMCs were isolated from blood samples collected at baseline, day 4, day 8, and, where applicable, on the day of IMV initiation and liberation. Frozen cell pellets were then submitted to Azenta Life Sciences/GENEWIZ (South Plainfield, NJ) for Whole-transcriptome RNA sequencing (Tracking #30-696137733). RNA extraction from frozen cell pellets, library preparation using poly-A selection for mRNA enrichment, and sequencing on an Illumina platform with paired-end 2 × 150 bp configuration targeting 20–30 million reads per sample were performed by GENEWIZ. RNA integrity and library quality were assessed by GENEWIZ using Agilent TapeStation, Qubit fluorometric quantification, and qPCR before sequencing. Samples with inadequate RNA integrity or sequencing quality metrics were excluded from downstream analyses. A total of 95 samples from 11 IMV and 25 non-IMV patients across available timepoints were included in the final analysis. A comparison of gene expressions between the groups of samples was performed using DESeq2. The Wald test was applied to generate p-values and log2 fold changes. Genes with an adjusted p-value <0.05 and absolute log2 fold change >1 were classified as differentially expressed genes (DEGs). Gene ontology (GO) analysis was performed using g:Profiler (version e107_eg54_p17_bf42210, database updated 15/09/2022) to visualize the enrichment information of the significant DEGs. Homo sapiens was chosen as a specific organism for analysis. GO molecular function (GO:MF), GO cellular component (GO:CC), and GO biological process (GO:BP) categories were used as Gene Ontology data sources. Pathway analysis used KEGG, Reactome, and WikiPathways databases. Analyses were restricted to annotated genes with significance thresholds set using g:SCS multiple testing correction method at α = 0.05. DEGs were further analyzed using gene set enrichment analysis via ShinyGO (9). Significant DEGs were imported into the ShinyGO V0.77 tool to perform KEGG pathway analysis and obtain a visual graph of KEGG pathways. When submitting the list of gene targets in the Gene Official Symbol format, the selected matching species were human, and the FDR cutoff was set to 0.05 (p ≤ 0.05). After running the analysis, the top five KEGG pathways that matched the conditions were sorted by number of genes and visualized in the graphs.

The GOnet web application was used to identify biological functions of significant DEGs at Day 4 (10). The following analysis parameters were used: species: homo sapiens; input proteins: 36; GO namespace: biological process; analysis type: GO-term enrichment; q-value threshold of 0.05. Enrichment analysis of significant DEGs at Day 4 was identified in DisGeNET, a database of gene-disease associations (11). All genes in the genome were used as the enrichment background. Terms with p < 0.01, a minimum count of 3, and an enrichment factor > 1.5 (the enrichment factor is the ratio between the observed counts and the counts expected by chance) were collected and grouped into clusters based on their membership similarities. The overlap and unique significant DEGs from Day 4 and Day 8 in the IMV group versus the non-IMV group were analyzed as Venn diagram using jvenn (12).

2.2 Cytokine and chemokine assay

Plasma concentrations of 38 analytes were measured by a Luminex multiplex immunoassay (Human Cytokine/Chemokine Magnetic Bead Panel, HCYTMAG-60 K-PX41, EMD Millipore, Burlington, MA). The panel included EGF, eotaxin, FGF-2, Flt-3L, fractalkine, G-CSF, GM-CSF, GRO, IFN-α2, IFN-γ, IL-1α, IL-1β, IL-1RA, IL-2, IL-3, IL-4, IL-5, IL-6, IL-7, IL-8, IL-9, IL-10, IL-12p40, IL-12p70, IL-13, IL-15, IL-17A, IP-10, MCP-1, MCP-3, MDC, sCD40L, MIP-1α, MIP-1β, TGF-α, TNF-α, TNF-β, and VEGF. Baseline, Day 4, Day 8, and/or Day-of-discharge samples were included in the longitudinal analysis, comprising 11 IMV and 25 non-IMV patients. For patients discharged before the scheduled Day 4 or Day 8 collection time points, the sample obtained closest to the day of discharge was assigned to the corresponding Day 4 or Day 8 time point for longitudinal analysis. Healthy donor controls were analyzed in parallel for cytokine/chemokine profiles. Outliers were identified with Grubbs’ test and excluded. Analyte concentrations were assessed using mixed-effects models with patient group (IMV vs. non-IMV) and time point (baseline, Day 4, Day 8) as fixed factors, without assuming sphericity. Significant interaction effects were further evaluated with Tukey’s post hoc tests. Diagnostic performance of differentially expressed analytes was assessed by receiver operating characteristic (ROC) curve analysis, with area under the curve (AUC) used to estimate discrimination between IMV and non-IMV patients. Statistical analyses were performed using GraphPad Prism (version 10.0.3; GraphPad Software, Boston, MA), with significance set at α = 0.05.

2.3 NETosis assay

Peripheral neutrophil extracellular trap (NET) formation was quantified in patient’s plasma by measuring myeloperoxidase (MPO)–DNA complexes using a commercial ELISA kit (11774425001; Roche, Mannheim, Germany). Age-matched samples from 10 IMV and 10 non-IMV patients were included in the analysis. To further evaluate NETosis, neutrophils isolated from healthy donors were incubated with plasma from either healthy controls or patients with COVID-19 at different disease stages. NET release was assessed by time-lapse imaging using Incucyte^® Netosis Assay (Sartorius, Göttingen, Germany) (13).

2.4 Mass cytometry

Cytometry by time-of-flight (CyTOF^®, Fluidigm) was performed for comprehensive immune profiling using patients’ PBMCs at an institutional core facility. A total of 30 PBMC samples from 10 patients (6 IMV and 4 non-IMV patients), collected at multiple time points, were analyzed using two CyTOF panels: a 36-marker myeloid panel (M-panel) and a 36-marker T cell panel (T-panel). Prior to staining, cryopreserved PBMCs were thawed and assessed for viability. Standardized quality control procedures were applied, including normalization using calibration beads, exclusion of debris and dead cells, and doublet discrimination before downstream analyses. Samples with inadequate viability or poor event quality were excluded. Data preprocessing and gating strategies were performed uniformly across all samples to minimize batch-related variability and ensure consistency in immune cell population identification. We analyzed differential abundance of immune cell populations between IMV and non-IMV groups at each time point using unpaired t-test or Welch’s t-test with Holm-Šídák correction for multiple comparisons. Within-group differences were evaluated using one-way mixed-effects models without assuming sphericity, followed by Tukey’s post hoc tests. Statistical analyses were performed using GraphPad Prism, with significance set at α = 0.05.

To address multiple testing concerns, corrections for multiple comparisons were applied separately for the transcriptomic, cytokine, and immunophenotyping analyses. In the RNA-seq analysis, differential gene expression was assessed using DESeq2 with Benjamini–Hochberg false discovery rate (FDR) adjustment, and genes with an adjusted p-value <0.05 were considered significant. For pathway enrichment analyses, significance thresholds were controlled using g:SCS correction in g:Profiler and FDR-adjusted p-values in ShinyGO. In the cytokine and chemokine analyses, mixed-effects models were followed by Tukey’s post hoc tests to account for multiple pairwise comparisons across groups and time points, and adjusted p-values are reported where applicable. For CyTOF immunophenotyping, differential immune cell abundance between IMV and non-IMV groups was evaluated using unpaired t-tests or Welch’s t-tests with Holm–Šídák correction for multiple comparisons, while longitudinal within-group analyses used mixed-effects models followed by Tukey’s post hoc testing. Adjusted p-values were used to determine statistical significance throughout these analyses.

3 Results

Our study population consisted of 36 patients, of whom 11 required IMV. The mean age was similar between the two groups (65 vs. 62; p = 0.49). 64% percent of the patients were male, with a mean BMI of 31 in the intubation group. Additionally, 18% of patients in the intubation group were ex-smokers. The mean daily steroid equivalent dosage was 9 mg in the IMV patients versus 34 mg in the non-IMV patients (p = 0.01). Table 1 describes the baseline characteristics of the study population.

45.5% of patients in the IMV group died after the 28th day, compared to 12.0% of patients who were in the non-IMV group. Additionally, the multivariate model showed that IMV group patients were eight times more likely to die by the 28th day compared to non-IMV group patients (aOR = 8.78; 95% CI [1.09–70.85]; p = 0.04). Table 2 presents the full model.

9% of patients in the IMV group were readmitted within 30 days, compared to 12% of patients in the non-IMV group. The multivariate model showed that IMV group patients were twice as likely to be readmitted within 30 days; however, this difference was not statistically significant (aOR = 2.77; 95% CI [0.09–86.38]; p = 0.56).

The mean length of hospital stay (in days) was higher in patients in IMV group compared to those in the non-IMV group (43.0 vs. 19.5; Δβ = 19.67; 95% CI [6.8–32.6]; p = 0.01). Table 3 presents the full model. The mean ICU length of stay (in days) was also higher in patients in IMV group compared to those in the non-IMV group (39.6 vs. 9.1; Δβ = 22.4; 95% CI [9.87–34.99]; p < 0.001). Table 4 presents the full model.

3.1 Transcriptomic profiling and analysis of PBMCs

To identify underlying molecular mechanisms of disease severity in patients with COVID-19, we isolated PBMCs and performed RNAseq and gene expression analysis, comparing transcriptomic profiles of PBMCs from the IMV group with those from the non-IMV group. Global transcriptional changes are visualized in Figure 1, where each data point in the volcano plot represents a gene, with log2-fold change on the x-axis and the -log10 adjusted p-value on the y-axis. Genes with adjusted p < 0.05 and log2-fold >1 (red) were considered upregulated, while those with adjusted p < 0.05 and log2-fold change <−1 (blue) were considered downregulated. When comparing patients in IMV with the non-IMV group across all time points, a total of 0, 36, and 21 genes were differentially expressed in the IMV group compared with the non-IMV group at baseline, Day 4, and Day 8, respectively. These correspond to 0.00, 0.20, and 0.11% of all genes assessed (18,222, 18,222, and 18,365 genes, respectively) (Figure 1). These DEGs comprised 0, 11, and 16 upregulated genes (accounting for 0.00, 30.56, and 76.19% of all significant DEGs) and 0, 25, and 5 downregulated genes (accounting for 0.00, 69.44, and 23.81% of all significant DEGs) at baseline, Day 4, and Day 8, respectively.

Figure 1

To explore the biological significance of these changes, we next performed pathway enrichment analysis. At Day 4, gene ontology molecular function (GO:MF) analysis identified enrichment of a single term (vascular endothelial growth factor receptor activity). In contrast, 16 terms were enriched in gene ontology biological processes (GO:BP). Many related to cellular response to stimuli and regulation of multicellular processes, including negative regulation of response to external stimuli, negative regulation of chronic inflammatory response, and response to carbon monoxide (Supplementary Figure 1a). KEGG pathway analysis further demonstrated enrichment of Renin secretion and Rap1 signaling pathways in PBMCs from IMV group compared to non-IMV (Supplementary Figure 1b). DisGeNET enrichment analysis linked differentially expressed genes to clinical phenotypes such as anoxia, essential thrombocythemia, neointima formation, visceral pain, and bronchopulmonary dysplasia (Supplementary Figure 1c).

We identified biological functions of significant DEGs at Day 4 to evaluate gene-term interactions. As shown in Figure 2a, the negative regulation of response to external stimulus pathway involved the downregulation of TFPI, PDGFRA, TNFAIP6, and KREMEN1 expression. Negative regulation of the chronic inflammatory response pathway included the downregulation of IL-10 and CYP19A1 expression, whereas response to carbon monoxide involved the up-regulation of KCNMA1 and the downregulation of IL-10 expression. We further identified uniquely significant DEGs from Day 4 and Day 8 of patients on IMV compared with non-IMV group as shown by Venn diagram (Figure 2b). Five DEGs were shared across both time points: three were upregulated (RNF165, KCNMA1, and AC245014.3) and two were downregulated (TFPI and MTARC2) in IMV patients compared with non-IMV patients (Table 5).

Figure 2

3.2 Immunoassays of peripheral cytokines

Because transcriptomics highlighted dysregulation of IL-10 and other inflammatory mediators, we next evaluated whether these changes were mirrored at the protein level. Plasma cytokine analysis revealed significantly higher levels of IL-10, IL-8, and IL-6 in patients on IMV compared with non-IMV (Figure 3a). A significant interaction effect between patient group and time point was observed for IL-15, MCP-1, eotaxin, GM-CSF, fractalkine, and MIP-1α (Figure 3b). Post hoc analyses revealed that the IMV group had significantly elevated IL-15 at Day 8 and higher MCP-1 at Day 4 compared to the non-IMV group. To assess diagnostic potential, receiver operating characteristic (ROC) analyses were performed for IL-10, IL-8, and IL-6 (Figure 4). All three cytokines demonstrated statistically significant AUC values, with IL-10 showing the strongest discrimination (AUC = 0.74, 95% CI 0.63–0.86, p = 0.0002), followed by IL-8 and IL-6. Additional cytokines such as IP-10, MCP-1 and MIP-1β also showed discriminatory potential, though less robust than IL-10 (Table 6). Collectively, these results indicate that patients on IMV had higher IL-10 levels than 75% of the non-IMV patients, highlighting IL-10 as the most promising biomarker. Flt-3L, IL-17A, and TGF-α were excluded from analysis because more than 50% of values were below assay detection limits.

Figure 3

Figure 4

3.3 NETosis analysis

We investigated whether the cytokine elevations resulting in systemic inflammation translated into enhanced neutrophil activation. However, analysis of patient plasma revealed no significant differences in MPO–DNA complex formation between IMV and non-IMV groups. In addition, NETosis was not detected in neutrophils incubated with plasma from either healthy controls or COVID-19 patients at different disease stages (Supplementary Figures 2a,b).

3.4 Immune cell phenotype analysis

We next performed immune cell phenotype analysis to evaluate differences in immune cell phenotypes by analyzing the proportion of CD45 + cells across immune cell subtypes. At baseline, non-IMV group had significantly higher percent of CD57⁻ natural killer (NK) cells compared with IMV group (Figures 5a,b). Since CD57 expression is a marker of NK cell senescence and reduced proliferative capacity, these findings suggest that patients on IMV may harbor a more senescent NK cell compartment, potentially impairing effective antiviral immunity and contributing to worse outcomes.

Figure 5

4 Discussion

In this study, we identified transcriptomic and immunological determinants of COVID-19 severity by comparing patients who required IMV with those managed without IMV. Our findings suggest the presence of altered mitochondrial signatures, dysregulated coagulation, cytokine imbalance, and immune exhaustion that may collectively drive progression to severe ARDS requiring invasive mechanical ventilation, consistent with emerging evidence implicating convergent immunometabolic dysfunction as a defining feature of critical COVID-19 (14–16). However, these observations should be interpreted cautiously, as the analyses were exploratory in nature, based on PBMCs, and do not establish causality or tissue-specific mechanisms underlying respiratory failure.

Transcriptomic analysis revealed a dynamic signature of disease progression. The number of differentially expressed genes was minimal at baseline, peaked on Day 4, and declined by Day 8, suggesting that systemic transcriptional dysregulation intensifies during early clinical deterioration. Profound downregulation of MTARC2, a gene that encodes a molybdenum enzyme, involved in drug metabolism and N-reductive activity, and its disruption impairs lipid and energy metabolism, suggesting a role in cellular energy regulation (17–19). The mechanism underlying reduced MTARC2 expression in PBMCs of IMV group remains unclear; however, one possibility is that impaired energy regulation, reflecting cellular hypoxemia particularly in respiratory muscles, may limit the ability to sustain spontaneous ventilation and require invasive support. Importantly, in a prior study showed that in COVID-19, the virus suppresses key mitochondrial genes needed for energy production and keeps them suppressed in critical organs like the heart, kidney, and liver, even after the virus clearance, highlighting a SARS-CoV-2-specific pattern of mitochondrial dysregulation (14). Concurrently, stress-response and ion-channel genes (RNF165, KCNMA1, AC245014.3) were upregulated, and may indicate maladaptive signaling networks in mitochondrial homeostasis and immune regulation. Of particular relevance, KCNMA1 encodes mitochondrial BK channels that regulate ROS handling under hypoxic stress, while RNF165, a RING-domain E3 ubiquitin ligase, has been shown to enhance BMP-Smad1/5/8 signaling in motor neurons through ubiquitin-mediated degradation of pathway inhibitors. Disruption of BK channel function may disturb the firing patterns of brainstem respiratory centers (such as the pre-Bötzinger complex, the key rhythm generator located in the medulla), impairing central respiratory rhythm generation and potentially contributing to respiratory failure (20). RNF165 was (16, 21) among the most consistently upregulated genes in IMV patients at both Day 4 and Day 8. Notably, RNF165-null mice exhibit shortening of phrenic nerve presynaptic branches within the diaphragm alongside evidence of impaired breathing, suggesting that dysregulated RNF165-mediated BMP signaling in phrenic motor neurons may compromise neuromuscular transmission to the diaphragm and contribute to ventilatory pump failure requiring invasive mechanical ventilation (16). The upregulation of AC245014.3, though of unknown function, may represent an unrecognized host–virus interaction. These findings position mitochondria not as passive injury targets but as active mediators of immunometabolic dysfunction.

Among the most significantly downregulated genes was TFPI, a regulator of the TF/FVIIa/FXa complex, limiting thrombin generation and exerting anti-inflammatory effects (22). Loss of TFPI has been linked to coagulopathy and mortality in COVID-19, consistent with our findings and with reports by Al-Tamimi et al. (23). This finding contrasts with systemic TFPI elevations in plasma reported during COVID-19 and may suggest compartment-specific transcriptional exhaustion under sustained inflammation (24). This paradox underscores a mitochondria–coagulation axis in which bioenergetic failure compromises vascular homeostasis and predisposes to thrombosis.

Pathway enrichment analyses further identified vascular injury and inflammatory signaling. KEGG identified enrichment of renin secretion and Rap1 signaling in patients on IMV, consistent with prior transcriptomic studies implicating dysregulation of the renin–angiotensin system and altered Rap1 signaling in severe COVID-19 (25, 26). These pathways may contribute to endothelial dysfunction, inflammation, and thrombosis associated with ARDS progression, however, these associations remain exploratory and should be validated in larger studies incorporating tissue-specific and functional analyses.

Cytokine profiling revealed elevated IL-10, IL-6, and IL-8 in patients on IMV, with IL-10 emerging as the strongest prognostic marker (AUC 0.74). While classically anti-inflammatory, IL-10 is significantly increased in patients with COVID-19 ARDS compared to those without ARDS, and this elevation is linked to disease severity, risk of pulmonary fibrosis, and adverse outcomes (27, 28). This paradoxical elevation reflects a failing compensatory mechanism or monocyte resistance to IL-10 with downstream T- and NK-cell hyperactivation. Taken together, IL-10 serves both as a marker of immune dysregulation and as a prognostic biomarker of disease severity.

A notable discrepancy emerged with IL-10 expression. Multiplex immunoassay revealed higher plasma IL-10 levels in IMV group, yet PBMC RNA-seq demonstrated downregulation of IL-10 transcripts. This likely reflects cellular composition, as IL-10 is predominantly secreted by monocytes and macrophages, which are underrepresented in PBMC isolates. Plasma IL-10 thus represents systemic monocyte/macrophage activation, while PBMC transcriptional profiles primarily capture T- and B-cell compartments.

Interestingly, we did not observe significant differences in NET formation between IMV and non-IMV groups. This contrasts with prior studies that associated excessive NET release with severe COVID-19, ARDS, and thromboinflammation (29), such as those by Melero et al. (30). Potential explanations include differences in patient populations, timing of sample acquisition, clinical management (e.g., anticoagulation use), and methodological limitations in NET quantification. While our results diverge from prior reports, they should be interpreted cautiously given the modest sample size, and they do not diminish the established role of NETs in COVID-19 pathophysiology.

Our immune phenotyping results extend these observations. Patients in the non-IMV group had a higher percentage of CD57-negative NK cells at baseline, representing a less terminally differentiated, more proliferative subset capable of expansion and functional response during viral infection. Patients requiring IMV, by contrast, exhibited enrichment of CD57-positive NK cells, consistent with a terminally differentiated phenotype characterized by reduced proliferative capacity and altered cytokine responsiveness, which may reflect impaired capacity for adaptive NK cell expansion against SARS-CoV-2. These findings align with studies by Zheng et al. and Varchetta et al., which described NK cell depletion, functional exhaustion, and skewing toward terminally differentiated subsets in severe COVID-19 (31, 32). Although these findings suggest that baseline NK cell differentiation status may be associated with disease severity and progression to respiratory failure requiring IMV, larger prospective studies are needed to determine their prognostic value and potential therapeutic implications.

4.1 Limitations

This study has a few limitations. First, there was a significant difference in corticosteroid exposure between the IMV and non-IMV groups, with the non-IMV cohort receiving higher daily steroid doses. Because corticosteroids are known to exert broad immunomodulatory effects, including suppression of inflammatory signaling pathways and alteration of gene transcription profiles, it is possible that some of the observed immune exhaustion signatures and transcriptomic differences were influenced, at least in part, by differences in steroid exposure rather than exclusively reflecting disease severity. Importantly, however, patients requiring IMV demonstrated evidence of immune dysregulation despite receiving lower mean steroid doses, suggesting that the observed alterations are unlikely to be solely attributable to corticosteroid therapy. Nevertheless, we recognize steroid use as a potential confounding factor that may have affected both immune phenotypes and gene expression patterns. Given the limited sample size, formal multivariable adjustment was constrained; however, this issue should be interpreted carefully when considering the biological significance of our findings. Future studies with larger cohorts and standardized corticosteroid protocols will be important to better distinguish the independent contributions of disease severity and steroid exposure to the immune and transcriptomic changes observed in critically ill patients. Second, A limitation of this study is the relatively small sample size across the different analytical platforms, including RNA sequencing, CyTOF, and cytokine profiling. The number of samples analyzed for each technique is reported in the corresponding methods. Although the longitudinal and multi-omics design provided important mechanistic insights, the limited cohort size may reduce statistical power, increase susceptibility to type I and type II errors, and raise the possibility of overfitting in multivariable analyses. Therefore, the findings should be interpreted cautiously and considered hypothesis-generating. Larger, independent validation cohorts will be necessary to confirm the robustness and generalizability of the observed immune and transcriptomic signatures. A major limitation of this study is that the transcriptomic and immunophenotyping analyses were performed using peripheral blood mononuclear cells (PBMCs) rather than tissue-derived samples from the lung, respiratory musculature, or central nervous system. As a result, the observed molecular signatures reflect systemic peripheral immune responses and may not fully represent tissue-specific mechanisms driving respiratory failure or ARDS progression in COVID-19. Future studies integrating tissue-based analyses, single-cell approaches, and functional validation experiments will be necessary to define the mechanistic relevance of these pathways in severe COVID-19 and the need for invasive mechanical ventilation.

4.2 Interpretation

In summary, this study identifies mitochondrial and immune dysfunction as defining molecular features of severe COVID-19–associated ARDS. First, mitochondrial gene expression patterns, along with plasma IL-10 levels and NK cell phenotypes, may serve as biomarkers to stratify patients at risk of respiratory failure and IMV. Second, therapeutic strategies aimed at restoring mitochondrial function whether by enhancing biogenesis, reducing oxidative stress, or modulating cytokine-driven metabolic shifts could address the convergent mechanisms of coagulopathy, immune exhaustion, and respiratory decompensation. Finally, the identification of TFPI dysregulation suggests that targeted anticoagulant strategies should account for mitochondrial-driven endothelial dysfunction. A deeper understanding of the mechanisms driving respiratory decompensation may ultimately improve risk stratification, guide clinical decision-making, and enhance patient outcomes.

References

• 1.

  YukiKFujiogiMKoutsogiannakiS. COVID-19 pathophysiology: a review. Clin Immunol. (2020) 215:108427. doi: 10.1016/j.clim.2020.108427,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 2.

  MatthayMAZemansRLZimmermanGAArabiYMBeitlerJRMercatAet al. Acute respiratory distress syndrome. Nat Rev Dis Primers. (2019) 5:18. doi: 10.1038/s41572-019-0069-0

  • CrossRef
  • Google Scholar

• 3.

  BrownGCBorutaiteV. Inhibition of mitochondrial respiratory complex I by nitric oxide, peroxynitrite and S-nitrosothiols. Biochim Biophys Acta. (2004) 1658:44–9. doi: 10.1016/j.bbabio.2004.03.016,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 4.

  ArcherSLDasguptaAChenKHWuDBaidKMamatisJEet al. SARS-CoV-2 mitochondriopathy in COVID-19 pneumonia exacerbates hypoxemia. Redox Biol. (2022) 58:102508. doi: 10.1016/j.redox.2022.102508,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 5.

  MiddletonEAHeX-YDenormeFCampbellRANgDSalvatoreSPet al. Neutrophil extracellular traps contribute to immunothrombosis in COVID-19 acute respiratory distress syndrome. Blood. (2020) 136:1169–79. doi: 10.1182/blood.2020007008,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 6.

  MurthySGomersallCDFowlerRA. Care for Critically ill Patients with COVID-19. JAMA. (2020) 323:1499–500. doi: 10.1001/jama.2020.3633,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 7.

  AbateSMAhmed AliSMantfardoBBasuB. Rate of intensive care unit admission and outcomes among patients with coronavirus: a systematic review and Meta-analysis. PLoS One. (2020) 15:e0235653. doi: 10.1371/journal.pone.0235653,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 8.

  ChenH-Mvan der TouwWWangYSKangKMaiSZhangJet al. Blocking immunoinhibitory receptor LILRB2 reprograms tumor-associated myeloid cells and promotes antitumor immunity. J Clin Invest. (2018) 128:5647–62. doi: 10.1172/JCI97570,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 9.

  GeSXJungDYaoR. ShinyGO: a graphical gene-set enrichment tool for animals and plants. Bioinformatics. (2020) 36:2628–9. doi: 10.1093/bioinformatics/btz931,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 10.

  PomaznoyMHaBPetersB. GOnet: a tool for interactive gene ontology analysis. BMC Bioinformatics. (2018) 19:470. doi: 10.1186/s12859-018-2533-3,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11.

  PiñeroJBravoÀQueralt-RosinachNGutiérrez-SacristánADeu-PonsJCentenoEet al. DisGeNET: a comprehensive platform integrating information on human disease-associated genes and variants. Nucleic Acids Res. (2017) 45:D833–d839. doi: 10.1093/nar/gkw943,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 12.

  BardouPMarietteJEscudiéFDjemielCKloppC. jvenn: an interactive Venn diagram viewer. BMC Bioinformatics. (2014) 15:293. doi: 10.1186/1471-2105-15-293,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 13.

  ZuoYYalavarthiSShiHGockmanKZuoMMadisonJAet al. Neutrophil extracellular traps in COVID-19. JCI Insight. (2020) 5:e138999. doi: 10.1172/jci.insight.138999,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 14.

  GuarnieriJWDybasJMFazeliniaHKimMSFrereJZhangYet al. Core mitochondrial genes are down-regulated during SARS-CoV-2 infection of rodent and human hosts. Sci Transl Med. (2023) 15:eabq1533. doi: 10.1126/scitranslmed.abq1533,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15.

  Schulte-SchreppingJReuschNPaclikDBaßlerKSchlickeiserSZhangBet al. Severe COVID-19 is marked by a dysregulated myeloid cell compartment. Cell. (2020) 182:1419–1440.e23. doi: 10.1016/j.cell.2020.08.001,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 16.

  KellyCEThymiakouEDixonJETanakaSGodwinJEpiskopouV. Rnf165/Ark2C enhances BMP-Smad signaling to mediate motor axon extension. PLoS Biol. (2013) 11:e1001538. doi: 10.1371/journal.pbio.1001538,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 17.

  RixenSHavemeyerATyl-BielickaAPysniakKGajewskaMKuleckaMet al. Mitochondrial amidoxime-reducing component 2 (MARC2) has a significant role in N-reductive activity and energy metabolism. J Biol Chem. (2019) 294:17593–602. doi: 10.1074/jbc.RA119.007606,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 18.

  JakobsHHMikulaMHavemeyerAStrzalkowskaABorowa-ChmielakMDzwonekAet al. The N-reductive system composed of mitochondrial amidoxime reducing component (mARC), cytochrome b5 (CYB5B) and cytochrome b5 reductase (CYB5R) is regulated by fasting and high fat diet in mice. PLoS One. (2014) 9:e105371. doi: 10.1371/journal.pone.0105371,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 19.

  HavemeyerALangJClementB. The fourth mammalian molybdenum enzyme mARC: current state of research. Drug Metab Rev. (2011) 43:524–39. doi: 10.3109/03602532.2011.608682,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 20.

  ZhangJChenLHeYDingYZhouHHuHet al. Large-conductance calcium-activated potassium channels in the neurons of pre-Bötzinger complex and their participation in the regulation of central respiratory activity in neonatal rats. Neurosci Lett. (2010) 481:159–63. doi: 10.1016/j.neulet.2010.06.077,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 21.

  SzteynKSinghH. BK(ca) channels as targets for cardioprotection. Antioxidants (Basel). (2020) 9:760. doi: 10.3390/antiox9080760,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 22.

  LwaleedBABassPS. Tissue factor pathway inhibitor: structure, biology and involvement in disease. J Pathol. (2006) 208:327–39. doi: 10.1002/path.1871,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 23.

  Al-TamimiAOYusufAMJayakumarMNAnsariAWElhassanMAbdulKarimFet al. SARS-CoV-2 infection induces soluble platelet activation markers and PAI-1 in the early moderate stage of COVID-19. Int J Lab Hematol. (2022) 44:712–21. doi: 10.1111/ijlh.13829,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 24.

  HassaniSSayyadiMAlmasi-HashianiA. Plasma tissue factor pathway inhibitor levels in coronavirus disease 2019 patients: a systematic review and meta-analysis. Blood Coagul Fibrinolysis. (2024) 35:196–205. doi: 10.1097/MBC.0000000000001301,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 25.

  LabrecqueMBrunet-RatnasinghamEHamiltonLKAuldDMontpetitARichardsBet al. Transcriptomic profiling of severe and critical COVID-19 patients reveals alterations in expression, splicing and polyadenylation. Sci Rep. (2025) 15:13469. doi: 10.1038/s41598-025-95905-y,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 26.

  MosharafMPRezaMSKibriaMKAhmedFFKabirMHHasanSet al. Computational identification of host genomic biomarkers highlighting their functions, pathways and regulators that influence SARS-CoV-2 infections and drug repurposing. Sci Rep. (2022) 12:4279. doi: 10.1038/s41598-022-08073-8,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 27.

  LuLZhangHDaupharsDJHeYW. A potential role of interleukin 10 in COVID-19 pathogenesis. Trends Immunol. (2021) 42:3–5. doi: 10.1016/j.it.2020.10.012,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 28.

  ZhaoYQinLZhangPLiKLiangLSunJet al. Longitudinal COVID-19 profiling associates IL-1RA and IL-10 with disease severity and RANTES with mild disease. JCI Insight. (2020) 5:e139834. doi: 10.1172/jci.insight.139834,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 29.

  MartinodKWagnerDD. Thrombosis: tangled up in NETs. Blood. (2014) 123:2768–76. doi: 10.1182/blood-2013-10-463646,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 30.

  MeleroIVillalba-EsparzaMRecalde-ZamaconaBJiménez-SánchezDTeijeiraÁArguetaAet al. Neutrophil extracellular traps, local IL-8 expression, and cytotoxic T-lymphocyte response in the lungs of patients with fatal COVID-19. Chest. (2022) 162:1006–16. doi: 10.1016/j.chest.2022.06.007,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 31.

  ZhengMGaoYWangGSongGLiuSSunDet al. Functional exhaustion of antiviral lymphocytes in COVID-19 patients. Cell Mol Immunol. (2020) 17:533–5. doi: 10.1038/s41423-020-0402-2,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 32.

  VarchettaSMeleDOlivieroBMantovaniSLudovisiSCerinoAet al. Unique immunological profile in patients with COVID-19. Cell Mol Immunol. (2021) 18:604–12. doi: 10.1038/s41423-020-00557-9,

  • Pubmed Abstract
  • CrossRef
  • Google Scholar