Identification of robust mortality-associated neutrophil gene programs and cytokine responses using large-scale ARDS endotracheal aspirate scRNA-seq data

Identification of robust mortality-associated neutrophil gene programs and cytokine responses using large-scale ARDS endotracheal aspirate scRNA-seq data | medRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-P4HH5NV'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search Identification of robust mortality-associated neutrophil gene programs and cytokine responses using large-scale ARDS endotracheal aspirate scRNA-seq data Owen Whitley , Tim Delemarre , Miguel Reyes , Kimberly Kajihara , Hannah Little-Hooy , Cynthia Chen , Marco De Simone , Min Xu , Haridha Shivram , Jason Hackney , Diana Chang , Xia Gao , Eugene Kim , Nandhini Ramamoorthi , Rafal K. Woycicki , Shoshana Zha , Charles R. Langelier , Sharookh B. Kapadia , Carrie M. Rosenberger , Huwenbo Shi doi: https://doi.org/10.1101/2025.08.21.25333911 Owen Whitley 1 Department of Data & Analytics Chapter – Computational Science , Roche Canada , Mississauga, ON L5L 4M1, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site Tim Delemarre 2 Department of Infectious Diseases, Genentech Inc. , South San Francisco, CA 94080, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Miguel Reyes 2 Department of Infectious Diseases, Genentech Inc. , South San Francisco, CA 94080, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Kimberly Kajihara 2 Department of Infectious Diseases, Genentech Inc. , South San Francisco, CA 94080, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Hannah Little-Hooy 3 Department of Translational Immunology, Genentech Inc. , South San Francisco, CA 94080, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Cynthia Chen 4 Department of Microchemistry, Proteomics, Lipidomics and Next-Generation Sequencing, Genentech Inc. , South San Francisco, CA 94080, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Marco De Simone 4 Department of Microchemistry, Proteomics, Lipidomics and Next-Generation Sequencing, Genentech Inc. , South San Francisco, CA 94080, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Min Xu 3 Department of Translational Immunology, Genentech Inc. , South San Francisco, CA 94080, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Haridha Shivram 5 Department of Computational Sciences, Genentech Inc. , South San Francisco, CA 94080, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Jason Hackney 5 Department of Computational Sciences, Genentech Inc. , South San Francisco, CA 94080, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Diana Chang 6 Department of Human Genetics, Genentech Inc. , South San Francisco, CA 94080, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Xia Gao 7 Department of Immunology & Regenerative Medicine , Genentech Inc., South San Francisco, CA 94080, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Eugene Kim 7 Department of Immunology & Regenerative Medicine , Genentech Inc., South San Francisco, CA 94080, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Nandhini Ramamoorthi 8 Pharma Biosample Services, Genentech Inc. , South San Francisco, CA 94080, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Rafal K. Woycicki 9 Department of Data & Analytics Chapter – Computational Science , Roche Poland , 02-679 Warszawa, Poland Find this author on Google Scholar Find this author on PubMed Search for this author on this site Shoshana Zha 10 School of Medicine, University of California , San Francisco, San Francisco, CA 94143, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Charles R. Langelier 10 School of Medicine, University of California , San Francisco, San Francisco, CA 94143, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Sharookh B. Kapadia 2 Department of Infectious Diseases, Genentech Inc. , South San Francisco, CA 94080, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Carrie M. Rosenberger 7 Department of Immunology & Regenerative Medicine , Genentech Inc., South San Francisco, CA 94080, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Huwenbo Shi 5 Department of Computational Sciences, Genentech Inc. , South San Francisco, CA 94080, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: shi.huwenbo.genentech{at}gmail.com Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF Abstract Myeloid dysregulation is involved in severe Acute Respiratory Distress Syndrome (ARDS), but the mechanistic understanding is still incomplete. We analyzed endotracheal aspirate (ETA) scRNA-seq data from 40 intubated patients (both COVID-19 and non-COVID-19) in the COMET cohort, to uncover cell types and gene programs associated with ARDS mortality. In cell-type-agnostic analyses, we identified neutrophils as a critical cell type prognostic for ARDS mortality. Subsequent analyses uncovered a mortality-associated neutrophil subpopulation characterized by upregulation of IFN-ɣ response and ferroptosis genes, and downregulation of TNF-ɑ response genes. These findings were validated across independent COVID-19 clinical cohorts. Analyses of in vitro data further suggested that gene programs active in the mortality-associated neutrophil subpopulation could be attributable to specific cytokine responses. Finally, we confirmed that this mortality-associated neutrophil subpopulation could be induced in rodents via influenza A or SARS-CoV-2 infections, providing insights for modeling this neutrophil phenotype for therapeutic development. Introduction Acute Respiratory Distress Syndrome (ARDS) is a severe and critical respiratory disorder, characterized by acute lung inflammation and injury, hypoxemia, edema, and significant decline in lung function. 1 – 3 The causes of ARDS are highly heterogeneous, including sepsis, bacterial pneumonia, and viral respiratory infection. 1 , 3 – 6 With a mortality rate of approximately 40%, 1 , 3 , 5 , 7 and a high financial burden on patients, 8 , 9 treatment of ARDS represents an unmet medical need in the United States. Standard of care treatments for ARDS typically include prone positioning, mechanical ventilation, and extracorporeal membrane oxygenation (ECMO). 1 , 6 , 10 , 11 Additional therapeutics, including the use of corticosteroids, have also been investigated for ARDS. 12 – 16 However, existing therapeutics for ARDS generally have limited efficacy on reducing ARDS mortality, due to patient heterogeneity and complex ARDS disease biology. 15 , 17 , 18 Thus, we need to develop new therapeutics that target specific ARDS pathophysiological mechanisms in well-defined patient populations. A key biological process involved in ARDS is dysregulated immune responses. 6 , 16 , 19 – 24 Recent work analyzing both COVID-19 and non-COVID-19 ARDS have implicated myeloid dysregulation as a key driver of ARDS onset and progression. 20 – 22 , 25 , 26 Most of these studies focused on monocytes and macrophages. Despite the critical role of neutrophils in ARDS, 6 , 27 – 30 the understanding of the impact of specific neutrophil subpopulations, functions, and phenotypes on patient outcome, particularly mortality, is currently limited. This is in part because assaying gene expression in neutrophils is challenging. 31 – 33 However, recent technological advancements and efforts during the COVID-19 pandemic have led to an increase in neutrophil single-cell RNA-seq (scRNA-seq) studies focused on ARDS. 23 , 34 , 35 Here, we analyzed endotracheal aspirate (ETA) scRNA-seq data (129K cells) from 40 intubated patients (both COVID-19 and non-COVID-19) in the COVID-19 Multiphenotyping for Effective Therapies (COMET) cohort. To understand the cell populations and gene programs involved in ARDS mortality in the lower airways, we performed both cell-type-agnostic and neutrophil-focused analyses using state-of-the-art computational methods, including CNA 36 and cNMF. 37 In the cell-type-agnostic analyses, we determined that airway neutrophils were a critical cell type prognostic for ARDS mortality. In the neutrophil-focused analyses, we identified a mortality-associated neutrophil subpopulation, in which up-regulated genes were involved in IFN-ɣ response (e.g., TXNIP , IFI30 , TRIM25 , STAT1 , RNF213 , and TMEM140 ), and ferroptosis (e.g., PIM1 , HMOX1 , TXNRD1 , GCLM , and FTL ), while down-regulated genes were involved in TNF-ɑ signaling (e.g., PLAUR , IL1B , CCL20 , GCH1 , ICAM1 , NFKBIA , TNFAIP6 , and TNFAIP3 ). We confirmed that this mortality-associated neutrophil subpopulation was associated with both COVID-19 and non-COVID-19 ARDS mortality, and was generally replicated in independent COVID-19 clinical data. Next, we determined that gene programs active in the mortality-associated neutrophil subpopulation could be recapitulated by cytokine stimulation in vitro . Finally, we confirmed that the mortality-associated neutrophil subpopulations could be induced in rodents in flu or SARS-CoV-2 infection models. Overall, our work provides additional insights into the role of neutrophils in ARDS. Results Cell-type agnostic analysis of 30-day mortality using the COMET ETA scRNA-seq data We investigated the cell types involved in 30-day mortality (whether a patient died within 30 days of hospital admission), a pertinent ARDS clinical endpoint. To that end, we analyzed the endotracheal aspirate (ETA) scRNA-seq data collected from 40 (N fatality =10, N survivor =30) intubated COVID-19 and non-COVID-19 ARDS patients (see Table 1 , Supplementary Table 1, Supplementary Figure 1). We applied Co-varying Neighborhood Analysis 36 (CNA) to analyze the baseline samples (i.e., earliest samples taken since intubation), which consisted of 128,822 cells across 5 broad cell types. On average, these ETA samples were collected 2.1 days (s.d. 2.1) post ICU admission; CNA is a computational method to detect disease-associated cell states/populations at cell neighborhood resolution, and is more powerful than cell-type/cell-cluster based approaches. 36 We elected to analyze the baseline samples to enrich for potential targetable biology and biomarkers for ARDS; later samples were more susceptible to the consequential effects of the disease and clinical interventions. We adjusted for known ARDS risk factors, e.g., age, sex, COVID-19 status, as well as other technical covariates, e.g., batch, time of sample collection in our analyses (see Methods). We elected not to adjust for corticosteroid use in our primary analyses, as this information was missing for 4 patients; and all but 1 of the remaining 36 patients received corticosteroid ( Table 1 a ). View this table: View inline View popup Table 1: Overview of the COMET ETA scRNA-seq data. ( a ) We report the number and proportion of the patients across demographic stratifications (e.g., age, sex, etc.), clinical outcomes (e.g., 30-day mortality status, max ordinal scale score), and treatment status (e.g., steroids and Tocilizumab). ( b , c ) We report the number and proportion of the major COMET ETA scRNA-seq cell types for all and non-COVID-19 patients, respectively. *Primary diagnoses of non-COVID-19 patients can be found in Supplementary Table 1. View this table: View inline View popup Download powerpoint Table 2: Results from gene set enrichment analysis of the top 100 genes with the highest spectra scores of each cNMF topic. We report the top 5 MSigDB Hallmark gene sets (FDR < 0.05) for each gene program, except for gene program 4, which yielded less than 5 significant results. We report the complete results in Supplementary Table 7. The results are reported in Figure 1 , Supplementary Figure 2–3. Overall, airway neutrophils and lymphocytes had significantly positive (p < 0.01) average CNA neighborhood correlations ( Figure 1 b, d ), suggesting an enrichment of these cell types in 30-day mortality cases. In contrast, other cell types (e.g., monocytes and epithelial cells) had significantly negative (p < 0.01) average CNA neighborhood correlation, indicating a depletion of these cell types in 30-day mortality cases. We note that although Red Blood Cells (RBCs) also had significantly positive (p < 0.01) average CNA neighborhood correlation, the proportion of RBCs was low in the scRNA-seq data ( Table 1 b , c ; Supplementary Figure 2). We did not observe statistically significant differences in immune cell type proportions across 30-day mortality cases vs. controls (Supplementary Figure 2a). Although we observed a significant differential RBC proportion, the representation of RBC in the scRNA-seq data was limited ( Table 1 b, c , Supplementary Figure 2a). Download figure Open in new tab Figure 1: Results from cell-type agnostic analysis of 30-day mortality using the COMET ETA scRNA-seq data. ( a ) Flowchart showing the workflow of the overall study. ( b , c ) UMAP plots showing the CNA neighborhood correlations for 30-day mortality obtained using all and non-COVID-19 patients, respectively. ( d , e ) Violin plots showing the distribution of CNA neighborhood correlations across cell types in ( b ) and ( d ). respectively. The white lines in the violin plots represent the median, with boxes representing the interquartile range. The areas covered by the violins are proportional to the number of cells in each cell type. We report permutation p-values (based on 1,000 permutations), testing whether the average CNA neighborhood correlations of each cell type is significantly different from zero, on the top of each violin plot. Next, we applied CNA specifically to 30,023 cells from the 12 non-COVID-19 patients (N fatality =3, N survivor =9), adjusting for age and sex, as well as technical covariates as above. Here, neutrophils had the highest average CNA neighborhood correlation ( Figure 1 c, e ). Interestingly, the average lymphocyte CNA neighborhood correlation was much lower compared with the all-patient analysis ( Figure 1 d , e ), while the results for other cell types were generally similar to the all-patient analysis ( Figure 1 d, e ). We also did not observe statistically significant differences in cell type proportions across 30-day mortality cases vs. controls in the non-COVID-19 analysis (Supplementary Figure 2b). We also performed a secondary analysis, in which corticosteroids use was included as an additional covariate to the ones adjusted in our primary analyses, using the 36 patients for whom this information is available. We observed similar disease-associated cell types and neutrophil subpopulations as those obtained without adjusting for corticosteroids use ( Figure 1 b , d , Supplementary Figure 3a, b ). This is likely because the results were dominated by the patients (35 out of 36), who were treated with corticosteroids. In summary, our CNA analysis of 30-day mortality using the COMET ETA scRNA-seq data suggested that neutrophils were strongly associated with severe ARDS outcomes in both COVID-19 and non-COVID-19 patients. ETA neutrophil-focused CNA analysis of 30-day mortality Next, we performed neutrophil-focused analyses to identify neutrophil-specific ARDS disease biology, motivated by the strong association of ETA neutrophils with 30-day mortality in the cell-type agnostic analysis regardless of disease etiology ( Figure 1 d , e ). We applied CNA to the 45,572 baseline ETA neutrophils across 37 (N fatality =8, N survivor =29) COVID-19 and non-COVID-19 patients post QC in the COMET scRNA-seq data ( Table 1 ), adjusting for similar covariates as in the cell-type agnostic analyses (see above, Methods). We also identified neutrophil-specific mortality-associated genes that were correlated with the CNA neighborhood correlations across the neutrophils (see Methods). We elected not to perform lymphocyte-focused analysis, as the association with 30-day mortality was much weaker in the non-COVID-19 vs. all-patient analysis ( Figure 1 d , e ); the relatively small number of lymphocytes ( Table 1 ) also limited the statistical power of lymphocyte-focused analysis. The results are reported in Figure 2 , Supplementary Figure 4–11, and Supplementary Tables 2–6. The association with 30-day mortality among the neutrophils exhibited a gradient pattern, with a subset of neutrophils exhibiting higher CNA neighborhood correlation than others ( Figure 2 ). As expected, the average CNA neighborhood correlation was higher for deceased vs. survived patients (Supplementary Figure 4). Download figure Open in new tab Figure 2: Results from the neutrophil-focused analysis of the COMET ETA scRNA-seq data using all patients. ( a ) We report the neutrophil-specific CNA neighborhood correlations for 30-day mortality obtained using all patients. ( b ) We report the correlation between the expression of each gene and CNA neighborhood correlations across the neutrophils. Genes with significant correlation (|R| > 0.1, FDR < 0.05) are highlighted as red and blue dots. ( c , d ) We report the top 5 MSigDB Hallmark (left) and WikiPathways (right) gene sets significantly (FDR < 0.05) overlapping genes positively and negatively correlated with the CNA neighborhood correlations, respectively. ( e , f ) We report the ssGSEA scores calculated using CNA mortality-increasing and decreasing genes, respectively, based on pseudobulked neutrophil expression from all patients (left panels) and non-COVID-19 patients (right panels). P-values were obtained using Wilcoxon rank-sum tests. ( g ) We report the top 20 genes with the highest spectra scores for gene program 1 (mortality-increasing genes) and program 2 (mortality-decreasing genes). The usage of the gene programs are shown in the UMAP plots. ( h ) We report the distribution of CNA neighborhood correlations in each quintile of the usages of gene program 1 and 2. The dark lines in the box plots represent the median, with the boxes representing the interquartile range (IQR), and whiskers representing 1.5× IQR. However, many survived patients also harbored neutrophils with high CNA neighborhood correlation, and vice versa. This is expected, as CNA is designed to identify subtle neutrophil states associated with mortality, which could exist in both survived and decreased patients. We identified 72 and 60 (out of 32,877) genes significantly positively (FDR 0.1) and negatively (FDR < 0.05, R < −0.1) correlated with the CNA neighborhood correlations across the cells (Supplementary Table 2), using similar approaches as introduced in ref. 36 . We refer to these genes as “CNA mortality-increasing genes” and “CNA mortality-decreasing genes” hereafter, respectively. Similar analyses using non-COVID-19 ARDS patients yielded smaller sets of mortality-associated genes that significantly overlapped those from the all-patient analysis (Supplementary Table 3, Supplementary Figure 5, 6 ). Thus, we elected to use the larger gene sets derived from the better powered all-patient analysis for downstream analyses. We performed several downstream analyses using the CNA mortality-associated genes. First, we performed over-representation analyses of the mortality-associated genes using GSEApy 38 , analyzing MSigDB Hallmark and WikiPathways gene sets. Both CNA mortality-increasing and mortality-decreasing genes significantly overlapped cytokine signaling gene sets, e.g., IFN-ɑ and IFN-ɣ response for the CNA mortality-increasing signature, and TNF-ɑ signaling via NF-κB for the CNA mortality-decreasing signature ( Figure 2 , Supplementary Table 4). We note the well-documented literature suggesting myeloid cells could be the cellular sources of IFN-ɣ 39 , 40 and TNF-ɑ. 41 , 42 Interestingly, the CNA mortality-increasing genes also significantly overlapped with reactive oxygen species and ferroptosis-related gene sets ( Figure 2 c , Supplementary Table 4 a , c ). Second, we quantified the overall activity of the mortality-associated genes in each patient by applying ssGSEA 43 , 44 to the pseudobulked neutrophil gene expression; higher ssGSEA scores indicate higher expression of the mortality-associated genes relative to other genes. As expected, ssGSEA scores for mortality-increasing genes were generally higher in 30-day mortality cases than controls, and vice versa for mortality-decreasing genes ( Figure 2 e , f ). This was the case in both non-COVID-19 patients and in all patients ( Figure 2 e , f ). We note that the lack of statistical significance in analysis of non-COVID-19 patients was likely attributable to a limited number of mortality cases. Third, we investigated the cell types in which the CNA mortality-associated genes were active. We calculated AUCell 45 scores for all the cells in the ETA scRNA-seq data using the CNA mortality-associated genes; higher AUCell scores represent higher gene activity. As expected, the activity of the CNA mortality-associated genes were specific to myeloid cells (Supplementary Figure 7). However, the activity of the CNA mortality-decreasing genes were generally specific to neutrophils, whereas the CNA mortality-increasing genes were also broadly active in monocytes, macrophages, and dendritic cells (Supplementary Figure 7). We also note that neutrophils with high AUCell scores for CNA mortality-increasing genes and high CNA neighborhood correlations did not fully overlap ( Figure 2 a , Supplementary Figure 7e, 8 ). This was likely because genes characterizing neutrophils with high CNA neighborhood correlations were not captured by the gene signature, e.g., due to limited statistical power. We also performed several additional secondary analyses to assess the robustness of the mortality-associated gene signatures. First, we calculated ssGSEA scores for each patient using mortality-associated genes derived from the remaining patients in a leave-one-out framework (Methods). The resulting ssGSEA scores for CNA mortality-increasing genes were generally higher in 30-day mortality cases than controls, and vice versa for mortality-decreasing genes (Supplementary Figure 9), consistent with the main results ( Figure 2 f ). Second, we performed differential expression (DE) analysis contrasting 30-day mortality cases vs. controls via DESeq2, 46 using the pseudobulked neutrophil gene expression (see Methods). As expected, the effects of the genes on mortality from pseudobulk DE analysis were generally consistent with those from CNA analysis (Supplementary Figure 10). Third, to identify proteins potentially driving the expression of the mortality-associated genes, we compared the ssGSEA scores for the CNA mortality-associated genes with the protein expression of IFN-ɣ and TNF-ɑ in ETA assayed using Olink 47 (“Methods”). We note that IFN-ɑ, which was also implicated in our over-representation analysis ( Figure 2 c ), is not included in the Olink protein assay panel. We observed weakly positive correlations between the ssGSEA scores and protein expression, with statistically significant result observed for the correlation between mortality-decreasing genes and TNF-ɑ (Supplementary Figure 11, Supplementary Table 6). In summary, we identified a robust ETA neutrophil subpopulation associated with 30-day mortality in both COVID-19 and non-COVID-19 patients, characterized by up-regulation of genes involved in IFN-ɑ, IFN-ɣ response pathways, and ferroptosis, and down-regulation of genes involved in TNF-ɑ response pathway. Gene program analysis of ETA neutrophils We investigated whether specific neutrophil functions and cellular processes were associated with 30-day mortality using consensus Nonnegative Matrix Factorization (cNMF). 37 Briefly, cNMF decomposes the gene expression matrix from a single-cell study into gene programs, a set of genes that collectively contribute to a specific biological function or cellular process, and inferring their overall activity (i.e., usage of gene programs) in each individual cell. Unlike CNA, cNMF identifies gene programs in an unsupervised fashion, without explicitly modeling the disease outcome. We applied cNMF to the 45,572 neutrophils from all 37 patients (baseline samples) in the COMET ETA scRNA-seq data using cNMF to extract consistent gene programs across COVID-19 and non-COVID-19 patients (“Methods”). The results of the cNMF analyses are reported in Figure 2 , Supplementary Figures 12–14, and Supplementary Tables 7–8. We set the cNMF k parameter, the pre-specified number of gene programs, to 4, based on the cNMF k selection plot, balancing between stability of the factorization and expression matrix reconstruction error (Supplementary Figure 12a). We determined that the median usages of gene programs 1 and 2 were significantly (p=0.001) positively and negatively, respectively, correlated with the median AUCell scores for the CNA mortality-increasing genes, across all and non-COVID-19 patients ( Figure 2 g , Supplementary Figure 14). Notably, gene programs 1 and 2 were also the gene programs with the highest usages (Supplementary Figure 12b, d ). Gene set enrichment analyses of the top 100 genes with the highest spectra scores for these gene programs revealed a set of shared pathways (e.g., IFN-ɣ response and TNF-ɑ signaling via NF-κB for both programs; Supplementary Table 8 a , b , e , f ), as well as a set of gene program specific pathways (e.g., ferroptosis and NRF2 pathways for gene program 1 only; Supplementary Table 8 e ). Despite the significant correlations, we note that the neutrophils with high CNA neighborhood correlations did not fully overlap those with high usage of gene program 1 and 2, suggesting that mortality-associated neutrophils likely involved multiple different cellular processes. We elected not to analyze gene programs 3 and 4 due to their relatively low usage across the ETA neutrophils (Supplementary Figure 12b, d ). In summary, our cNMF analysis of the ETA neutrophils suggested that ETA neutrophils expression could be captured by gene programs involving IFN-ɣ and TNF-ɑ signaling pathways, and that these gene programs partially explained the cell populations associated with 30-day mortality identified by CNA. Analysis of mortality-associated genes in independent COVID-19 clinical scRNA-seq data We assessed the replicability of the CNA mortality-associated genes derived from the COMET ETA scRNA-seq data, in the publicly available COVID-19 scRNA-seq data from Bost et al. 48 The Bost et al. data included 26,336 bronchoalveolar lavage fluid (BALF) and 21,652 blood neutrophils from 21 severe COVID-19 patients (N fatility =8, N survivor =12 for BALF; N fatility =7, N survivor =13 for blood; Supplementary Table 9). For each patient in the Bost et al. BALF or blood data, we calculated aggregated ssGSEA 43 , 44 scores for CNA mortality-increasing and decreasing genes, based on the pseudobulked neutrophil gene expression. The results for the analyses of the Bost et al. data are reported in Figure 3 , Supplementary Figures 15–17, Supplementary Tables 10–11. In the analysis of the BALF neutrophils, the ssGSEA scores for CNA mortality-increasing were significantly positively associated with the SOFA scores (Spearman’s ⍴=0.49 (p=0.029), Figure 3 a ), and remained significant after adjusting for age and sex (Supplementary Table 10 a ). The association with mortality was also positive, but not statistically significant (p=0.067, Supplementary Figure 17a). We did not detect significant association between the ssGSEA scores for CNA mortality-decreasing genes and SOFA scores or mortality ( Figure 3 b , Supplementary Figure 16b, 17 b , f ). In the analysis of the blood neutrophils, we did not detect a significant association between the ssGSEA scores for CNA mortality-associated genes and SOFA scores or mortality ( Figure 3 c , d , Supplementary Figure 16c, d , 17 c , d , g , h ). Download figure Open in new tab Figure 3: Results from the analysis of CNA mortality-associated genes in the neutrophils from the Bost et al. COVID-19 patient scRNA-seq data. ( a , b ) We report the Spearman’s correlation between patient SOFA clinical severity scores and ssGSEA scores for CNA mortality-increasing and decreasing genes, respectively, for BALF; we report analogous results for blood neutrophils in ( c , d ). ssGSEA scores were calculated based on log-normalized pseudobulk gene expression. Solid blue lines in the scatter plots are regression lines fitting the SOFA scores against the ssGSEA scores; shaded regions represent 95% confidence intervals on both sides. As a secondary analysis, we also analyzed the CNA mortality-associated genes in a smaller publicly available COVID-19 scRNA-seq data from Liao et al., 34 which included 3,628 BALF neutrophils from 9 patients (Supplementary Table 12, Supplementary Figure 18). We calculated aggregated scores for CNA mortality-increasing and decreasing genes for each cell using AUCell. 45 The aggregated CNA mortality-increasing and decreasing gene signature scores were both significantly higher in severe vs. mild patients (Supplementary Figure 19b, d ). In summary, the CNA mortality-increasing gene signatures derived from the analysis of 30-day mortality using the COMET ETA neutrophil scRNA-seq data were generally replicated in BALF, but not blood, neutrophils from independent COVID-19 scRNA-seq data. However, replication of the CNA mortality-decreasing gene signatures were generally weaker. We discuss ramifications of these results in Discussion. Analysis of mortality-associated genes using in vitro models We next assessed the effects of TNF-ɑ, IFN-ɣ, dexamethasone, and lipopolysaccharide (LPS), on the expression of the CNA mortality-associated genes in human blood neutrophils in in vitro experiments; TNF-ɑ and IFN-ɣ are cytokines implicated in our pathway analyses, whose protein expression correlated with the expression of CNA mortality-associated genes in ETA (see subsection above); dexamethasone is a standard of care treatment for hospitalized COVID-19 patients and showed efficacy in reducing mortality; 13 LPS is a pro-inflammatory stimulus for inducing sepsis in animal models. 49 In the neutrophils treated with IFN-ɣ, 19 of the 72 CNA mortality-increasing genes were significantly (|log 2 FC| > 0.5, FDR < 0.05) up-regulated and 13 down-regulated relative to medium ( Figure 4 a , c ). The up-regulated CNA mortality-increasing genes exhibited significantly (p=0.03) stronger magnitude of log 2 FC than down-regulated genes (Supplementary Figure 20a). In contrast, the magnitude of the log 2 FC was similar for the 14 up-regulated and 18 down-regulated (out of 60) CNA mortality-decreasing genes (Supplementary Figure 20e) Download figure Open in new tab Figure 4: Results from the analysis of CNA mortality-associated genes in the in vitro blood neutrophil experiments. ( a, b ) We report the distributions of the log 2 fold change (FC) of the expression of CNA mortality-increasing and decreasing genes relative to media, respectively, in neutrophils stimulated with dexamethasone (dex), TNF-ɑ, lipopolysaccharide (LPS), and IFN-ɣ. ( c , d ) We report the values of the log 2 FC across the 4 contrasts, for the CNA mortality-increasing and decreasing genes, respectively. Red and blue dots in ( a ) and ( b ) represent significantly (|log 2 FC| > 0.5, FDR < 0.05) up and down regulated genes, respectively; black dots and error bars in the center of the violin plots in ( a , b ) represent the average log 2 FC and 95% confidence intervals (calculated using 1,000 bootstraped samples), respectively. Significant differentially expressed genes are marked by stars (“★”) in ( c ) and ( d ). Numerical results are reported in Supplementary Table 13. In the neutrophils treated with TNF-ɑ, LPS, and dexamethasone, 9, 10, and 10 CNA mortality-increasing genes were up-regulated, respectively ( Figure 4 a , c ). Importantly, these were generally distinct from those up-regulated by IFN-ɣ; 14, 6, and 15 genes were down-regulated, respectively ( Figure 4 a , c ). We did not observe significant differences in the magnitude of log 2 FC between up-regulated and down-regulated CNA mortality-increasing genes in the neutrophils treated with these stimuli(Supplementary Figure 20b, c , d ). Among the 60 CNA mortality-decreasing genes, 53 (0), 50 (0), and 6 (38) genes were up-regulated (down-regulated), respectively, in neutrophils treated with TNF-ɑ, LPS, and dexamethasone ( Figure 4 b , d , Supplementary Figure 20f, g , h ). We concluded that CNA mortality-increasing genes were likely induced by a combination of cytokines, whereas CNA mortality-decreasing genes could be induced by TNF-ɑ alone. Effects of IFN-ɣ and TNF-ɑ on neutrophil function and cytokine secretion Motivated by the results from the analyses of transcriptomics data, we next performed in vitro experiments to assess the effects of IFN-ɣ, TNF-ɑ, and LPS on human neutrophils functions. We stimulated human blood neutrophils from 6 donors with TNF-ɑ, LPS, and IFN-ɣ, as well as combinations of IFN-ɣ with TNF-ɑ or LPS to mimic the complex cytokine environments in ARDS patients. We assessed cytokine release, and various immunological functions, including activation (indicated by reduced CD62L and CD16 ), degranulation (indicated by increased CD63 and CD66b ), and extracellular elastase activity. The results are reported in Figure 5 . We observed that IFN-ɣ alone had limited effect on neutrophil activation or degranulation ( Figure 5 a – c ), while TNF-ɑ and LPS substantially induced both neutrophil functions ( Figure 5 a – c ). However, when combined with either TNF-ɑ or LPS, IFN-ɣ enhanced neutrophil activation and degranulation ( Figure 5 a , b ), indicating a synergistic effect. Observations of increased degranulation were also consistent with the elevated (although not statistically significant) activity of extracellular elastase, a byproduct of degranulation, 50 , 51 in these conditions ( Figure 5 d ). Download figure Open in new tab Figure 5: Functional and cytokine profiles of stimulated human primary blood neutrophils. ( a ) We report the proportions of neutrophils that were CD16 - (left) and CD62L - (right), markers for neutrophil activation. ( b ) We report the proportions of neutrophils that were CD63 + (left) and CD66b + (right), markers for neutrophil degranulation. ( c ) Flow cytometry plot showing neutrophil activation ( CD62L - , x-axis) vs. degranulation ( CD63 + , y-axis,) for IFN-ɣ and TNF-ɑ stimulated neutrophils. ( d ) We report the fold change in measured elastase in stimulated vs. unstimulated (medium) neutrophils. ( e ) We report the log 2 fold change in measured cytokine secretion by stimulated vs. unstimulated (medium) neutrophils. ( f , g , h ) We report the concentration (in picograms per milliliter) of CCL3 , CCL4 , and IL-8 , respectively, in neutrophils stimulated with TNF-ɑ vs. IFN-ɣ + TNF-ɑ (left panels) and LPS vs. IFN-ɣ + LPS (right panels). Gray squares in the strip plots represent the average values; error bars represent one standard deviation on both sides. All p-values were obtained via the Wilcoxon signed-rank tests. Statistically significant results (p < 0.05) are denoted by “★”. We also observed up-regulation of CCL3 , CCL4 , and IL-8 , cytokines among the CNA mortality-decreasing genes, in neutrophils stimulated with TNF-ɑ or LPS, and down-regulation of these cytokines in neutrophils stimulated with IFN-ɣ ( Figure 5 e ), consistent with our analysis of the bulk RNA-seq data (see subsection above). Notably, the expression of these cytokines were lower in neutrophils co-stimulated with IFN-ɣ and TNF-ɑ or LPS, compared to those stimulated with TNF-ɑ or LPS alone ( Figure 5 f – h ), suggesting that IFN-ɣ could downregulate the production of these cytokines despite its synergistic effects on activation or degranulation. Finally, although IFN-ɣ had limited effect on activation, degranulation, or elastase release ( Figure 5 a – c ), it up-regulated several cytokines, including CCL2 , CCL7 , CCL22 , IL-17a , and IL-1 ( Figure 5 e ), indicating an impact on inflammatory signaling and further neutrophil recruitment. 52 – 55 These results, along with our bulk RNA-seq results, suggest that IFN-ɣ and TNF-ɑ could synergistically influence neutrophil functions and gene programs, leading to a complex response that surpasses the effects of either cytokine alone. IFN-ɣ increased and TNF-ɑ decreased neutrophil ferroptosis We investigated the role of IFN-ɣ and TNF-ɑ on neutrophil ferroptosis, a significant biological process implicated by pathway analysis of the CNA mortality-increasing signature (Supplementary Table 4). We analyzed the expression of both ferroptosis promoting genes, specifically, LPCAT3 , ACSL4 , ALOX15 , ALOX5 , SAT1 , 56 – 61 and ferroptosis inhibiting genes, specifically, GPX4 and GCLC , 62 in the bulk RNA-seq data for IFN-ɣ and TNF-ɑ treated neutrophils. The results are reported in Figure 6 and Supplementary Table 13. IFN-ɣ significantly (FDR < 0.05) up-regulated ferroptosis promoting genes and down-regulated ferroptosis inhibiting genes, whereas TNF-ɑ had the opposite effect ( Figure 6 a ). We note that the up/down regulation of ferroptosis promoting/inhibiting genes only entailed a ferroptosis-prone state, and did not confirm that the cells were actively undergoing ferroptosis. Download figure Open in new tab Figure 6: Transcriptional and functional profiles of ferroptosis in stimulated human primary blood neutrophils. ( a ) We report the results from DE analysis of bulk RNA-seq data of neutrophils stimulated with IFN-ɣ vs. TNF-ɑ. Significant DE genes (|log 2 FC| > 0.5, FDR < 0.05) are highlighted in colors. Red and blue dots represent significant DE ferroptosis-inducing and inhibiting genes, respectively. ( b, c ) We report the membrane lipid peroxidization rates and relative viability (assayed via PI staining) of stimulated neutrophils. Each dot represents a donor. Gray squares in ( b , c ) represent the average values; error bars represent one standard deviation on both sides. Statistically significant differences are denoted by “★” (for p < 0.05) and “★★” (for p < 0.01). To further investigate the role of IFN-ɣ and TNF-ɑ in ferroptosis, we measured oxidized lipids, a marker of ferroptosis, on the membrane of stimulated human primary neutrophils, across 6 biological replicates from different donors via flow cytometry. We observed that membrane lipid peroxidation was significantly increased (p<0.05) in IFN-ɣ stimulated neutrophils, and decreased (p<0.05) in TNF-ɑ stimulated neutrophils ( Figure 6 b ). The membrane lipid peroxidation in IFN-ɣ stimulated neutrophils was reduced, after ferrostatin-1, a specific inhibitor of ferroptosis, 63 was added ( Figure 6 b ), confirming that the increase in membrane lipid peroxidation was indeed due to ferroptosis. We also observed that IFN-ɣ significantly decreased neutrophil viability, whereas TNF-ɑ did not ( Figure 6 c ). However, adding ferrostatin-1 in IFN-ɣ treated neutrophils did not completely prevent ferroptotic cell death ( Figure 6 c ), suggesting that ferrostatin-1 was not sufficient to prevent cytotoxic effects of IFN-ɣ. In summary, our results suggested that IFN-ɣ and TNF-ɑ could be potential inducers and inhibitors, respectively, of ferroptosis in human primary neutrophils. Analysis of mortality-associated genes using in vivo models We investigated whether the mortality-associated genes derived from the human data could be recapitulated using in vivo models; being able to model the mortality-associated gene signatures in vivo can facilitate the development of therapeutics for ARDS. We analyzed publicly available BALF scRNA-seq data from Peidli et al., 64 which includes 12 hamsters infected with SARS-CoV-2 of varying doses (high vs. low) across 2 time points (3 biological replicates per (time point, dosage) combination; Supplementary Table 14). We also analyzed BALF and blood scRNA-seq data from a mouse influenza A (PR8) model, which includes 18 mice (3 replicate per condition, pooled into 1 sequencing batch) infected with 0, 10 3 , and 10 4 plaque-forming unit (PFU) of flu virus, with samples collected at day 3 and day 7 (Supplementary Table 14, 15 , Supplementary Figure 21; Methods). For the Peidli et al. scRNA-seq data, we calculated ssGSEA 43 , 44 scores for the CNA mortality-associated genes based on pseudobulked neutrophil gene expression of each hamster sample. For the mouse flu model scRNA-seq data, we calculated AUCell scores using the CNA mortality-associated genes for each individual cell; we were not able to calculate aggregated scores for each mouse sample as we did with the Peidli et al. data, as cells from multiple mice were pooled together during sequencing. The results are reported in Figure 7 and Supplementary Figures 21–24. In the analysis of the Peidli et al. data, we observed significantly higher ssGSEA scores for CNA mortality-increasing genes in hamsters infected with high vs. low doses of SARS-CoV-2 at both day 2 and day 3 ( Figure 7 a ); we did not observe significant changes in the ssGSEA scores for CNA mortality-decreasing genes across SARS-COV-2 doses at either time points ( Figure 7 b ). Download figure Open in new tab Figure 7: Results from the analysis of CNA mortality-associated genes in BALF and blood neutrophils from in vivo model scRNA-seq data. ( a , b ) We report the ssGSEA scores calculated using the CNA mortality-increasing and decreasing genes (mapped to mouse orthologs), respectively, across the 4 (dosage, timepoint) pairs, based on the pseudobulked BALF neutrophil gene expression, derived from the Peidli et al. hamster SARS-CoV-2 model scRNA-seq data. P-values were obtained via two-sided t-tests. Each dot represents a hamster replicate. Gray squares represent the average ssGSEA scores; error bars represent one standard deviation on both sides. ( c , d ) We report the AUCell scores calculated using the CNA mortality-increasing and decreasing genes, respectively, across 4 the (dosage, timepoint) pairs, based on the expression of each BALF neutrophil in the mouse flu model scRNA-seq data; analogous results obtained based on blood neutrophils are shown in ( e , f ). P-values were obtained via 1,000 permutations. Each dot represents a cell. The boxes within the violin plots represent the interquartile ranges; white solid lines in the boxes represent the median. In the analysis of BALF neutrophils from our mouse influenza A model (Supplementary Figure 22), we observed that the AUCell scores for CNA mortality-increasing genes significantly increased with viral doses (p < 0.001) at both time points ( Figure 7 c , Supplementary Figure 23a, 24 a ), with stronger effect at day 7 (Hedges’ g = 0.89) than day 3 (Hedges’ g = 0.22). We also observed significantly lower AUCell scores for CNA mortality-decreasing genes (p 0.1) ( Figure 7 d , Supplementary Figure 23b, 24 b ). We elected to exclude non-infected samples in our analyses of BALF data, as they yielded less than 50 neutrophils. In contrast, in the analysis of CNA mortality-increasing genes in blood neutrophils, we observed much weaker effect sizes in samples infected with high vs. low flu dose (Hedges’ g = 0.18 and −0.034 at day 7 and day 3, respectively; Figure 7 e , Supplementary Figure 23c, 24 c ). These results suggest that the association between mortality and CNA mortality-increasing genes may manifest in a tissue-specific manner (also see “Discussion”). Interestingly, we also observed significant lower AUCell scores for CNA mortality-decreasing genes in mice infected with high flu dose ( Figure 7 f , Supplementary Figure 23d, 24 d ). In summary, the mortality-associated gene signatures derived from analysis of human data were generally recapitulated in in vivo models, which can be used in future studies to test therapeutic hypotheses. Discussion In this work, we highlight airway neutrophils as a key cell type associated with 30-day ARDS mortality, and identify mortality-associated neutrophil gene programs and signatures in the lower airway using CNA and cNMF. Notably, the mortality-associated neutrophil gene signatures were shared in ARDS patients of diverse etiology and not specific to COVID-19, and were also observed in independent COVID-19 patient cohorts and in rodents infected with SARS-CoV-2 or influenza virus. Pathway analyses implicated IFN-ɣ response and ferroptosis for CNA mortality-increasing genes, and TNF-ɑ response for mortality-decreasing genes. In in vitro experiments and functional assays of primary human neutrophils, IFN-ɣ upregulated a subset of the CNA mortality-increasing and ferroptosis-inducing genes, increased membrane lipid peroxidation (a marker of ferroptosis), while downregulating ferroptosis-inhibiting genes and decreasing neutrophil viability. Conversely, TNF-ɑ upregulated a subset of the CNA mortality-decreasing genes, while downregulating ferroptosis-inducing genes. Our in vitro experiments also suggest that IFN-ɣ and TNF-ɑ may interact to enhance the production of several chemokines, including CCL3 , CCL4 , and CXCL8, which could amplify recruitment of inflammatory cells to the lung. Our study has several implications. First, our work implicates airway neutrophil ferroptosis, which could be induced by IFN-ɣ, as a critical component of ARDS disease biology based on analyses of human, in vitro, and in vivo experimental data. We postulate that reactive oxygen species (ROS) and peroxidized lipids, markers of ferroptosis 65 that are damaging to the lung epithelium and endothelium, 66 , 67 contribute to ARDS disease severity. We note that ferroptosis could also induce neutrophil recruitment. 67 The increased abundance of neutrophils could subsequently elevate neutrophil extracellular traps (NETs) and proteases, 66 , 67 both damaging to lung tissues, 28 , 68 , 69 leading to more severe ARDS. Indeed, several cytokines involved in the recruitment of neutrophils, including IL-17 , 70 IL-1ɑ , 71 CCL3 , 72 and CXCL8 , 73 were upregulated in IFN-ɣ stimulated neutrophils. We note, however, that additional studies (e.g., using in vitro and in vivo models) are necessary to assess the relative impact of different neutrophil biologies on ARDS severity. Second, our work also creates avenues for testing ARDS therapeutic strategies; we showed that the ARDS mortality-associated airway neutrophil phenotype in human data could be induced in vitro via IFN-ɣ stimulation and in rodent models in vivo via flu/SARS-CoV-2 infection. Third, our analyses suggest that the mortality-associated neutrophil phenotype is likely specific to the lung: an association between CNA mortality-increasing gene signature and disease severity/influenza infection status was only detected in BALF and not blood, in both the analysis of the Bost et al. human data and the analysis of the mouse flu model data. Fourth, our study suggests that complex cytokine interactions could be involved in ARDS severity/mortality. Specific neutrophil functions (e.g., activation and degranulation) were much stronger when neutrophils were co-stimulated with IFN-ɣ and TNF-ɑ vs. stimulated with IFN-ɣ alone; distinct cytokine secretion profiles were also observed when neutrophils were stimulated with IFN-ɣ and TNF-ɑ together vs. IFN-ɣ alone. Future work investigating the synergistic effects of cytokines on ARDS is warranted. Finally, we discuss several limitations and caveats of this study. First, our all-patient analyses may suffer from reduced statistical power to detect mortality-associated disease biology, due to heterogeneity across COVID-19 and non-COVID-19 ARDS. Indeed, we observed different prioritizations of cell types (specifically between neutrophils and lymphocytes) in the all-patient vs. non-COVID-19 only analyses ( Figure 1 d , e ). We note that our choice to analyze all patients was in part driven by the limited number of patients, and that the neutrophil-specific mortality-associated disease biology was similar across COVID-19 and non-COVID-19 ARDS. Second, the relatively smaller number of non-COVID-19 vs. COVID-19 ARDS in our data could bias the identification of the mortality-associated genes and disease biology towards COVID-19 ARDS. We mitigated this imbalance by regressing out COVID-19 status from 30-day mortality in our all-patient analyses. However, unadjusted residual heterogeneity in ARDS disease biology for COVID-19 vs. non-COVID-19 patients may still bias our results. Third, we did not account for race and ethnicity, which have been shown to be associated with ARDS severity, 74 – 76 in our analyses, due to the limited number of sample sizes of specific minority populations in our study. Fourth, we also did not account for corticosteroids use, which has been shown to reduce ARDS mortality, 77 , 78 in our analyses of the COMET ETA scRNA-seq data. However, since the majority of the patients in the COMET cohort were treated with corticosteroids, adjusting for this as a covariate had minimal impact on the cell types and cell subpopulations identified. Fifth, our pathway analyses of CNA mortality-decreasing genes implicated TNF-ɑ ( Figure 2 d , Supplementary Table 4), previously reported to promote ferroptosis, 67 , 79 seemingly contradicting our conclusion that ferroptosis increases mortality. We note that the interaction between TNF-ɑ and IFN-ɣ more strongly induced ferroptosis than TNF-ɑ alone in our in vitro experiments. And it is possible that there was indeed stronger interaction between TNF-ɑ and IFN-ɣ, even though TNF-ɑ was comparably lower, in mortality cases vs. controls. Sixth, we note that while the effects of TNF-ɑ and IFN-ɣ on the CNA mortality-associated gene signatures were observed in our in vitro models, their effects remain uninvestigated in in vivo models, where the cytokine environments could be more complex. We also note that investigating the effects of other cytokines, e.g., IFN-ɑ and IL-6, etc., which were also implicated in our analysis ( Figure 2 c ), may provide additional biological insights into the CNA mortality-associated genes. Seventh, we note that the CNA mortality-decreasing neutrophil gene signature derived from COMET ETA scRNA-seq data was significantly positively associated with COVID-19 severity in the analysis of BALF neutrophils from Liao et al. 34 (Supplementary Figure 19d). We attribute this discrepancy to differences in the patient populations: relative expression of the gene signature could differ between severe vs. mild patients and mortality cases vs. controls among already severe patients. Eighth, we also note that replication was generally weaker for CNA mortality-decreasing vs. increasing gene signatures in the analysis of independent COVID-19 data ( Figure 3 ) or in vivo model data ( Figure 7 ). This suggests that mechanisms for survival/recovery could be more heterogeneous than those for mortality/deterioration, leading to overall weaker associations for the gene signature, although the sample sizes were also limited in these analyses. Ninth, we note that only a limited range of viral dosages were tested in the in vivo rodent models; and the resulting phenotype may better reflect ARDS severity instead of mortality. Tenth, we did not distinguish between subtypes of ARDS (e.g., hyperinflammatory vs. hypoinflammatory) in this work, other than COVID-19 vs. non-COVID-19 ARDS. Further reducing the heterogeneity in the ARDS patient populations could increase statistical power in detecting disease-associated biologies and lead to more effective treatments. 80 – 82 Eleventh, we note that various forms of cell death (e.g., ferroptosis, pyroptosis, and necrosis) are interconnected and involve shared biological processes, e.g., production of ROS. 67 Although our pathway analyses highlighted ferroptosis, other forms of cell death may also play critical roles in ARDS. 67 Twelfth, we have not investigated the interactions between mortality-associated neutrophils with other cell types. Understanding the cross talk between these neutrophils and other cell types (e.g., epithelial and endothelial cells) may further facilitate the development of therapeutics for ARDS. 83 , 84 Despite the limitations and caveats mentioned here, our study provides additional mechanistic insights into the role of airway neutrophils in ARDS mortality. Methods Overview of the COMET cohort The COVID-19 Multi-Phenotyping for Effective Therapies (COMET) collected samples longitudinally from hospitalized patients presenting with COVID-19 symptoms, further described in ref. 25 . Clinical severity metrics include 30 day all-cause mortality and Sequential Organ Failure Assessment (SOFA) scores that were assessed at study enrollment (Day 0), when baseline samples were collected. The COMET study was approved by the Institutional Review Board: UCSF Human Research Protection Program (HRPP) IRB# 20-30497 and informed consent was obtained for all patients. COMET endotracheal aspirate (ETA) scRNA-seq data The COMET ETA scRNA-seq data was collected as described in ref. 23 (see Table 1 for the patient characteristics of the data). We started from the raw read counts from this dataset, removing all cells with fewer than 200 total reads or 200 detected genes. We also removed cells with percent mitochondrial reads greater than 20%. After the QC step, 304K cells remained. We then log-normalized the raw read count data using the scater package, and performed PCA using the 2,000 most variable genes. We corrected for batch effects using mnnCorrect, 85 a batch effect correlation method based on mutual nearest neighbor (MNN). We performed shared nearest neighbor (SNN) based clustering using the 25 MNN dimensions; specifically, we used the specFindClusters function from Seurat v4.3.0.1, 86 with the R parameter set to 2.0 and all other parameters set as the default. We manually annotated the cell types based on the expression of the cell type marker genes (Supplementary Figure 1). Prior to our main analyses, we subset for the 129K cells corresponding to the earliest sample for each patient. Identifying cell subpopulations associated with 30-day mortality We performed co-varying neighborhood analysis 36 (CNA) to identify cell subpopulations associated with 30-day mortality. Briefly, CNA determines that a cell neighborhood, defined as a group of cells with similar transcriptomics profile (e.g., within some random walks on the k-nearest-neighbor graph of the scRNA-seq data) is associated with severity, if the cell neighborhood is enriched for cells from more severe patients. Since CNA defines a cell neighborhood anchored at each cell, it enables the identification of severity-associated cell subpopulations at individual cell resolution. In detail, we constructed k-nearest-neighborhood (k-NN) graphs for the cells in the scRNA-seq data based on the expression profiles of the top 2,000 highly variable genes (HVGs). We used the default settings of CNA to construct the neighborhood abundance matrices based on the k-NN graphs, and adjusted for appropriate covariates including age, sex at birth, COVID-19 test status, and batch. We did not adjust for corticosteroid use in our primary analyses, as this information was not available for 4 (out of 40) patients ( Table 1 a ). However, we note that analyses with and without adjusting for corticosteroid use yielded similar results (see Results). We performed CNA at both all cell-type levels, identifying the most critical broad cell types for ARDS severity, and at neutrophil-specific level, identifying neutrophil subpopulations associated with ARDS. We also performed CNA analyses with all (both COVID-19 and non-COVID-1) patients included and analyses restricting to non-COVID-19 patients. To identify genes indicative of a cell neighborhood’s association with mortality in neutrophils, we correlated the expression of each gene with the CNA neighborhood correlations across the cells. We defined CNA mortality-increasing and decreasing genes as the set of genes with R > 0.1 and R < 0.01, with two-tailed permutation based FDR < 0.05. We performed gene set enrichment analysis using the CNA mortality-associated genes using GSEApy 38 to identify the biological processes enriched in mortality-associated cell populations. Leave-one-out CNA analysis In the leave one-out framework, CNA was rerun as described above, but with one patient removed, and a CNA mortality-increasing signature and CNA mortality-decreasing signature derived through correlation of gene expression with the 30-day mortality CNA association statistic (FDR 0.1 for mortality-increasing, Pearson’s R < −0.1 for mortality-decreasing). The AUCell scores for each signature were then calculated for all patients. Once this process was repeated through all patients, the median leave-one-out AUCell score across all derived signatures was reported in Supplementary Figure 9. Inferring gene programs in ETA neutrophils We inferred gene programs for ETA neutrophils using cNMF. 37 Briefly, cNMF infers the sets of genes that collectively contribute to specific biological functions or cellular processes, and quantifies their overall activity (i.e., usages of the gene programs) in each individual cell based on non-negative factorization of the gene expression matrix. cNMF also requires a pre-specified number of factors, k, as the input. We applied cNMF to the count data from our preprocessed ETA neutrophil scRNA-seq data (see subsection below), which included 45,572 neutrophils (baseline samples) from all 38 patients, using the top 2,000 highly variable genes to infer the non-negative factor. We determined the optimal value for k, by incrementally testing k from 2 to 10. For each k, we ran cNMF with 30 random starts with a maximum of 1,000 iterations per random start. We elected to set k to 4, as it yielded moderate reconstruction error of the gene expression matrix with high stability. We applied cNMF to all patients to infer a consistent set of gene programs across COVID-19 and non-COVID-19 patients. To glean biological insights into each gene program, we performed gene set overrepresentation analysis of the top 100 genes with the highest spectra scores of each gene program using GSEApy. 38 Differential gene expression analysis of 30-day mortality using pseudobulked ETA neutrophils We performed differential expression analysis of 30-day mortality using pseudobulked ETA neutrophil gene expression. We restricted our analyses to the 31 patients with a minimum of 100 neutrophils in the samples of the earliest timepoint. For each patient, we obtained the pseudobulk expression as the sum of the raw read counts. We then used the PyDESeq2 package 87 to perform normalization and differential expression analysis of the pseudobulked gene expression. In the analysis of 30-day mortality of all patients (N fatality =7, N survivor =24), we adjusted for age, sex, COVID-19 status, and timepoint as the covariates; in the analysis of non-COVID-19 patients (N fatality =2, N survivor =8), we adjusted for age, sex, and timepoint as covariates. Comparing IFN-γ and TNF-α protein expression vs. the expression of CNA mortality-associated genes We obtained ETA protein expression data for a panel of 2,944 proteins from 17 patients using the Olink technology; 47 15 of these patients also had corresponding ETA scRNA-seq data (see Supplementary Table 6 a ). We obtained normalized protein expression (NPX) from the raw data using the approach as described in ref. 88 . We retained samples with more than 75% proteins above limit of detection (LOD) and with 2,000 proteins passing Olink’s default QC. We calculated aggregated scores for CNA mortality-increasing and decreasing genes, respectively, based on the log-normalized pseudobulked ETA neutrophil gene expression, using ssGSEA 44 for each of the 15 patients with both ETA protein expression and neutrophil scRNA-seq data. We compared the aggregated ssGSEA scores of CNA mortality-increasing and decreasing genes with the NPX of IFN-γ and IFN-α, respectively; IFN-γ and IFN-α signaling pathways were implicated in the gene set enrichment analysis of CNA mortality-increasing and decreasing genes, respectively (Results). We assessed the association between the ssGSEA scores and NPX both through Spearman’s rank correlation, and through linear models adjusting for age and sex. Mouse influenza A model lung and blood scRNA-seq data We obtained lung and blood scRNA-seq data from influenza A (PR8) infected mice with varied viral doses (10 3 and 10 4 PFU) at day 3 and 7 post treatment; cells from 3 mice across both the 2 tissues and time points were pooled into the same pool using HTO multiplexing. A full table of HTO ids and corresponding samples is provided in Supplementary Table 15. All 6 blood sample pools were run through the same sequencing run, and all BALF samples were run together in a separate run. Reads were aligned and processed with the CellRanger v6.1.1 pipeline and bcl2fastq2-v2.20.0.422., with GRCm38 used as a reference and with ENSEMBL version 90 annotations. Ambient RNA correction was applied to the data with CellBender (v0.2.0) 89 with the following settings: FPR = 0.01, epochs = 80, expected cells = 37k, expected droplets = 52,000 for each sequencing pool. Following this, HTO demultiplexing was performed with HTODemux (Seurat v4.3.0.1); and droplets were removed if assigned as a doublet or if the difference between the expression score for the first and second assignment was less than 0.5 (i.e., an ambiguous droplet). Cells were then filtered so that each cell had at least 500 UMI counts and 400 detected genes, and with mitochondrial proportions below 3 median absolute deviation above median. scDblFinder (v1.14.0) 90 was then used with the following parameters: dbr.sd = 1.0 (renders prior estimated doublet rate meaningless), nfeatures = 2500. scDblFinder was run separately for each sequencing run. Cell types were then labeled by marker gene expression (Supplementary Figure 21). We then calculated AUCell 45 scores for each neutrophils using the mouse orthologs of the CNA mortality-associated genes, and performed t-tests to compare tissue/dose/timepoint conditions, particularly between 10 3 PFU and 10 4 PFU for a fixed tissue and timepoint. in vitro human blood neutrophil bulk RNA-seq data Neutrophils were isolated from peripheral blood of healthy donors (n=3) by negative selection via the StemCell Neutrophil Isolation Kit. These neutrophils were stimulated per 10 6 cells with vehicle, TNF-ɑ (0.1μg/ml), LPS (0.1μg/ml), IFN-ɣ (0.1-1.0μg/ml) and combinations of IFN-ɣ with TNF-ɑ or LPS for 4 hours at 37 °C. All stimulations were performed in technical triplicates for each donor. RNA was isolated using the Qiagen RNeasy mini kit with additional DNase digestion. RNA was pooled from technical triplicate wells for each of 3 biological replicates (human donors). RNA purity and concentrations were analyzed via NaNoDrop. Bulk RNAseq was performed on all pooled samples. For each donor, RNA was isolated from 10 6 freshly isolated neutrophils to compare changes in gene expression following stimulation, since pilot experiments show that gene expression and cell phenotypes are altered after culture at 37 °C. RNA-seq data was processed as follows. First, reads with low nucleotide qualities (70% of bases with quality less than 23) or matches to rRNA and adapter sequences were removed using the HTSeqGenie software. 91 Remaining sequencing reads were mapped to the reference human genome (GRCh38), using the GSNAP short read aligner. 92 Expression counts per gene were quantified using HTSeqGenie, 91 using ENSEMBL version 90 features. Neutrophil flow cytometry activation and degranulation analysis Human peripheral neutrophils were isolated from healthy donors (n=6) as mentioned above. For each donor, 0.25*10 6 neutrophils were stimulated with vehicles, TNF-ɑ (0.1μg/ml), LPS (0.1μg/ml), IFN-ɣ (0.1-1.0μg/ml) and combinations of IFN-ɣ with TNF-ɑ or LPS. After 4 hours of incubation, neutrophils were labeled with near-infrared live-dead dye (Invitrogen), SuperBright600 anti-human CD14, NovaFluorBlue 610-30s anti-human CD16, and BrilliantViolet711 anti-human CD62L, PacificBlue anti-human CD63. This antibody panel was verified with the EasyPanel software. Viable neutrophils were gated by subsequent gating for single cells (forward scatter area/forward scatter height), viable cells (live-dead staining). Activated neutrophils were expressed as a percentage of CD62L - and CD14 - cells, degranulating neutrophils are expressed as percentage for CD63 + and CD66b + cells. Luminex immunoassay for measuring cytokine release by stimulated neutrophils A multiplexed Luminex immunoassay was used to map the cytokine release profile of the different neutrophil subpopulations. Human peripheral neutrophils were isolated from healthy donors (n=6) as mentioned above. For each donor, 2.5×10 5 neutrophils were stimulated with vehicles, TNF-ɑ (0.1μg/ml), LPS (0.1μg/ml), IFN-ɣ (0.1-1.0μg/ml) and combinations of IFN-ɣ with TNF-ɑ or LPS. After 4 hours of incubation, supernatants were taken for cytokine measurement via Millipore Human 30-plex Luminex. Concentrations below detection limit were considered negative and were given a value one-half of the detection limit. Neutrophil elastase activity assay Human peripheral neutrophils were isolated from healthy donors (n=6) as mentioned above. For each donor, 2.5×10 5 neutrophils were stimulated with vehicles, TNF-ɑ, LPS, IFN-ɣ or combinations of these as described in the previous section. A fluorometric neutrophil elastase activity assay kit (Abcam) was used as instructed by the manufacturers. Duplicates and background controls were included for all samples in both assays. Values were normalized to vehicle conditions for each donor. Lipid peroxidation and viability assays to analyze ferroptosis in vitro Human peripheral neutrophils were isolated from healthy donors (n=6) as mentioned above. For each donor, 2.5×10 5 neutrophils were stimulated with vehicles, TNF-ɑ (0.1 μg/ml) or IFN-ɣ (0.1μg/ml). Erastin (2μM), a known ferroptosis inducer, was added as positive control. Neutrophils were stimulated with each of these stimuli with and without the addition of ferrostatin-1 (1μM). Ferrostatin-1 is a synthetic anti-oxidant that inhibits ferroptosis in a two fold manner. It scavenging the alkoxyl radicals of oxidized lipids while oxidizing ferrous iron back to its less reactive form ferric iron. At the end of 4 hours incubation, lipid peroxidation was analysed through oxidation of BODIPY 581/591 C11 reagent via the Image-iT™ Lipid Peroxidation Kit (Thermofisher). The ratio of lipid peroxidation was calculated as the ratio of oxidized lipids (MFI FITC) to both oxidized and reduced lipids (MFI Texas Red), analysed via flow cytometry. Values of lipid peroxidation were normalized to the vehicle conditions for each donor. Neutrophil viability was determined after 24h in culture with the above mentioned stimuli via PI staining. The number of viable cells was analysed via automatic counting chamber and the ratio of neutrophil viability was calculated after normalizing to vehicle conditions for each donor. Bost et al. COVID-19 bronchoalveolar lavage fluid (BALF) and blood scRNA-seq data We obtained publicly available BALF and blood COVID-19 scRNA-seq data from 21 severe COVID-19 patients (N fatility =8, N survivor =12 for BALF; N fatility =7, N survivor =13 for blood; Supplementary Table 9) through Gene Expression Omnibus (GSE number: GSE157344) from ref. 48 . We performed log-normalization, PCA, and clustering 93 on the raw read count data. We manually annotated the cells using lung cell type marker genes (Supplementary Figure 15). We then subset for the neutrophils, and obtained pseudobulked gene expression for patients with more than 25 neutrophils, yielding 20 patients for blood and BALF samples in total (19 overlapping patients between blood and BALF). We performed ssGSEA scoring 44 (gsva v1.50.0) 94 on the pseudobulk gene expression data using the CNA mortality-increasing and decreasing genes, and assessed association between SOFA scores and ssGSEA scores through Spearman’s correlation (p-values obtained via permutation) and through modeling SOFA scores as a function of ssGSEA scores in linear models adjusting for age or age and sex (p-values obtained analytically). We also assessed association between the ssGSEA scores and death through comparing the mean ssGSEA scores of survived vs. deceased patients (p-values obtained via permutation) and through modeling mortality status as a function of the the ssGSEA scores in logistic regression models adjusting for age or for age and sex as covariates (p-values obtained analytically). Liao et al. COVID-19 BALF scRNA-seq data We obtained publicly available BALF COVID-19 scRNA-seq data from 9 donors (5 severe, 3 mild, and 1 healthy control) through Gene Expression Omnibus (GSE number: GSE145926) from ref. 34 . We started from the raw read count data, filtering out cells that had fewer than 500 UMI counts or fewer than 200 detected genes, and setting upper thresholds for mitochondrial read fraction (above which droplets were excluded) at the minimum of 3 median absolute deviations or 0.25 (Supplementary Figure 18a–c). Once this initial filtering was done, we applied scDblFinder (v1.14.0) with the following parameters: dbr =NULL (i.e. doublet rate assumed to be 1% doublets for every 1,000 droplets in a sample), dbr.sd =0.0, nfeatures =2500, and batch set to sample ID (i.e., run scDblFinder fitting model for each batch separately) (Supplementary Figure 18d). After the QC and doublet filtering, 65,187 singlets remained. We then log normalized the count data, ran PCA, corrected PCA coordinates via mutual nearest neighbors (MNN), 85 and performed clustering based on a shared nearest graph in the MNN corrected space. 93 We calculated AUCell scores (AUCell v 1.2.4.0) 45 for a set of lung cell type markers (Supplementary Figure 18f), and examined the expression of FCGR3B . We annotated cluster 13 as neutrophils. We then calculated AUCell scores for each of the neutrophils 13 using CNA mortality-associated genes. We compared the mean AUCell scores across different severity, and obtained test statistics using permutation. Peidli et al. hamster SARS-CoV-9 model BALF scRNA-seq data We obtained publicly available hamster SARS-CoV-9 model BALF scRNA-seq data from 12 hamsters across multiple infection dosages and time points through Gene Expression Omnibus (GSE number: GSE241133) from ref. 64 . We subset for neutrophils based on the cell type labels provided by the authors, and obtained pseudobulk neutrophil gene expression of the dataset, aggregating by tissue, dose, timepoint, and replicate (corresponding to one experimental run per set of aggregated cells). We performed ssGSEA 44 scoring for the CNA mortality-increasing and decreasing genes on the pseudobulk data using the gsva (v1.50.0) package. 94 We compared the mean ssGSEA scores across different infection dosages for each timepoint, and obtained test statistics using two-sided t-tests. Data Availability All data produced in the present study are available upon reasonable request to the authors Author contributions O.W., T.D., C.R., and H.S. designed the experiments. T.D. performed the experiments. O.W. and H.S. analyzed the data. O.W., T.D., and H.S. wrote the paper with assistance from all authors. Competing interests O.W., T.D., M.R., K.K, H.L., C.C, M.D.S., M.X., H.S., J.H., D.C., X.G., E.K., N.R., S.B.K., C.M.R. and H.S. are employees of Genentech, Inc. at the time of this study and own stocks in Roche. Data and code availability The Bost et al. COVID-19 bronchoalveolar lavage fluid (BALF) and blood scRNA-seq data is available on GEO with accession number GSE157344. The Liao et al. COVID-19 BALF scRNA-seq data is available on GEO with accession number GSE145926. data The Peidli et al. hamster SARS-CoV-9 infection model BALF scRNA-seq data is available on GEO with accession number GSE241133. Analysis scripts and additional experimental data generated in this study will be available upon reasonable request to the authors. Acknowledgements The authors are grateful to K. Kiani, Y. Yang, W. Mu, and L. Orozco for helpful discussions. Footnotes Author list updated to include co-author Rafal K. Woycicki. Author affiliation updated to include the affiliation for Rafal K. Woycicki: 9 Department of Data & Analytics Chapter, Computational Science, Roche Poland, 02-679 Warszawa, Poland References 1. ↵ Matthay , M.A. , Zemans , R.L. , Zimmerman , G.A. , Arabi , Y.M. , Beitler , J.R. , Mercat , A. , Herridge , M. , Randolph , A.G. , and Calfee , C.S . ( 2019 ). Acute respiratory distress syndrome . Nature Reviews Disease Primers 5 , 1 – 22 . OpenUrl CrossRef PubMed 2. Ashbaugh , D.G. , Bigelow , D.B. , Petty , T.L. , and Levine , B.E . ( 1967 ). Acute respiratory distress in adults . Lancet 2 , 319 – 323 . OpenUrl CrossRef PubMed Web of Science 3. ↵ ARDS Definition Task Force , Ranieri , V.M. , Rubenfeld , G.D. , Thompson , B.T. , Ferguson , N.D. , Caldwell , E. , Fan , E. , Camporota , L. , and Slutsky , A.S. ( 2012 ). Acute respiratory distress syndrome: the Berlin Definition . JAMA 307 , 2526 – 2533 . OpenUrl CrossRef PubMed Web of Science 4. ↵ Swenson , K.E. , and Swenson , E.R . ( 2021 ). Pathophysiology of Acute Respiratory Distress Syndrome and COVID-19 Lung Injury . Crit Care Clin 37 , 749 – 776 . OpenUrl CrossRef PubMed 5. ↵ Bellani , G. , Laffey , J.G. , Pham , T. , Fan , E. , Brochard , L. , Esteban , A. , Gattinoni , L. , van Haren , F. , Larsson , A. , McAuley , D.F. , et al. ( 2016 ). Epidemiology, Patterns of Care, and Mortality for Patients With Acute Respiratory Distress Syndrome in Intensive Care Units in 50 Countries . JAMA 315 , 788 – 800 . OpenUrl CrossRef PubMed 6. ↵ Matthay , M.A. , and Zemans , R.L . ( 2011 ). The acute respiratory distress syndrome: pathogenesis and treatment . Annu Rev Pathol 6 , 147 – 163 . OpenUrl CrossRef PubMed 7. ↵ Diamond , M. , Peniston , H.L. , Sanghavi , D.K. , and Mahapatra , S. ( 2024 ). Acute Respiratory Distress Syndrome . In StatPearls ( StatPearls Publishing ). 8. ↵ Ruhl , A.P. , Huang , M. , Colantuoni , E. , Karmarkar , T. , Dinglas , V.D. , Hopkins , R.O. , Needham , D.M. , and With the National Institutes of Health, National Heart, Lung, and Blood Institute Acute Respiratory Distress Syndrome Network ( 2017 ). Healthcare utilization and costs in ARDS survivors: a 1-year longitudinal national US multicenter study . Intensive Care Med 43 , 980 – 991 . OpenUrl PubMed 9. ↵ Robles , A.J. , Kornblith , L.Z. , Hendrickson , C.M. , Howard , B.M. , Conroy , A.S. , Moazed , F. , Calfee , C.S. , Cohen , M.J. , and Callcut , R.A . ( 2018 ). Health care utilization and the cost of posttraumatic acute respiratory distress syndrome care . J Trauma Acute Care Surg 85 , 148 – 154 . OpenUrl PubMed 10. ↵ Acute Respiratory Distress Syndrome Network , Brower , R.G. , Matthay , M.A. , Morris , A. , Schoenfeld , D. , Thompson , B.T. , and Wheeler , A. ( 2000 ). Ventilation with lower tidal volumes as compared with traditional tidal volumes for acute lung injury and the acute respiratory distress syndrome . N Engl J Med 342 , 1301 – 1308 . OpenUrl CrossRef PubMed Web of Science 11. ↵ Combes , A. , Hajage , D. , Capellier , G. , Demoule , A. , Lavoué , S. , Guervilly , C. , Da Silva , D. , Zafrani , L. , Tirot , P. , Veber , B. , et al. ( 2018 ). Extracorporeal Membrane Oxygenation for Severe Acute Respiratory Distress Syndrome . N Engl J Med 378 , 1965 – 1975 . OpenUrl CrossRef PubMed 12. ↵ Qadir , N. , Sahetya , S. , Munshi , L. , Summers , C. , Abrams , D. , Beitler , J. , Bellani , G. , Brower , R.G. , Burry , L. , Chen , J.-T. , et al. ( 2024 ). An Update on Management of Adult Patients with Acute Respiratory Distress Syndrome: An Official American Thoracic Society Clinical Practice Guideline . Am J Respir Crit Care Med 209 , 24 – 36 . OpenUrl CrossRef PubMed 13. ↵ RECOVERY Collaborative Group , Horby , P. , Lim , W.S. , Emberson , J.R. , Mafham , M. , Bell , J.L. , Linsell , L. , Staplin , N. , Brightling , C. , Ustianowski , A. , et al. ( 2021 ). Dexamethasone in Hospitalized Patients with Covid-19 . N Engl J Med 384 , 693 – 704 . OpenUrl CrossRef PubMed 14. Khilnani , G.C. , and Hadda , V . ( 2011 ). Corticosteroids and ARDS: A review of treatment and prevention evidence . Lung India 28 , 114 – 119 . OpenUrl CrossRef PubMed 15. ↵ Boyle , A.J. , Mac Sweeney , R. , and McAuley , D.F . ( 2013 ). Pharmacological treatments in ARDS; a state-of-the-art update . BMC Med 11 , 166 . OpenUrl CrossRef PubMed 16. ↵ Han , S. , and Mallampalli , R.K . ( 2015 ). The acute respiratory distress syndrome: from mechanism to translation . J Immunol 194 , 855 – 860 . OpenUrl Abstract / FREE Full Text 17. ↵ Donahoe , M . ( 2011 ). Acute respiratory distress syndrome: A clinical review . Pulm Circ 1 , 192 – 211 . OpenUrl CrossRef PubMed 18. ↵ Wick , K.D. , McAuley , D.F. , Levitt , J.E. , Beitler , J.R. , Annane , D. , Riviello , E.D. , Calfee , C.S. , and Matthay , M.A . ( 2021 ). Promises and challenges of personalized medicine to guide ARDS therapy . Crit Care 25 , 404 . OpenUrl CrossRef PubMed 19. ↵ Hu , B. , Huang , S. , and Yin , L . ( 2021 ). The cytokine storm and COVID-19 . J Med Virol 93 , 250 – 256 . OpenUrl CrossRef PubMed 20. ↵ D’Alessio , F.R. , and Heller , N.M . ( 2020 ). COVID-19 and myeloid cells: complex interplay correlates with lung severity . J Clin Invest 130 , 6214 – 6217 . OpenUrl CrossRef PubMed 21. Jiang , Y. , Rosborough , B.R. , Chen , J. , Das , S. , Kitsios , G.D. , McVerry , B.J. , Mallampalli , R.K. , Lee , J.S. , Ray , A. , Chen , W. , et al. ( 2020 ). Single cell RNA sequencing identifies an early monocyte gene signature in acute respiratory distress syndrome . JCI Insight 5 . doi: 10.1172/jci.insight.135678 . OpenUrl CrossRef 22. ↵ Dang , W. , Tao , Y. , Xu , X. , Zhao , H. , Zou , L. , and Li , Y . ( 2022 ). The role of lung macrophages in acute respiratory distress syndrome . Inflamm Res 71 , 1417 – 1432 . OpenUrl CrossRef PubMed 23. ↵ Sarma , A. , Christenson , S.A. , Byrne , A. , Mick , E. , Pisco , A.O. , DeVoe , C. , Deiss , T. , Ghale , R. , Zha , B.S. , Tsitsiklis , A. , et al. ( 2021 ). Tracheal aspirate RNA sequencing identifies distinct immunological features of COVID-19 ARDS . Nat Commun 12 , 5152 . OpenUrl CrossRef PubMed 24. ↵ Schulte-Schrepping , J. , Reusch , N. , Paclik , D. , Baßler , K. , Schlickeiser , S. , Zhang , B. , Krämer , B. , Krammer , T. , Brumhard , S. , Bonaguro , L. , et al. ( 2020 ). Severe COVID-19 Is Marked by a Dysregulated Myeloid Cell Compartment . Cell 182 , 1419 – 1440.e23 . OpenUrl CrossRef PubMed 25. ↵ Hackney , J.A. , Shivram , H. , Vander Heiden , J. , Overall , C. , Orozco , L. , Gao , X. , Kim , E. , West , N. , Qamra , A. , Chang , D. , et al. ( 2023 ). A myeloid program associated with COVID-19 severity is decreased by therapeutic blockade of IL-6 signaling . iScience 26 , 107813 . OpenUrl PubMed 26. ↵ Reyes , M. , Filbin , M.R. , Bhattacharyya , R.P. , Sonny , A. , Mehta , A. , Billman , K. , Kays , K.R. , Pinilla-Vera , M. , Benson , M.E. , Cosimi , L.A. , et al. ( 2021 ). Plasma from patients with bacterial sepsis or severe COVID-19 induces suppressive myeloid cell production from hematopoietic progenitors in vitro . Science Translational Medicine . doi: 10.1126/scitranslmed.abe9599 . OpenUrl Abstract / FREE Full Text 27. ↵ Williams , A.E. , and Chambers , R.C . ( 2014 ). The mercurial nature of neutrophils: still an enigma in ARDS? Am J Physiol Lung Cell Mol Physiol 306 , L217 – 30 . OpenUrl CrossRef PubMed Web of Science 28. ↵ Yang , S.-C. , Tsai , Y.-F. , Pan , Y.-L. , and Hwang , T.-L . ( 2021 ). Understanding the role of neutrophils in acute respiratory distress syndrome . Biomed J 44 , 439 – 446 . OpenUrl PubMed 29. Juss , J.K. , House , D. , Amour , A. , Begg , M. , Herre , J. , Storisteanu , D.M.L. , Hoenderdos , K. , Bradley , G. , Lennon , M. , Summers , C. , et al. ( 2016 ). Acute Respiratory Distress Syndrome Neutrophils Have a Distinct Phenotype and Are Resistant to Phosphoinositide 3-Kinase Inhibition . Am J Respir Crit Care Med 194 , 961 – 973 . OpenUrl CrossRef PubMed 30. ↵ Hsieh , L.L. , Thompson , E.A. , Jairam , N.P. , Roznik , K. , Figueroa , A. , Aytenfisu , T. , Zhou , W. , Gour , N. , Chao , K.-H. , Milstone , A.M. , et al. ( 2025 ). SARS-CoV-2 induces neutrophil degranulation and differentiation into myeloid-derived suppressor cells associated with severe COVID-19 . Science Translational Medicine . doi: 10.1126/scitranslmed.adn7527 . OpenUrl CrossRef 31. ↵ Ng , L.G. , Ostuni , R. , and Hidalgo , A . ( 2019 ). Heterogeneity of neutrophils . Nat Rev Immunol 19 , 255 – 265 . OpenUrl CrossRef PubMed 32. Wigerblad , G. , Cao , Q. , Brooks , S. , Naz , F. , Gadkari , M. , Jiang , K. , Gupta , S. , O’Neil , L. , Dell’Orso , S. , Kaplan , M.J. , et al. ( 2022 ). Single-Cell Analysis Reveals the Range of Transcriptional States of Circulating Human Neutrophils . J Immunol 209 , 772 – 782 . OpenUrl Abstract / FREE Full Text 33. ↵ Kuhns , D.B. , Priel , D.A.L. , Chu , J. , and Zarember , K.A . ( 2015 ). Isolation and Functional Analysis of Human Neutrophils . Curr Protoc Immunol 111 , 7 .23.1–7.23.16. OpenUrl 34. ↵ Liao , M. , Liu , Y. , Yuan , J. , Wen , Y. , Xu , G. , Zhao , J. , Cheng , L. , Li , J. , Wang , X. , Wang , F. , et al. ( 2020 ). Single-cell landscape of bronchoalveolar immune cells in patients with COVID-19 . Nat. Med . 26 , 842 – 844 . OpenUrl CrossRef PubMed 35. ↵ Yoshida , M. , Worlock , K.B. , Huang , N. , Lindeboom , R.G.H. , Butler , C.R. , Kumasaka , N. , Dominguez Conde , C. , Mamanova , L. , Bolt , L. , Richardson , L. , et al. ( 2022 ). Local and systemic responses to SARS-CoV-2 infection in children and adults . Nature 602 , 321 – 327 . OpenUrl CrossRef PubMed 36. ↵ Reshef , Y.A. , Rumker , L. , Kang , J.B. , Nathan , A. , Korsunsky , I. , Asgari , S. , Murray , M.B. , Moody , D.B. , and Raychaudhuri , S . ( 2021 ). Co-varying neighborhood analysis identifies cell populations associated with phenotypes of interest from single-cell transcriptomics . Nat. Biotechnol . 40 , 355 – 363 . OpenUrl PubMed 37. ↵ Kotliar , D. , Veres , A. , Aurel Nagy , M. , Tabrizi , S. , Hodis , E. , Melton , D.A. , and Sabeti , P.C. ( 2019 ). Identifying gene expression programs of cell-type identity and cellular activity with single-cell RNA-Seq . doi: 10.7554/eLife.43803 . OpenUrl CrossRef PubMed 38. ↵ Fang , Z. , Liu , X. , and Peltz , G . ( 2022 ). GSEApy: a comprehensive package for performing gene set enrichment analysis in Python . Bioinformatics 39 , btac757 . OpenUrl 39. ↵ Gomez , J.C. , Yamada , M. , Martin , J.R. , Dang , H. , Brickey , W.J. , Bergmeier , W. , Dinauer , M.C. , and Doerschuk , C.M . ( 2015 ). Mechanisms of interferon-γ production by neutrophils and its function during Streptococcus pneumoniae pneumonia . Am J Respir Cell Mol Biol 52 , 349 – 364 . OpenUrl CrossRef PubMed 40. ↵ Darwich , L. , Coma , G. , Peña , R. , Bellido , R. , Blanco , E.J.J. , Este , J.A. , Borras , F.E. , Clotet , B. , Ruiz , L. , Rosell , A. , et al. ( 2009 ). Secretion of interferon-gamma by human macrophages demonstrated at the single-cell level after costimulation with interleukin (IL)-12 plus IL-18 . Immunology 126 , 386 – 393 . OpenUrl CrossRef PubMed Web of Science 41. ↵ Feiken , E. , Rømer , J. , Eriksen , J. , and Lund , L.R . ( 1995 ). Neutrophils express tumor necrosis factor-alpha during mouse skin wound healing . J Invest Dermatol 105 , 120 – 123 . OpenUrl CrossRef PubMed Web of Science 42. ↵ Parameswaran , N. , and Patial , S . ( 2010 ). Tumor necrosis factor-α signaling in macrophages . Crit Rev Eukaryot Gene Expr 20 , 87 – 103 . OpenUrl CrossRef PubMed 43. ↵ Krug , K. , Mertins , P. , Zhang , B. , Hornbeck , P. , Raju , R. , Ahmad , R. , Szucs , M. , Mundt , F. , Forestier , D. , Jane-Valbuena , J. , et al. ( 2019 ). A Curated Resource for Phosphosite-specific Signature Analysis . Mol Cell Proteomics 18 , 576 – 593 . OpenUrl Abstract / FREE Full Text 44. ↵ Barbie , D.A. , Tamayo , P. , Boehm , J.S. , Kim , S.Y. , Moody , S.E. , Dunn , I.F. , Schinzel , A.C. , Sandy , P. , Meylan , E. , Scholl , C. , et al. ( 2009 ). Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1 . Nature 462 , 108 – 112 . OpenUrl CrossRef PubMed Web of Science 45. ↵ Aibar , S. , González-Blas , C.B. , Moerman , T. , Huynh-Thu , V.A. , Imrichova , H. , Hulselmans , G. , Rambow , F. , Marine , J.-C. , Geurts , P. , Aerts , J. , et al. ( 2017 ). SCENIC: single-cell regulatory network inference and clustering . Nat. Methods 14 , 1083 – 1086 . OpenUrl CrossRef PubMed 46. ↵ Love , M.I. , Huber , W. , and Anders , S . ( 2014 ). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 . Genome Biol . 15 , 1 – 21 . OpenUrl CrossRef PubMed 47. ↵ Wik , L. , Nordberg , N. , Broberg , J. , Björkesten , J. , Assarsson , E. , Henriksson , S. , Grundberg , I. , Pettersson , E. , Westerberg , C. , Liljeroth , E. , et al. ( 2021 ). Proximity Extension Assay in Combination with Next-Generation Sequencing for High-throughput Proteome-wide Analysis . Molecular & Cellular Proteomics: MCP 20 , 100168 . OpenUrl PubMed 48. ↵ Bost , P. , De Sanctis , F. , Canè , S. , Ugel , S. , Donadello , K. , Castellucci , M. , Eyal , D. , Fiore , A. , Anselmi , C. , Barouni , R.M. , et al. ( 2021 ). Deciphering the state of immune silence in fatal COVID-19 patients . Nat Commun 12 , 1428 . OpenUrl CrossRef PubMed 49. ↵ Buras , J.A. , Holzmann , B. , and Sitkovsky , M . ( 2005 ). Animal models of sepsis: setting the stage . Nat Rev Drug Discov 4 , 854 – 865 . OpenUrl CrossRef PubMed Web of Science 50. ↵ Lacy , P . ( 2006 ). Mechanisms of degranulation in neutrophils . Allergy Asthma Clin Immunol 2 , 98 – 108 . OpenUrl CrossRef PubMed 51. ↵ Wright , D.G . ( 1988 ). Human neutrophil degranulation . Methods Enzymol 162 , 538 – 551 . OpenUrl CrossRef PubMed Web of Science 52. ↵ Williams , A.E. , José , R.J. , Mercer , P.F. , Brealey , D. , Parekh , D. , Thickett , D.R. , O’Kane , C. , McAuley , D.F. , and Chambers , R.C . ( 2017 ). Evidence for chemokine synergy during neutrophil migration in ARDS . Thorax 72 , 66 – 73 . OpenUrl Abstract / FREE Full Text 53. Mikacenic , C. , Hansen , E.E. , Radella , F. , Gharib , S.A. , Stapleton , R.D. , and Wurfel , M.M . ( 2016 ). Interleukin-17A Is Associated With Alveolar Inflammation and Poor Outcomes in Acute Respiratory Distress Syndrome . Crit Care Med 44 , 496 – 502 . OpenUrl CrossRef PubMed 54. Aggarwal , B.B. , Sung , B. , and Gupta , S.C . ( 2014 ). Inflammation and Cancer (Springer ). 55. ↵ Sugiyama , M. , Kinoshita , N. , Ide , S. , Nomoto , H. , Nakamoto , T. , Saito , S. , Ishikane , M. , Kutsuna , S. , Hayakawa , K. , Hashimoto , M. , et al. ( 2021 ). Serum CCL17 level becomes a predictive marker to distinguish between mild/moderate and severe/critical disease in patients with COVID-19 . Gene 766 , 145145 . OpenUrl PubMed 56. ↵ Reed , A. , Ichu , T.-A. , Milosevich , N. , Melillo , B. , Schafroth , M.A. , Otsuka , Y. , Scampavia , L. , Spicer , T.P. , and Cravatt , B.F . ( 2022 ). LPCAT3 Inhibitors Remodel the Polyunsaturated Phospholipid Content of Human Cells and Protect from Ferroptosis . ACS Chem Biol 17 , 1607 – 1618 . OpenUrl CrossRef PubMed 57. Doll , S. , Proneth , B. , Tyurina , Y.Y. , Panzilius , E. , Kobayashi , S. , Ingold , I. , Irmler , M. , Beckers , J. , Aichler , M. , Walch , A. , et al. ( 2017 ). ACSL4 dictates ferroptosis sensitivity by shaping cellular lipid composition . Nat Chem Biol 13 , 91 – 98 . OpenUrl CrossRef PubMed 58. Ma , X.-H. , Liu , J.-H.-Z. , Liu , C.-Y. , Sun , W.-Y. , Duan , W.-J. , Wang , G. , Kurihara , H. , He , R.-R. , Li , Y.-F. , Chen , Y. , et al. ( 2022 ). ALOX15-launched PUFA-phospholipids peroxidation increases the susceptibility of ferroptosis in ischemia-induced myocardial damage . Signal Transduct Target Ther 7 , 288 . OpenUrl PubMed 59. Zhang , W. , Huang , F. , Ding , X. , Qin , J. , Wang , W. , and Luo , L . ( 2024 ). Identifying ALOX15-initiated lipid peroxidation increases susceptibility to ferroptosis in asthma epithelial cells . Biochim Biophys Acta Mol Basis Dis 1870 , 167176 . OpenUrl 60. Tang , D. , Chen , X. , Kang , R. , and Kroemer , G . ( 2021 ). Ferroptosis: molecular mechanisms and health implications . Cell Res 31 , 107 – 125 . OpenUrl CrossRef PubMed 61. ↵ Ou , Y. , Wang , S.-J. , Li , D. , Chu , B. , and Gu , W . ( 2016 ). Activation of SAT1 engages polyamine metabolism with p53-mediated ferroptotic responses . Proc Natl Acad Sci U S A 113 , E6806 – E6812 . OpenUrl Abstract / FREE Full Text 62. ↵ Zhou , Q. , Meng , Y. , Le , J. , Sun , Y. , Dian , Y. , Yao , L. , Xiong , Y. , Zeng , F. , Chen , X. , and Deng , G . ( 2024 ). Ferroptosis: mechanisms and therapeutic targets . MedComm (2020) 5 , e70010 . OpenUrl 63. ↵ Miotto , G. , Rossetto , M. , Di Paolo , M.L. , Orian , L. , Venerando , R. , Roveri , A. , Vučković , A.-M. , Bosello Travain , V. , Zaccarin , M. , Zennaro , L. , et al. ( 2020 ). Insight into the mechanism of ferroptosis inhibition by ferrostatin-1 . Redox Biol 28 , 101328 . OpenUrl PubMed 64. ↵ Peidli , S. , Nouailles , G. , Wyler , E. , Adler , J.M. , Kunder , S. , Voß , A. , Kazmierski , J. , Pott , F. , Pennitz , P. , Postmus , D. , et al. ( 2024 ). Single-cell-resolved interspecies comparison shows a shared inflammatory axis and a dominant neutrophil-endothelial program in severe COVID-19 . Cell Rep 43 , 114328 . OpenUrl PubMed 65. ↵ Li , J. , Cao , F. , Yin , H.-L. , Huang , Z.-J. , Lin , Z.-T. , Mao , N. , Sun , B. , and Wang , G . ( 2020 ). Ferroptosis: past, present and future . Cell Death Dis 11 , 88 . OpenUrl CrossRef PubMed 66. ↵ Qu , M. , Zhang , H. , Chen , Z. , Sun , X. , Zhu , S. , Nan , K. , Chen , W. , and Miao , C . ( 2021 ). The Role of Ferroptosis in Acute Respiratory Distress Syndrome . Front. Med . 8 , 651552 . OpenUrl 67. ↵ Zheng , Y. , Huang , Y. , Xu , Y. , Sang , L. , Liu , X. , and Li , Y . ( 2023 ). Ferroptosis, pyroptosis and necroptosis in acute respiratory distress syndrome . Cell Death Discov 9 , 91 . OpenUrl PubMed 68. ↵ Neutrophil extracellular traps induce persistent lung tissue damage via thromboinflammation without altering virus resolution in a mouse coronavirus model ( 2024 ). Journal of Thrombosis and Haemostasis 22 , 188 – 198 . OpenUrl 69. ↵ Lefrançais , E. , Mallavia , B. , Zhuo , H. , Calfee , C.S. , and Looney , M.R . ( 2018 ). Maladaptive role of neutrophil extracellular traps in pathogen-induced lung injury . JCI Insight 3 . doi: 10.1172/jci.insight.98178 . OpenUrl CrossRef PubMed 70. ↵ Zenobia , C. , and Hajishengallis , G . ( 2015 ). Basic biology and role of interleukin-17 in immunity and inflammation . Periodontol 2000 69 , 142 – 159 . OpenUrl 71. ↵ Malik , A. , and Kanneganti , T.-D . ( 2018 ). Function and regulation of IL-1α in inflammatory diseases and cancer . Immunol Rev 281 , 124 – 137 . OpenUrl CrossRef PubMed 72. ↵ Bonville , C.A. , Percopo , C.M. , Dyer , K.D. , Gao , J. , Prussin , C. , Foster , B. , Rosenberg , H.F. , and Domachowske , J.B . ( 2009 ). Interferon-gamma coordinates CCL3-mediated neutrophil recruitment in vivo . BMC Immunol 10 , 14 . OpenUrl CrossRef PubMed 73. ↵ Park , J.H. , and Lee , H.K . ( 2020 ). Re-analysis of Single Cell Transcriptome Reveals That the NR3C1-CXCL8-Neutrophil Axis Determines the Severity of COVID-19 . Front. Immunol . 11 , 573611 . OpenUrl 74. ↵ Price-Haywood , E.G. , Burton , J. , Fort , D. , and Seoane , L . ( 2020 ). Hospitalization and Mortality among Black Patients and White Patients with Covid-19 . N Engl J Med 382 , 2534 – 2543 . OpenUrl CrossRef PubMed 75. Yancy , C.W . ( 2020 ). COVID-19 and African Americans . JAMA 323 , 1891 – 1892 . OpenUrl CrossRef PubMed 76. ↵ Webb Hooper , M. , Nápoles , A.M. , and Pérez-Stable , E.J. ( 2020 ). COVID-19 and Racial/Ethnic Disparities . JAMA 323 , 2466 – 2467 . OpenUrl CrossRef PubMed 77. ↵ Chaudhuri , D. , Sasaki , K. , Karkar , A. , Sharif , S. , Lewis , K. , Mammen , M.J. , Alexander , P. , Ye , Z. , Lozano , L.E.C. , Munch , M.W. , et al. ( 2021 ). Corticosteroids in COVID-19 and non-COVID-19 ARDS: a systematic review and meta-analysis . Intensive Care Med 47 , 521 – 537 . OpenUrl CrossRef PubMed 78. ↵ Chang , X. , Li , S. , Fu , Y. , Dang , H. , and Liu , C . ( 2022 ). Safety and efficacy of corticosteroids in ARDS patients: a systematic review and meta-analysis of RCT data . Respir Res 23 , 301 . OpenUrl PubMed 79. ↵ Wen , Q. , Liu , J. , Kang , R. , Zhou , B. , and Tang , D . ( 2019 ). The release and activity of HMGB1 in ferroptosis . Biochem Biophys Res Commun 510 , 278 – 283 . OpenUrl CrossRef PubMed 80. ↵ Wilson , J.G. , and Calfee , C.S . ( 2020 ). ARDS Subphenotypes: Understanding a Heterogeneous Syndrome . Crit Care 24 , 102 . OpenUrl CrossRef PubMed 81. Sinha , P. , and Calfee , C.S . ( 2019 ). Phenotypes in acute respiratory distress syndrome: moving towards precision medicine . Curr Opin Crit Care 25 , 12 – 20 . OpenUrl CrossRef PubMed 82. ↵ Famous , K.R. , Delucchi , K. , Ware , L.B. , Kangelaris , K.N. , Liu , K.D. , Thompson , B.T. , Calfee , C.S. , and ARDS Network ( 2017 ). Acute Respiratory Distress Syndrome Subphenotypes Respond Differently to Randomized Fluid Management Strategy . Am J Respir Crit Care Med 195 , 331 – 338 . OpenUrl CrossRef PubMed 83. ↵ Herold , S. , Gabrielli , N.M. , and Vadász , I . ( 2013 ). Novel concepts of acute lung injury and alveolar-capillary barrier dysfunction . Am J Physiol Lung Cell Mol Physiol 305 , L665 – 81 . OpenUrl CrossRef PubMed Web of Science 84. ↵ Zhou , H. , Fan , E.K. , and Fan , J . ( 2021 ). Cell-Cell Interaction Mechanisms in Acute Lung Injury . Shock 55 , 167 – 176 . OpenUrl PubMed 85. ↵ Haghverdi , L. , Lun , A.T.L. , Morgan , M.D. , and Marioni , J.C . ( 2018 ). Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors . Nat Biotechnol 36 , 421 – 427 . OpenUrl CrossRef PubMed 86. ↵ Stuart , T. , Butler , A. , Hoffman , P. , Hafemeister , C. , Papalexi , E. , Mauck , W.M ., 3rd . , Hao , Y. , Stoeckius , M. , Smibert , P. , and Satija , R. ( 2019 ). Comprehensive Integration of Single-Cell Data . Cell 177 , 1888 – 1902.e21 . OpenUrl CrossRef PubMed 87. ↵ Muzellec , B. , Teleńczuk , M. , Cabeli , V. , and Andreux , M . ( 2023 ). PyDESeq2: a python package for bulk RNA-seq differential expression analysis . Bioinformatics 39 . doi: 10.1093/bioinformatics/btad547 . OpenUrl CrossRef PubMed 88. ↵ Shivram , H. , Hackney , J.A. , Rosenberger , C.M. , Teterina , A. , Qamra , A. , Onabajo , O. , McBride , J. , Cai , F. , Bao , M. , Tsai , L. , et al. ( 2023 ). Transcriptomic and proteomic assessment of tocilizumab response in a randomized controlled trial of patients hospitalized with COVID-19 . iScience 26 , 107597 . OpenUrl PubMed 89. ↵ Fleming , S.J. , Chaffin , M.D. , Arduini , A. , Akkad , A.-D. , Banks , E. , Marioni , J.C. , Philippakis , A.A. , Ellinor , P.T. , and Babadi , M . ( 2023 ). Unsupervised removal of systematic background noise from droplet-based single-cell experiments using CellBender . Nat Methods 20 , 1323 – 1335 . OpenUrl CrossRef PubMed 90. ↵ Germain , P.-L. , Lun , A. , Garcia Meixide , C. , Macnair , W. , and Robinson , M.D . ( 2021 ). Doublet identification in single-cell sequencing data using . F1000Res 10 , 979 . OpenUrl CrossRef 91. ↵ Website Gregoire Pau , J.R. ( 2017 ). HTSeqGenie (Bioconductor) . doi: 10.18129/B9.BIOC.HTSEQGENIE . OpenUrl CrossRef 92. ↵ Wu , T.D. , and Nacu , S . ( 2010 ). Fast and SNP-tolerant detection of complex variants and splicing in short reads . Bioinformatics 26 , 873 – 881 . OpenUrl CrossRef PubMed Web of Science 93. ↵ A multilevel clustering technique for community detection ( 2021 ). Neurocomputing 441 , 64 – 78 . OpenUrl 94. ↵ Hänzelmann , S. , Castelo , R. , and Guinney , J . ( 2013 ). GSVA: gene set variation analysis for microarray and RNA-seq data . BMC Bioinformatics 14 , 7 . OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted September 10, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about medRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Identification of robust mortality-associated neutrophil gene programs and cytokine responses using large-scale ARDS endotracheal aspirate scRNA-seq data Message Subject (Your Name) has forwarded a page to you from medRxiv Message Body (Your Name) thought you would like to see this page from the medRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Identification of robust mortality-associated neutrophil gene programs and cytokine responses using large-scale ARDS endotracheal aspirate scRNA-seq data Owen Whitley , Tim Delemarre , Miguel Reyes , Kimberly Kajihara , Hannah Little-Hooy , Cynthia Chen , Marco De Simone , Min Xu , Haridha Shivram , Jason Hackney , Diana Chang , Xia Gao , Eugene Kim , Nandhini Ramamoorthi , Rafal K. Woycicki , Shoshana Zha , Charles R. Langelier , Sharookh B. Kapadia , Carrie M. Rosenberger , Huwenbo Shi medRxiv 2025.08.21.25333911; doi: https://doi.org/10.1101/2025.08.21.25333911 Share This Article: Copy Citation Tools Identification of robust mortality-associated neutrophil gene programs and cytokine responses using large-scale ARDS endotracheal aspirate scRNA-seq data Owen Whitley , Tim Delemarre , Miguel Reyes , Kimberly Kajihara , Hannah Little-Hooy , Cynthia Chen , Marco De Simone , Min Xu , Haridha Shivram , Jason Hackney , Diana Chang , Xia Gao , Eugene Kim , Nandhini Ramamoorthi , Rafal K. Woycicki , Shoshana Zha , Charles R. Langelier , Sharookh B. Kapadia , Carrie M. Rosenberger , Huwenbo Shi medRxiv 2025.08.21.25333911; doi: https://doi.org/10.1101/2025.08.21.25333911 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Intensive Care and Critical Care Medicine Subject Areas All Articles Addiction Medicine (570) Allergy and Immunology (864) Anesthesia (302) Cardiovascular Medicine (4445) Dentistry and Oral Medicine (444) Dermatology (383) Emergency Medicine (609) Endocrinology (including Diabetes Mellitus and Metabolic Disease) (1515) Epidemiology (15236) Forensic Medicine (30) Gastroenterology (1127) Genetic and Genomic Medicine (6610) Geriatric Medicine (669) Health Economics (1000) Health Informatics (4549) Health Policy (1370) Health Systems and Quality Improvement (1613) Hematology (543) HIV/AIDS (1266) Infectious Diseases (except HIV/AIDS) (15926) Intensive Care and Critical Care Medicine (1104) Medical Education (623) Medical Ethics (147) Nephrology (668) Neurology (6613) Nursing (346) Nutrition (999) Obstetrics and Gynecology (1147) Occupational and Environmental Health (957) Oncology (3341) Ophthalmology (975) Orthopedics (369) Otolaryngology (420) Pain Medicine (436) Palliative Medicine (130) Pathology (665) Pediatrics (1694) Pharmacology and Therapeutics (693) Primary Care Research (714) Psychiatry and Clinical Psychology (5458) Public and Global Health (9244) Radiology and Imaging (2205) Rehabilitation Medicine and Physical Therapy (1370) Respiratory Medicine (1197) Rheumatology (596) Sexual and Reproductive Health (715) Sports Medicine (530) Surgery (713) Toxicology (99) Transplantation (289) Urology (265) (function(){function c(){var b=a.contentDocument||a.contentWindow.document;if(b){var d=b.createElement('script');d.innerHTML="window.__CF$cv$params={r:'a024a14f89d29839',t:'MTc3OTg4MDgyNA=='};var a=document.createElement('script');a.src='/cdn-cgi/challenge-platform/scripts/jsd/main.js';document.getElementsByTagName('head')[0].appendChild(a);";b.getElementsByTagName('head')[0].appendChild(d)}}if(document.body){var a=document.createElement('iframe');a.height=1;a.width=1;a.style.position='absolute';a.style.top=0;a.style.left=0;a.style.border='none';a.style.visibility='hidden';document.body.appendChild(a);if('loading'!==document.readyState)c();else if(window.addEventListener)document.addEventListener('DOMContentLoaded',c);else{var e=document.onreadystatechange||function(){};document.onreadystatechange=function(b){e(b);'loading'!==document.readyState&&(document.onreadystatechange=e,c())}}}})();