key: cord-0825119-03ybjkdv
authors: Sinha, Sarthak; Rosin, Nicole L.; Arora, Rohit; Labit, Elodie; Jaffer, Arzina; Cao, Leslie; Farias, Raquel; Nguyen, Angela P.; de Almeida, Luiz G. N.; Dufour, Antoine; Bromley, Amy; McDonald, Braedon; Gillrie, Mark; Fritzler, Marvin J.; Yipp, Bryan; Biernaskie, Jeff
title: An Immune Cell Atlas Reveals Dynamic COVID-19 Specific Neutrophil Programming Amenable to Dexamethasone Therapy
date: 2021-07-05
journal: bioRxiv
DOI: 10.1101/2021.04.18.440366
sha: 33285e0ffd8547a3b3c9da264f33b2aef96bffa9
doc_id: 825119
cord_uid: 03ybjkdv

SARS-CoV-2 is a novel coronavirus that causes acute respiratory distress syndrome (ARDS), death and long-term sequelae. Innate immune cells are critical for host defense but are also the primary drivers of ARDS. The relationships between innate cellular responses in ARDS resulting from COVID-19 compared to other causes of ARDS, such as bacterial sepsis is unclear. Moreover, the beneficial effects of dexamethasone therapy during severe COVID-19 remain speculative, but understanding the mechanistic effects could improve evidence-based therapeutic interventions. To interrogate these relationships, we developed an scRNA-Seq and plasma proteomics atlas (biernaskielab.ca/COVID_neutrophil). We discovered that compared to bacterial ARDS, COVID-19 was associated with distinct neutrophil polarization characterized by either interferon (IFN) or prostaglandin (PG) active states. Neutrophils from bacterial ARDS had higher expression of antibacterial molecules such as PLAC8 and CD83. Dexamethasone therapy in COVID patients rapidly altered the IFNactive state, downregulated interferon responsive genes, and activated IL1R2+ve neutrophils. Dexamethasone also induced the emergence of immature neutrophils expressing immunosuppressive molecules ARG1 and ANXA1, which were not present in healthy controls. Moreover, dexamethasone remodeled global cellular interactions by changing neutrophils from information receivers into information providers. Importantly, male patients had higher proportions of IFNactive neutrophils, a greater degree of steroid-induced immature neutrophil expansion, and increased mortality benefit compared to females in the dexamethasone era. Indeed, the highest proportion of IFNactive neutrophils was associated with mortality. These results define neutrophil states unique to COVID-19 when contextualized to other life-threatening infections, thereby enhancing the relevance of our findings at the bedside. Furthermore, the molecular benefits of dexamethasone therapy are also defined, and the identified pathways and plasma proteins can now be targeted to develop improved therapeutics.

SARS-CoV-2 is a novel coronavirus that causes acute respiratory distress syndrome (ARDS), death and 23 long-term sequelae. Innate immune cells are critical for host defense but are also the primary drivers of 24 ARDS. The relationships between innate cellular responses in ARDS resulting from COVID-19 25 compared to other causes of ARDS, such as bacterial sepsis is unclear. Moreover, the beneficial effects of 26 dexamethasone therapy during severe COVID-19 remain speculative, but understanding the mechanistic 27 effects could improve evidence-based therapeutic interventions. To interrogate these relationships, we 28 developed an scRNA-Seq and plasma proteomics atlas (biernaskielab.ca/COVID_neutrophil). We 29 discovered that compared to bacterial ARDS, COVID-19 was associated with distinct neutrophil 30 polarization characterized by either interferon (IFN) or prostaglandin (PG) active states. Neutrophils from 31 bacterial ARDS had higher expression of antibacterial molecules such as PLAC8 and CD83. 32

Dexamethasone therapy in COVID patients rapidly altered the IFN active state, downregulated interferon 33 responsive genes, and activated IL1R2 +ve neutrophils. Dexamethasone also induced the emergence of 34 immature neutrophils expressing immunosuppressive molecules ARG1 and ANXA1, which were not 35 present in healthy controls. Moreover, dexamethasone remodeled global cellular interactions by changing 36 neutrophils from information receivers into information providers. Importantly, male patients had higher 37

proportions of IFN active neutrophils, a greater degree of steroid-induced immature neutrophil expansion, 38 and increased mortality benefit compared to females in the dexamethasone era. Indeed, the highest 39 proportion of IFN active neutrophils was associated with mortality. These results define neutrophil states 40 unique to COVID-19 when contextualized to other life-threatening infections, thereby enhancing the 41 relevance of our findings at the bedside. Furthermore, the molecular benefits of dexamethasone therapy 42 are also defined, and the identified pathways and plasma proteins can now be targeted to develop 43 improved therapeutics. The online companion atlas (biernaskielab.ca/COVID_neutrophil) contains accessible scRNAseq data 78 performed on freshly obtained whole blood at timepoint 1 (t1, <72h after ICU admission) and at timepoint 79 2 (t2, 7 days after t1) (Fig 1a) . Cellular identity was mapped to 30 immune cell types/states using a 80 UMAP projection from 21 patients and 86,935 cells (Fig 1b, Extended Figure 3a ). Global magnitude of 81 gene expression was directly compared between COVID-19 and bacterial ARDS patients (Extended Data 82 Table 3 ), which revealed a more globally altered distribution of differential expression at t1 than at t2. 83

Altered regulation of genes was most pronounced in neutrophils at t1, with lower neutrophil gene 84 expression in COVID-19 compared to bacterial ARDS (Fig 1c; Extended Fig 3b-c) . At t2, the global 85 alterations in gene expression when comparing COVID-19 to bacterial ARDS were most pronounced in 86 plasmablasts (Fig 1d; Extended Fig 3d-e) . We further compared and quantified the proportions of known 87 peripheral blood cellular constituents, which highlighted significant differences in CD4 T cells, CD8 T 88 cells and NK cells (Extended Fig 3f) . These data highlight that significant global differences in immune 89 cell gene expression exist between COVID-19 ARDS and bacterial ARDS. 90

Neutrophils are a primary participant in the development of ARDS 18 ; yet despite similar severity of 92 ARDS between our bacterial and COVID-19 cohorts, the numbers of circulating neutrophils from clinical 93 cell counts were significantly different (Extended Fig 2d) . We hypothesized that neutrophil qualitative 94 states may be important determinants of disease. Neutrophils were subjected to velocity analysis 19,20 to 95 reconstruct maturation dynamics. Louvain clusters (Fig 1e) , clinical cohort, individual patient, and 96 velocity length were overlayed on velocity vector fields ( Extended Fig 4a-d) . The proportions of distinct 97 neutrophil states were compared at t1 and this revealed a divergent expansion of IFN active neutrophils 98 (clusters 2, 4 and 5) marked by IFITM1 expression in COVID-19, which became similar to bacterial 99 ARDS at t7 (Fig 1f-h) . Expression of IFITM1 in neutrophils from COVID-19 patients at t1 was 100 confirmed by immunofluorescent staining for IFITM1, colocalized with S100A8/9 and typical neutrophil 101 nuclear morphology. Relative to healthy donors, the IFN active population in both COVID-19 and bacterial 102 ARDS patients were elevated (Extended Fig 4h- states. h. Immunocytochemistry for S100A8/A9 (red) and IFITM1 (green) expression on leukocyte-rich 119 preparation from COVID-19 donor at tl. i-k. Transcriptional kinetics driving expansion of IFN active (i), 120

Bacterial ARDS-enriched (j), and PG active (k) neutrophils. Latent time distribution of trajectory-associated 121 louvain clusters (left), phase portraits with equilibrium slopes of spliced-unspliced ratios (center), and 122

RNA velocity and gene expression (right) of selected genes driving divergent maturation trajectories. 123

Phase portraits are coloured by clinical cohort. 124

Classically, peripheral neutrophils are considered terminally differentiated and non-dividing, however the 125 increase in velocity length suggested the ability to alter phenotypic states once in circulation along 126 specific paths or 'lineages'. COVID-19 neutrophils followed unique maturation paths compared to 127 bacterial ARDS, culminating in three distinct terminal states: Interferon active (IFN active ), prostaglandin 128 active (PG active ) or bacterial ARDS enriched (Fig 1e-g; Extended Fig 4e) . Interestingly, the apex of this 129 trajectory was marked by high velocity lengths, characteristic of cells undergoing differentiation 130 (Extended Fig 4c, with a uniform proinflammatory capacity. Global neutrophil expression aligned with neutrophil state 147 specific markers, such as interferon response genes (IFITM1, RSAD2, IFI6, and ISG10), being more 148 highly expressed in COVID-19 neutrophils (Fig 2a; Extended Fig 4f) . The inverse was the case for anti-149 bacterial proteins like PLAC8 (Fig 2a; Online Atlas). However, the discovery of differential neutrophil 150 states prompted further exploration of the factors driving neutrophil state polarization. Gene regulatory 151 network reconstruction using SCENIC analysis 26 revealed differentially activated transcription factors 152 STAT1, IRF2 and PRDM1 in COVID-19 (Fig 2b) , while bacterial ARDS neutrophils had increased 153 prototypical granulocyte transcription factors such as CEBPA, CEBPB, STAT5B and less defined factors 154 such as NFE2 (Fig 2b, Online Atlas). PRDM1 activation was most pronounced in the IFN active neutrophil 155 population and was likely responsible for driving expression of interferon response elements (IFIT1, 156 ISG15, IFI6) and antiviral signaling, such as RSAD2 and STAT1 (Fig 2c; Online Atlas). A hallmark of 157 PG active neutrophil polarization was the activation of an E2F4 pathway (Fig 2d) , while neutrophil 158 programming during bacterial ARDS included STAT5B (Fig 2e) . To summarize, in response to COVID-159 19, neutrophils were polarized by unique transcriptional regulation towards one of two main populations, 160 either an IFN active population or a PG active population (Fig 2f) . 161 was associated with an increase in cytotoxic CD4 T cells, naïve B cells, plasmablasts, and decreased 187 proliferating NK cells, and CD4 effector memory cells (Extended Fig 5g) . By t2, dexamethasone was 188 associated with suppressed neutrophil proportions in circulation compared to untreated controls (13% vs 189 41%, Extended Fig 5g) . Plasma proteomics from the same cohort revealed that dexamethasone 190 suppressed 10 host proteins (S100A8, S100A9, SERPINA1, SERPINA3, ORM1, LBP, VWF, PIGR, 191 AZGP1, CRP) that others have previously identified as biomarkers distinguishing severe COVID-19 192 cases from mild to moderate counterparts (full host proteome quarriable via Online Atlas; Extended Table  193 2) 28-31 . Suppression of calprotectin (S100A8/S100A9) and neutrophil serine proteases (SERPINA1 and 194 SERPINA3), paired with depletion of neutrophil proportions, implicates the modulation of neutrophil-195 related inflammatory processes as a method of action for dexamethasone treatment. Interestingly, we also identified IL7R +ve neutrophils (comprising roughly 8% of total neutrophils) whose 221 trajectories remained completely separate ( administration at t1 halted dynamic state changes in IFN active and IL7R +ve neutrophils, followed by 227 preferential depletion of IFN active subsets (Fig 3 g) . Indeed, dexamethasone was associated with a 228 reduction in IFN active neutrophils to a proportion more similar to that detected in healthy controls (9% 229 post-Dex at t2 versus 10% in healthy controls) (Fig. 4a, Extended Fig 4h-k) . Although collection of 230 airway samples (i.e. bronchoalveolar lavage fluid; BALF) was not feasible at our institution, we leveraged 231 two recent BALF scRNA-Seq datasets 11,32 to assess whether IFN active neutrophils dominate the 232 bronchoalveolar landscape during severe COVID-19. Projection of CSF3R + S100A8 + S100A9 + BALF 233 neutrophils onto our reference revealed: a. 1.5 FC expansion of IFN active neutrophils in severe COVID-19 234 relative to moderate disease (77% vs 52%, Extended Fig 7a-b) , b. preferential activation of IFN-235 stimulated genes such as IFITM1, IFITM2, IFI6, IRF7, and ISG20 in severe COVID-19 neutrophils 236 (Extended Fig 7c) , and c. 4.7 FC greater IFN active neutrophils in COVID-19 relative to bacterial 237 pneumonia patients (14% vs 3%, Extended Fig 7d-f ). Albeit anecdotal, in our whole blood cohort, the 238 IFN active neutrophil state was dominant in patient S7 32 , an 80-year-old male with remarkably high viral 239 titers who succumbed to COVID-19 complications within 3-4 days of sampling (Extended Fig 7f) . 240

Consensus DEG analysis highlighted that upregulation of IL1R2, a decoy receptor that sequesters IL-1, 241

and downregulation of IFITM1 were the most prominent discriminating features of treatment with 242 steroids (Fig. 3h) . Additionally, dexamethasone attenuated neutrophil expression of IFN pathways more 243 broadly, including the reduction of IFITM1-3, IFIT1, ISG15 and RSAD2 (Fig 3h) . Examination of 244 were driven by differential splicing kinetics. 247

Corticosteroid therapy shifted neutrophil state compositions. While IFN active neutrophils were significantly 249 depleted by seven days of therapy, there was >2-fold expansion in immature neutrophils relative to non-250 treated COVID-19 controls (Fig 4a; Extended Fig 6 h, i) , which were absent in the healthy controls. 251

Albeit anecdotal, the dominance of IFN active neutrophils at t1 in the patient who succumbed to COVID-19 252 in the non-dexamethasone cohort further supports depletion of IFN active neutrophils as a mechanism by 253 which dexamethasone is protective ( Extended Fig 8 g-j) . Assessment of gene regulatory networks 254 demonstrated that IRF7 and MEF2A exhibited opposing activation patterns, with IRF7 being the most 255 suppressed and MEF2A the most enhanced transcription factors identified with dexamethasone, which 256 correlates with the emergence of PG active and IL1R2 +ve states and attenuation of the IFN active neutrophil 257 states (Fig 4b, Extended Fig 6k-m) . To assess the generalizability of the dexamethasone regulated DEGs 258 identified in our cohort, we asked whether they accurately predicted mortality due to COVID-19 in a 259 larger validation cohort. By leveraging a whole blood bulk RNA-Seq dataset from 103 COVID-19 260 patients 33 34 , we scored each sample by the aggregated expression of dexamethasone suppressed DEGs at 261 t1 and t2 (Extended Data Table 3 ). Interestingly, suppressed DEGs at t2 (but not t1) proved to be a far 262

superior predictor of 28-day mortality (AUC: 0.78, CI: 0.67 -0.89) compared to clinical severity scales 263 such as sequential organ failure assessment (SOFA) (AUC: 0.67, CI: 0.51-0.82) across all classification 264 thresholds (Fig 4c) . 265

Unexpectedly, steroid administration was associated with an increase in circulating immature neutrophils, 266 which highly expressed TOP2A, and activated ATF4 and JDP2, transcription factors seen in 267 undifferentiated cells or those undergoing nuclear reprogramming (Extended Fig 6h) . Interestingly, these 268 immature neutrophils expressed high levels of ARG1, ANXA1 (Fig 4d) , and CD24 (both mRNA and Fig 9c, d) . Interestingly, while neutrophils were depleted in both sexes post-297 dexamethasone, this was particularly pronounced in males (1.9 FC higher in males at t1 and 3.4 FC 298 higher in males at t2, Extended Fig 9e) . Of the two salient neutrophil state alterations, an immature 299 (ARG1 +ve immunosuppressive) state was preferentially expanded with dexamethasone in males (Extended 300 Fig. 9e) , whereas ISGs were preferentially suppressed (Extended Fig. 9f ) and IFN active states were 301 depleted in females (Extended Fig. 9g-h) at both t1 and t2 (Fig 4h, i) . Sexually dimorphic effects of 302 dexamethasone on neutrophil maturation kinetics may in part explain these state alterations. Dynamo-303 reconstructed vector dynamics revealed that dexamethasone slowed IFN active transitions (Extended Fig. 9i ) 304

whilst accelerating immature (ARG1 +ve immunosuppressive) neutrophil differentiation in females 305 (Extended Fig. 9j ) ultimately leading to a diminished immature neutrophil progenitor pool. 306 TFs activated or suppressed post-dexamethasone in at least 3 of 6 patients at t1 and predicted activity of 311 MEF2A and IRF7, two of the most differentially regulated TFs post-dexamethasone. c. Receiver 312 operating characteristic (ROC) curves assessing the discriminatory capacity of dexamethasone suppressed 313

DEGs at t1, t2, and sequential organ failure assessment (SOFA) scores for predicting 28-day mortality in 314 a validation cohort of 103 bulk whole blood RNA-Seq samples where 17 cases were fatal. d. Immature  315 and IL1R2 +ve neutrophil subsets express high levels of immunosuppressive neutrophil marker ARG1 and 316

ANXA1. e. Neutrophil-driven signaling pathways induced post-dexamethasone, identified using CellChat 317

(MHC-I signalling filtered out). f, g. Topology of annexin signalling without (e) and with dexamethasone 318

(f) treatment (edges filtered to those where neutrophils function as senders or recipients of annexin 319 signals). h. Neutrophil state composition separated by sex and dexamethasone status at t1 and t2. i. 320

Schematic summarizing the effects of dexamethasone on neutrophil fates and function in COVID-19 321

following dexamethasone treatment. The research teams did not participate in clinical decisions. Study inclusion required a minimal age of 18, 362 the ability to provide consent, or for most participants, the ability of a surrogate decision maker to provide 363 regained capacity consent. All participants required an arterial catheter for blood draws, but the insertion 

The serum of COVID-19 patients (COVID-19 = 9, dexamethasone-treated = 4) and bacterial ARDS 412 controls (N = 6) were collected and subjected to quantitative proteomics. The total protein concentrations 413

were determined by Pierce™ BCA Protein Assay Kit (23225, ThermoFisher). A trichloroacetic acid 414 (TCA)/acetone protocol was used to pellet 100µg of proteins per sample. Samples were subjected to a 415 quantitative proteomics workflow as per supplier (Thermo Fisher) recommendations. Samples were 416 reduced in 200mM tris(2-carboxyethyl)phosphine (TCEP), for 1h at 55°C, reduced cysteines were 417 alkylated by incubation with iodoacetamide solution (50mM) for 20min at room temperature. Samples 418 were precipitated by acetone/methanol, and 600µL ice-cold acetone was added followed by incubation at 419 -20°C overnight. A protein pellet was obtained by centrifugation (8,000 , 10min, 4°C) followed by 420 acetone drying (2min). Precipitated pellet was resuspended in100 µL of 50mM triethylammonium 421 bicarbonate (TEAB) buffer followed by tryptase digestion (5µg trypsin per 100µg of protein) overnight at 422 37°C. TMT-6plex™ Isobaric Labeling Reagents (90061, Thermo Fisher) were resuspended in anhydrous 423 acetonitrile and added to each sample (41µL TMT-6plex™ per 100µL sample) and incubated at room 424 temperature for 1h. The TMT labeling reaction was quenched by 2.5% hydroxylamine for 15min at room 425 temperature. TMT labeled samples were combined and acidified in 100% trifluoroacetic acid to pH < 3.0 426 and subjected to C18 chromatography (Sep-Pak) according to manufacturer recommendations. Samples 427 were stored at -80°C before lyophilization, followed by resuspension in 1% formic acid before liquid 428 chromatography and tandem mass spectrometry analysis. of 120,000 FWHM to detect the precursor ion having a m/z between 375 and 1,575 and a +2 to +4 charge. 442

The Orbitrap AGC (Auto Gain Control) and the maximum injection time were set at 4 x 10 5 and 50ms, 443 respectively. The Orbitrap was operated using the top speed mode with a 3 second cycle time for 444 precursor selection. The most intense precursor ions presenting a peptidic isotopic profile and having an 445 intensity threshold of at least 2 x 10 4 were isolated using the quadrupole (Isolation window (m/z) of 0.7) 446 and fragmented using HCD (38% collision energy) in the ion routing multipole. The fragment ions (MS2) 447

were analyzed in the Orbitrap at a resolution of 15,000. The AGC and the maximum injection time were 448 set at 1 x 10 5 and 105ms, respectively. The first mass for the MS2 was set at 100 to acquire the TMT 449 reporter ions. Dynamic exclusion was enabled for 45 seconds to avoid of the acquisition of same 450 precursor ion having a similar m/z (plus or minus 10ppm). 451

Spectral data acquired from the mass spectrometer were matched to peptide sequences using MaxQuant 453 software (v.1.6.14) 49 . Due to a lack of direct compatibility with Maxquant, spectra generated using the 454 rehydrated. Slides were permeabilized and blocked with 10% normal donkey serum in PBS (with 0.5% 497 triton X-100), primary antibodies (S100A8/9 Abcam ab22506; IFITM1 Abcam ab233545) were 498 incubated at 4 o C overnight, followed by incubation with donkey anti-rabbit-Alexa488 (Invitrogen 499 A32790) or anti-mouse-Alexa555 (Invitrogen A31570) for 1h at room temperature (RT). Cytospun slides 500

were sequentially stained with CD24 (Abcam ab202073) on the same slides for 1h at RT, followed by 501 donkey anti-rabbit-Alexa647 (Invitrogen A31573). Imaging was done using a VS-120 slide scanner 502 (Olympus) and high resolution image imaging was done using an SP8 spectral confocal microscope 503 (Leica). Image processing was completed in Fiji 52 . 504

Single-cell RNA-Seq library construction, alignment, and quality control. A total of 15,000 single 505 cells (containing an equal proportion of leukocytes and lymphocytes) were loaded for partitioning using 506 10X Genomics NextGEM Gel Bead emulsions (Version 3.1). All samples were processed as per 507 manufacturer's protocol (with both PCR amplification steps run 12X). Quality control of resulting 508 libararies and quantification was performed using TapeStation D1000 ScreenTape assay (Agilent). 509

Sequencing was performed using Illumina NovaSeq S2 and SP 100 cycle dual lane flow cells over 510 multiple rounds to ensure each sample received approximately 32,000 reads per cell. Sequencing reads 511 were aligned using CellRanger 3.1.0 pipeline 53 to the standard pre-built GRCh38 reference genome. 512

Samples that passed alignment QC were aggregated into single datasets using CellRanger aggr with 513 between-sample normalization to ensure each sample received an equal number of mapped reads per cell. 514 Single-cell RNA-Seq computational analyses and workflows. Filtered feature-barcode HDF5 matrices 520 from aggregated datasets were imported into the R package Seurat v.3.9 for normalization, scaling, 521 integration, multi-modal reference mapping, louvain clustering, dimensionality reduction, differential 522 expression analysis, and visualization 54 . Briefly, cells with abnormal transcriptional complexity (fewer 523 than 500 UMIs, greater than 25,000 UMIs, or greater than 25% of mitochondrial reads) were considered 524 artifacts and were removed from subsequent analysis. Since granulocytes have relatively low RNA 525 content (due to high levels of RNases), QC thresholds were informed by 8 as they recently defined several 526 rodent and human neutrophil subsets from scRNA-Seq samples. Cell identity was classified by mapping 527 single cell profiles to the recently published PBMC single-cell joint RNA/CITE-Seq multi-omic reference 528 55 . 529

Annotation of neutrophil states. Since no published reference automates granulocyte annotations, 530 neutrophil clusters were manually annotated by querying known markers (i.e. CSF3R, S100A8, S100A9, 531 MMP8, MMP9, ELANE, MPO) 56 and were corroborated using the R package SingleR 57 . Neutrophil 532 states were defined by grouping unsupervised (louvain at default resolution) subclusters based on two 533 overlapping criteria: scVelo-inferred neutrophil maturity, and 2. by corroborating gene expression and 534 SCENIC-inferred GRN signatures with previous human and rodent neutrophil scRNA-Seq studies. 535

Immature neutrophils were defined as CD24 + ARG1 + ELANE + MPO + ATF4 GRN-active JDP2 GRN-active 536 neutrophils 7,8,47,58 that were reproducibly assigned as 'root cells' in scVelo-based latent time pseudo-537 ordering. IFN active neutrophils were defined by preferential mRNA splicing (positive velocity) and 538 expression of IFN-stimulated genes such as IFITM1/2, IFIT1/2/3, ISG15/20, and IFI6/27/44/44L 6,44,59 . 539 PG active neutrophils were distinguished by preferential splicing of PTGS2/COX2 (as well as expression for 540 prostaglandin transport LST1) 44 and included a subset that expressed high levels of IL1β decoy receptor 541 IL1R2 33 . Lastly, IL7R + neutrophils (a small but distinct subset that maybe of thymic origin 60 expressed 542 high levels of ribosomal subunit genes (e.g. RPL5/7A/8/13/18/19/23/24/27/P0) that are highly 543 reminiscent of 'ribosomal hi -specific cluster 7' identified previously 47 .. 544

Statistical approach for comparing cell proportions. To test whether cell composition was changed 545 due to infection type (COVID-19 versus Bacterial ARDS) or treatment group (dexamethasone versus 546 non-dexamethasone), a generalized linear mixed-effects model was employed where infection type and 547 treatment group were considered fixed and individual patients were considered random effect. Fitting was 548 done with Laplace approximation using the 'glmer' function in the 'lme4' R package 61 and p-values were 549 calculated using the R package 'car'. Boxplots comparing cell type composition were generated using the 550 ggplot2 package. Since a subset of patients sampled at t1 were discharged from ICU prior to t2 collection 551 (non-random or non-ignorable missing data), we limit statistical comparisons to between group 552 comparisons within one time point (e.g., COVID-19 72h vs Bacterial ARDS 72hr, dexamethasone-treated 553 72h vs non-dexamethasone-treated 72h) and do not estimate temporal differences across t1 and t2. reproducible changes across patients (>3 for 72-hour comparisons, > 2 for 7-day comparisons). Gene Set 566

Enrichment analyses of consensus DEGs were performed using gProfiler's g:GOSt (p-value cutoff 567 <0.05). A cell state-specific 'perturbation score' was calculated to reflect the magnitude of response 568 elicited by factoring in number and cumulative FC of consensus DEGs. Perturbation scores were 569 visualized using Nebulosa-generated density plots 64 . 570

Constructing cellular trajectories using RNA velocity. Analysis of neutrophil trajectories was 571 performed by realigning CellRanger count-generated BAMs with RNA velocity command-line tool 20 572 using the run10x command and human (GRCh38) annotations. The output loom files containing spliced 573 and unspliced counts were combined to compare neutrophils in COVID-19 with Bacterial ARDS controls 574 and dexamethasone-treated with non-treated COVID-19 patients. For both analyses, combined looms 575 were imported into Seurat v.3.9 using the ReadVelocity function in SeuratWrappers v.0.2.0, normalized 576 using SCTransform v.0.3.2 65 , reduced and projected onto a UMAP, and exported as a .h5 file using the 577 SaveH5Seurat function. Counts stored in H5 files were imported, filtered, and normalized as 578 recommended in the scVelo v.0.2.1 workflow 19 . RNA velocities were estimated using stochastic and 579 dynamical models. Since both models yielded comparable results, stochastic model was used as default 580 for all subsequent analyses. Calculations stored in AnnData's metadata were exported as CSVs and kernel 581 density lines depicting Velocity-inferred latent time distribution were plotted with ggplot2 v.3.1.1. 582

Gene Regulatory Network reconstruction. Single-cell regulatory network inference and clustering 583 (SCENIC) 26 was employed to infer regulatory interactions between transcription factors (TFs) and their 584 targetome by calculating and pruning co-expression modules. Briefly, neutrophils were subsetted from 585 scVelo-realigned Seurat object and processed using default and recommended parameters specified in 586 SCENIC's vignette (https://github.com/aertslab/SCENIC) using the hg19 RcisTarget reference. Regulon 587 activity scores (in '3.4_regulonAUC.Rds', an output of the SCENIC workflow) were added to scVelo 588 object (using CreateAssayObject function) to jointly project trajectory and TF activity onto the same 589 UMAP embeddings. Consensus stacked bars showing cumulative logFC of AUCell scores for each TF 590

(colored by individual sample contributions) were generated by modifying the constructConsensus 591 function 7 for SCENIC assay. Targetome of TFs predicted as drivers of neutrophil states (stored in 592 '2.6_regulons_asGeneSet.Rds') was profiled using g:Profiler's functional enrichment analysis and genes 593 intersecting with the Interferon pathway were plotted using iRegulon (Cytoscape plugin) 66 . 594

Comparing scRNA-Seq findings with published datasets. To test whether dexamethasone-suppressed 595 neutrophil genes at t1 and t2 (Extended Data Table 4 ) predicted COVID-19 mortality, we repurposed 596 methods described in 33 and employed whole blood bulk RNA-Seq datasets generated by 34 as a validation 597 cohort of 103 samples (where 17 were fatal). Briefly, each of the 103 samples were scored by the 598 aggregated expression of dexamethasone-suppressed neutrophil consensus genes at t1 and t2 using 599

Seurat's AddModuleScore(). Dexamethasone-suppressed module scores were used as the predictor 600 variable and 28-day mortality was used as the response variable to construct an ROC curve using pROC's 601 roc() function. To infer bronchoalveolar neutrophil composition in severe and moderate COVID-19 11 and 602 across bacterial pneumonia and COVID-19 32 , neutrophils (CSF3R + , S100A8 + , S100A9 + ) captured in 603 BALF scRNA-Seq datasets were projected onto our peripheral blood reference using mutual nearest 604 neighbor anchoring (FindTransferAnchors) and identity transferring (TransferData and AddMetaData) 605 strategy implemented in Seurat v4 54 . 606 COVID Neutrophil Atlas. To enable intuitive exploration of single-cell datasets, a web portal 607 (http://biernaskielab.ca/covid_neutrophil or http://biernaskielab.com/covid_neutrophil) was built using 608 

Seq datasets employed as an independent validation cohort were downloaded from GSE157103. BALF 615 scRNA-Seq datasets from severe and moderate COVID-19 were downloaded from GSE145926

Processed BALF scRNA-Seq objects from patients with bacterial pneumonia and COVID-19

GSE167118) were downloaded from authors

Mass spectrometry datasets will be available via PRIDE Archive

pride/), it has been submitted 621 (submission #: 1-20210702-114055) and is pending accessioning

All analyses were performed using publicly available software as described in the 623 methods section

Acknowledgements: This work was funded by a FastGrant from the Thistledown Foundation

BY) and Calgary Firefighters Burn Treatment Society (JB). S Sinha received CIHR Vanier

Killam doctoral scholarships. E.L received an Alberta Children's Hospital Research 628

Y is a tier II Canada Research Chair in Pulmonary Immunology

We acknowledge the assistance of the nurse practitioners

Kirsten Deemer and Robert Ralph as well as the healthcare teams from the Calgary Adult 631

ICU's at South Health Campus, Rockyview General Hospital

We thank Dr. Kirsten Fiest and the ICU study coordinators Cassidy Codan

We acknowledge Dan Jones, Cathy Curr and the eCritical team

Canada) for their help in data acquisition and extraction via eCritical 635 databases. Mortality predictions using dexamethasone-suppressed gene signatures were completed by 636 repurposing computational workflows kindly shared by

Acute respiratory distress syndrome

Clinical Characteristics of 138 Hospitalized Patients With

A Novel Coronavirus from Patients with Pneumonia in China

Neutrophil extracellular traps contribute to immunothrombosis in COVID-656 19 acute respiratory distress syndrome

SARS-CoV-2-triggered neutrophil extracellular traps mediate COVID-19 659 pathology

Global absence and targeting of protective immune states in severe COVID-661 19

A single-cell atlas of the peripheral immune response in patients with severe 663 COVID-19

Single-cell transcriptome profiling reveals neutrophil heterogeneity in homeostasis 665 and infection

Immunophenotyping of COVID-19 and influenza highlights the role of type I 667 interferons in development of severe COVID-19

The differential immune responses to COVID-19 in peripheral and lung revealed by 670 single-cell RNA sequencing

Single-cell landscape of bronchoalveolar immune cells in patients with COVID-19

Discriminating mild from critical COVID-19 by innate and adaptive immune 674 single-cell profiling of bronchoalveolar lavages

SARS-CoV-2 infection of circulating immune cells is not responsible for virus 677 dissemination in severe COVID-19 patients. bioRxiv

Acute respiratory distress syndrome: the Berlin Definition

COVID-19: consider cytokine storm syndromes and immunosuppression

Early prediction of disease progression in COVID-19 pneumonia patients with 684 chest CT and clinical characteristics

Risk Factors Associated With Acute Respiratory Distress Syndrome and Death in 687 Patients With Coronavirus Disease

Targeting potential drivers of COVID-19: Neutrophil extracellular traps

Generalizing RNA velocity to transient 692 cell states through dynamical modeling

RNA velocity of single cells

Non-steroidal anti-inflammatory drugs dampen the cytokine and antibody 697 response to SARS-CoV-2 infection

Neutrophil plasticity: acquisition of phenotype and functionality of 700 antigen-presenting cell

Impaired Host Defense in Mice Lacking ONZIN. The 702

Inborn errors of type I IFN immunity in patients with life-threatening

SCENIC: single-cell regulatory network inference and clustering

Dexamethasone in Hospitalized Patients with Covid-19 -Preliminary Report

Circulating Calprotectin as a 711 Biomarker of COVID-19 Severity

A time-resolved proteomic and prognostic map of COVID-19

In-depth blood proteome profiling analysis revealed distinct functional 716 characteristics of plasma proteins between severe and non-severe COVID-19 patients

Plasma Proteomics Identify Biomarkers and Pathogenesis of COVID-19

Clonal expansion and activation of tissue-resident memory-like T<sub>H</sub>17 721 cells expressing GM-CSF in the lungs of patients with severe COVID-19

Multi-omic profiling reveals widespread dysregulation of innate immunity and 724 hematopoiesis in COVID-19

Large-Scale Multi-omic Analysis of COVID-19 Severity

Annexin A1 modulates natural and glucocorticoid-induced resolution of 729 inflammation by enhancing neutrophil apoptosis

Annexin A1 Is Involved in the Resolution of Inflammatory Responses during 732 <em>Leishmania braziliensis</em> Infection

Early Expansion of Circulating Granulocytic Myeloid-derived Suppressor Cells 735 Predicts Development of Nosocomial Infections in Patients with Sepsis

Arg1 expression defines immunosuppressive subsets of tumor-associated 738 macrophages

Arginase 1 (Arg1) as an Up-Regulated Gene in COVID-19 Patients: A 740 Promising Marker in COVID-19 Immunopathy

ER Stress Regulates Immunosuppressive Function of Myeloid Derived 742

Suppressor Cells in Leprosy that Can Be Overcome in the Presence of IFN-γ. iScience 23

Molecular regulation of urea cycle function by the liver glucocorticoid receptor

Male sex identified by global COVID-19 meta-analysis as a risk factor for death 747 and ITU admission

Type I and Type III Interferons -Induction, Signaling, Evasion, and 749 Application to Combat COVID-19

Sex differences in neutrophil biology modulate response to type I interferons and 752 immunometabolism

Safety and efficacy of inhaled nebulised interferon beta-1a (SNG001) for 755 treatment of SARS-CoV-2 infection: a randomised, double-blind, placebo-controlled, phase 2 756 trial

Retrospective Multicenter Cohort Study Shows Early Interferon Therapy Is 758 Associated with Favorable Clinical Responses in COVID-19 Patients

Severe COVID-19 Is Marked by a Dysregulated Myeloid Cell 761 Compartment

Secondary EMR data for quality improvement and research: A 763 comparison of manual and electronic data collection from an integrated critical care electronic 764 medical record system

MaxQuant enables high peptide identification rates, individualized p.p.b.-766 range mass accuracies and proteome-wide protein quantification

Andromeda: a peptide search engine integrated into the MaxQuant environment

MSstatsTMT: Statistical Detection of Differentially Abundant Proteins in 771 Experiments with Isobaric Labeling and Multiple Mixtures

Fiji: an open-source platform for biological-image analysis

Massively parallel digital transcriptional profiling of single cells

Comprehensive Integration of Single-Cell Data

Integrated analysis of multimodal single-cell data. bioRxiv, 2020

Single-Cell Transcriptomics of Human and Mouse Lung Cancers Reveals 782

Reference-based analysis of lung single-cell sequencing reveals a transitional 785 profibrotic macrophage

Neutrophils in COVID-19

Single-Cell Transcriptome Analysis 790 Highlights a Role for Neutrophils and Inflammatory Macrophages in the Pathogenesis of Severe 791 COVID-19

Fate Mapping Reveals Separate Origins of T Cells and Myeloid Lineages in 793 the Thymus

Fitting Linear Mixed-Effects Models Using lme4

<em>Connectome</em>: computation and visualization of cell-cell 797 signaling topologies in single-cell systems data. bioRxiv

Inference and analysis of cell-cell communication using CellChat

Nebulosa recovers single-cell gene expression signals by 802 kernel density estimation

Normalization and variance stabilization of single-cell RNA-seq data 804 using regularized negative binomial regression

iRegulon: From a Gene List to a Gene Regulatory Network Using Large Motif 807 and Track Collections

The authors have no competing interests. 646