key: cord-1047253-u29717en
authors: Ong, Eugenia Z; Kalimuddin, Shirin; Chia, Wen Chong; Ooi, Sarah H; Koh, Clara WT; Tan, Hwee Cheng; Zhang, Summer L; Low, Jenny G; Ooi, Eng Eong; Chan, Kuan Rong
title: Temporal dynamics of the host molecular responses underlying severe COVID-19 progression and disease resolution
date: 2021-03-07
journal: EBioMedicine
DOI: 10.1016/j.ebiom.2021.103262
sha: ea7da3acf62baa14d0650bae5f0da0f59f953519
doc_id: 1047253
cord_uid: u29717en

BACKGROUND: The coronavirus disease-19 (COVID-19) pandemic has cost lives and economic hardships globally. Various studies have found a number of different factors, such as hyperinflammation and exhausted/suppressed T cell responses to the etiological SARS coronavirus-2 (SARS-CoV-2), being associated with severe COVID-19. However, sieving the causative from associative factors of respiratory dysfunction has remained rudimentary. METHODS: We postulated that the host responses causative of respiratory dysfunction would track most closely with disease progression and resolution and thus be differentiated from other factors that are statistically associated with but not causative of severe COVID-19. To track the temporal dynamics of the host responses involved, we examined the changes in gene expression in whole blood of 6 severe and 4 non-severe COVID-19 patients across 15 different timepoints spanning the nadir of respiratory function. FINDINGS: We found that neutrophil activation but not type I interferon signaling transcripts tracked most closely with disease progression and resolution. Moreover, transcripts encoding for protein phosphorylation, particularly the serine-threonine kinases, many of which have known T cell proliferation and activation functions, were increased after and may thus contribute to the upswing of respiratory function. Notably, these associative genes were targeted by dexamethasone, but not methylprednisolone, which is consistent with efficacy outcomes in clinical trials. INTERPRETATION: Our findings suggest neutrophil activation as a critical factor of respiratory dysfunction in COVID-19. Drugs that target this pathway could be potentially repurposed for the treatment of severe COVID-19. FUNDING: This study was sponsored in part by a generous gift from The Hour Glass. EEO and JGL are funded by the National Medical Research Council of Singapore, through the Clinician Scientist Awards awarded by the National Research Foundation of Singapore.

The novel coronavirus, SARS-CoV-2, has emerged to threaten global health and economic security. Following viral infection, individuals can be asymptomatic or present with mild symptoms such as fever, dry cough and lethargy. In some cases, patients develop severe pulmonary dysfunction, including acute respiratory distress syndrome (ARDS), which is associated with poor prognosis. Clinical management of patients that progress to respiratory failure involves intensive care management and reliance on ventilators, which can impose a significant burden on healthcare systems that are already overstretched [1] . Pathogenic mechanism guided approach to expand the extremely limited therapeutic options for COVID-19 could mitigate the present negative impact of this pandemic on global health and economies.

Rapid progress in research has identified several factors that are associated with severe COVID-19. Increased inflammatory monocytes and neutrophils, accompanied by a sharp decline in lymphocytes, including T cell count and activation, have been found to be associated with severe COVID-19 [2À5] . Severe COVID-19 patients have also shown heightened production of proinflammatory cytokines (TNF, IL-6 and IL-1b) and interferon (IFN) [5] . However, these responses are highly dynamic and changes daily following illness onset [6] . Their statistical association with severe COVID-19 may be heavily influenced by sampling points in the course of acute illness. Consequently, efficacy outcome of trials that tested inhibitors of hyperinflammation have been, at best, mixed. For instance, despite the consistent findings of elevated IL-6 expression in COVID-19 patients, IL-6 inhibitors (tocilizumab and sarilumab) failed to show efficacy in preventing severe COVID-19 [7, 8] . IFN therapy reduced mortality if administered early, but late administration potentially worsens disease outcome [9] . IL-1 receptor antagonist (Anakinra) and TNF antagonists are promising candidates [10, 11] but their effectiveness have yet to be evaluated in randomised, placebo-controlled trials. Among the licensed corticosteroids, only dexamethasone has shown efficacy in reducing COVID-19 mortality [12] . Moreover, T cell related therapies, such as the use of checkpoint inhibitors could be useful in preventing respiratory dysfunction given the findings of exhausted T cell phenotype in COVID-19 patients [13] . However, overcoming T cell exhaustion risks exacerbating the hyperinflammatory state. The mixed clinical trial outcomes of anti-inflammatory therapy and the uncertain clinical risk-benefit of T cell related therapies collectively underscore a lack of resolution on the causative factors of severe COVID-19. Sieving the causative from the associative factors of pulmonary dysfunction thus remains an unmet need.

To define the principal drivers of pulmonary dysfunction in COVID-19 patients, we focused our attention on the host response to SARS-CoV-2 infection around the period of the nadir of respiratory function. We hypothesised that the principal drivers of COVID-19 severity must track tightly with respiratory function and must resolve in sync with respiratory improvement. Thus, instead of time interval analysis [14À16], we performed daily transcriptomic profiling to obtain a detailed continuum of the host immune response in COVID-19 patients before, during and after the nadir of respiratory function. We found that the neutrophil activation pathway, but not type-I IFN signaling pathways, showed a direct inverted relationship with respiratory function À genes in this pathway increased and peaked during the nadir of respiratory function and then declined during the recovery phase. On the other hand, transcripts associated with protein phosphorylation pathways, particularly the serine threonine kinases involved in cell cycle progression and T cell development increased at convalescence. The understanding of the temporal dynamics involved in severe COVID-19 progression and recovery provided insights into the potential drug candidates that can be considered to reduce the disease burden caused by SARS-CoV-2 infection. This study was approved by the SingHealth Combined Institutional Review Board (2013/397/F) and the patients were enrolled for the study at Singapore General Hospital. Healthy controls consisted of baseline samples obtained from participants of a measles, mumps and rubella re-vaccination study [17] . Ten healthy control individuals (7 females, 3 males) were included, with average age of 44.1 years (s.d.: 11.8, range: 30À61 years). All participants gave written informed consent, and study approval was obtained from the Sing-Health Combined Institutional Review Board (CIRB 2017/2374).

Patients hospitalised at Singapore General Hospital (SGH) testing PCR-positive for SARS-CoV-2 were approached to seek informed consent to participate in the study. Samples were collected from recruited patients between January 2020 to April 2020. Clinical data was obtained from the patients' medical records and clinical charts.

RNA was isolated from whole blood collected in Tempus Blood RNA tubes (Applied Biosystems Cat# 4342792) using the Tempus Spin RNA Isolation Kit (Life Technologies Cat# 4380204) according to manufacturer's instructions. Bioanalyzer assessment was performed on RNA samples before microarray was performed using the Affymetrix GeneChip Human Gene 2.0 ST Array (Applied Biosystems Cat# 902113) at the Duke-NUS Genome Biology Core Facility. Briefly, 50 ng of each RNA sample was processed using Affymetrix GeneChip Whole Transcript (WT) PLUS Reagent Kit (Applied Biosystems Cat# 902280), according to manufacturer's instructions. 15 mg of cRNA was used for the second cycle cDNA reaction and 5.5 mg of ss-cDNA was used for fragmentation. DNA fragment was end-labeled with biotin using terminal deoxynucleotidyl transferase before hybridisation on the arrays. Hybridisation cocktails containing fragmented, endlabeled cDNA were prepared and applied to GeneChip Human Gene 2.0 ST arrays. Hybridisation was performed at 60 rpm for 16h at 45°C using the FS450_0001 fluidics protocol. Arrays were scanned using Affymetrix GeneChip Command Console Software (AGCC) to produce .CEL intensity files. The raw .CEL intensity files and the processed files are deposited in ArrayExpress.

Evidence before this study

The COVID-19 pandemic has caused significant mortality and morbidity worldwide, but the principal drivers leading to severe disease progression and resolution remain poorly understood. Studies on severe COVID-19 patients have indicated that increased dysregulation and hyperactivation of the innate and acquired immune systems are associated with severity of the disease. However, it is unclear which of these responses are causative or merely associative.

Through daily transcriptomic profiling of severe and nonsevere COVID-19 patients spanning the nadir of respiratory function, this study reveals that neutrophil activation but not type I interferon signaling transcripts, tracked most closely with severe disease progression and resolution. This study thus highlights neutrophil activation as a critical factor of respiratory dysfunction in COVID-19.

Our findings can guide the future development of therapeutics against COVID-19. As neutrophil activation is functionally involved in COVID-19 resolution, drugs that target this pathway could be repurposed for the treatment of severe COVID-19. Moreover, our data may provide a valuable resource towards our understanding of the host-pathogen dynamics involved in severe COVID-19 disease progression and resolution.

NanoString profiling was performed using the nCounter Human Immunology v2 Panel (NanoString Technologies, Inc. Cat# XT-CSO-HIM2-12). 50 ng of RNA was mixed with 8 ul of reporter probeset and 2ul of capture probeset, and incubated at 65°C for 24 h. Hybridised samples were loaded onto an nCounter Sprint cartridge and imaged on a nCounter Sprint system. Raw and normalised counts were derived with the NanoString nSolver 4.0 software as per manufacturer's protocols. nCounter probe sets included negative and positive controls that were used for background thresholding, and normalizing samples for differences in sample input or hybridisation.

Plaque reduction neutralisation test (PRNT) was carried out on Vero E6 cells (RRID:CVCL_0574). Briefly, serial two-fold dilution of serum samples in DMEM media (Gibco Cat# 11995-040) containing 2% FCS (Hyclone Cat# SH30396.83) was incubated with 40 pfu of SARS-CoV-2 clinical isolate in equal volumes for 1 h. Thereafter, 100 mL of this serum/virus mix was added to Vero E6 cells in 24-well plates and incubated for 1 hour at 37°C, rocking every 15 min. After 1 h, serum/virus mix was removed and cells overlaid with 500ul of 0.8% methylcellulose (Calbiochem Cat# 17851) in DMEM media (DMEM, 2% FCS, 25 mM HEPES, 1x penicillin/streptomycin). After incubating at 37°C for 4 days, the cells were fixed with 4% formaldehyde for 30 min to inactivate the virus. Plates were rinsed with water and stained with 0.5% crystal violet for 5 min to visualise plaques. Excess crystal violet was removed and plates were rinsed with water, and allowed to dry before counting plaques.

The CIBERSORT tool was used to digitally quantify the relative frequencies of different immune cell types across time (https://ciber sort.stanford.edu) [18] . To assess if the upregulated proteins in the protein phosphorylation pathway had altered phosphorylation levels, genes from the protein phosphorylation pathway were cross-referenced with the phosphoproteome study described by Bouhaddou et al. [19] . Differential levels of phosphorylation modulated by SARS-CoV-2 infection in Vero E6 were extracted from the data analysis by Bouhaddou et al. For the analysis of the potential drugs that can target module A genes, we used the Drug Gene Interaction Database (DGIdb) to query the genes represented in module A. As clusters 3 and 4 are most associated with severe disease progression and resolution, the drugs that exclusively target only clusters 1 or 2 were filtered out. Subsequently, the top drug candidates were ranked based on the total number of genes targeted in module A [20] . The clinical usage information was obtained from Health Sciences Authority (HSA) at the website:

(https://eservice.hsa.gov.sg/prism/common/enquirepublic/Search DRBProduct.do?action=load&_ga=2.183810082.563179921. 1554083187-551332391.1551944793). Other information pertaining to the drugs such as route of administration and inflammatory properties were obtained via UptoDate and clue.io.

The Transcriptome Analysis Console (Thermo Fisher) was used to analyse transcriptomic changes between the different days relative to respiratory nadir function (Day 0). Differentially expressed genes (DEGs) were defined based on fold-change > 1.5 with respect to day 0, with false discovery rate adjusted p-value < 0.05 using the Benjamini-Hochberg Step-Up FDR-controlling Procedure (a=0.05). For pathway analysis, the identified DEGs were used as input data and analysed against the Gene Ontology (GO) Biological Processes database using the Enrichr tool [21] . Adjusted p-values for each enriched pathway were obtained from the Enrichr tool analysis, and REVIGO was used to summarise the non-redundant GO terms [22] . Volcano plots and heatmaps were constructed using the Prism 8.4.3 software.

Temporal gene expression was analysed using EDGE on log2 intensity counts [23] . A total of 4,831 genes were temporally altered in severe COVID-19 patients, p-value < 0.05 and q-value < 0.05. Coclustering of significant genes was analysed using Self-Organizing Maps (SOM) from Partek Ò Genomics Suite Ò based on the tabulated Least Square Means (LSMeans). 8 main clusters were identified which can be classified into module A and module B, representing 99.7% of the total significant genes. We excluded cluster 9 from further analysis as it contained only 13 genes. The genes in the different SOM clusters were then queried using the Enrichr tool to identify significantly enriched pathways (based on adjusted p-value < 0.05). To visualise the expression changes of pathway genes across time in severe and mild COVID-19 patients, normalised expression was tabulated by taking the average LSMeans values of all significant genes in the indicated pathways, such as cellular response to type I IFN pathway, regulation of immune effector process, neutrophil mediated immunity and protein phosphorylation pathways. Graphs and heatmaps were constructed using the Prism 8.4.3 software. To visualise the functional network and interacting partners involved in the protein phosphorylation pathway, the RegPhos database [24] and Ingenuity Pathway Analysis were used. For Principal Component Analysis (PCA), normalised gene transcript counts profiled using the Nano-String Human Immunology Panel v2 were used for analysis, where the R packages used are factoextra, ggplot2 and FactoMineR. Pearson correlation analysis was performed using Prism 8.4.3 software.

The raw data for microarray profiling is available at Array Express (E-MTAB-9721). The raw data and normalised counts for NanoString profiling of immune responses is available at Array Express (E-MTAB-9719).

The funders have no role in study design, data collection, data analyses, interpretation or writing of report.

We prospectively enrolled 10 SARS-CoV-2 confirmed patients admitted to the Singapore General Hospital after written informed consent. Full blood count, C-reactive protein (CRP) and lactate dehydrogenase (LDH) were measured as part of routine COVID-19 clinical management. These data are shown in Supplemental Figure 1a -d, with normal reference range for healthy subjects shaded in grey. The severe and moderately severe COVID-19 patients displayed lymphopenia despite normal white blood cell (WBC) counts, and exhibited higher CRP and LDH levels compared to the mild COVID-19 patients. These laboratory parameters were tracked longitudinally, and were consistent with previously reported associations with severe COVID-19 [25, 26] . Supplemental oxygen was administered if arterial oxygen saturation dropped below 95%. We used oxygen requirement to stratify patients into mild, moderately severe or severe COVID-19. Of the 10 patients, 4 were classified as mild as they did not require supplemental oxygen over the course of hospitalisation. Of the remaining 6 patients, all patients required supplemental oxygen; 5 patients needed more than 2 liters/minute of supplemental oxygen, requiring either non-invasive or invasive mechanical ventilation (MV) while the remaining patient received not more than 2 liters/minute of supplemental oxygen throughout hospitalisation (Supplemental Figure 1e ). All 6 patients were considered to have moderately severe or severe COVID-19, and were considered as a group in our analysis. In addition, patients who required supplemental oxygen were tracked using a 4-point ordinal scale, providing an overview of the disease progression and recovery over time (Fig. 1a) . For the 5 patients that needed non-invasive ventilation (NIV) or MV, the day of respiratory function nadir was defined as the day the patient was initiated on NIV or MV. For P6, the patient with moderately severe COVID-19, day of respiratory function nadir was defined as the day when arterial oxygen saturation was lowest at 93% (Supplemental Figure 1e ).

To interrogate the transcriptomic changes spanning the nadir of respiratory function nadir, we sampled daily, where possible, whole Waterfall plot representing the number of differentially upregulated and downregulated genes (Fold change > 1.5, adjusted p-value < 0.05; Benjamini-Hochberg step-up procedure) relative to nadir (day 0) in the severe patients (n=6). (d) Volcano plot displaying genes that were differentially expressed at day 7 post-nadir relative to day 0 in the severe patients (n = 6). Genes that were most differentially regulated were annotated on the volcano plot. (e) PCA of the NanoString data for healthy controls (black, n = 9), mild (blue, n = 3), moderately severe (orange, n = 1) and severe (red, n = 5) COVID-19 patients at the various days relative to respiratory function nadir. The top contributing genes in the respective PC is shown in Supplemental Figure 2b . blood from COVID-19 patients ranging from -4 days to 13 days, with the nadir being day 0 (Fig. 1b, Supplemental Figure 2a ). Healthy controls consisted of baseline samples obtained from participants of a measles, mumps and rubella re-vaccination study (CIRB 2017/2374), as previously reported [17] .

We first sought to attain an overview of whole genome expression changes before, during and after the nadir of respiratory function. Measurement of 29,405 annotated gene transcripts by microarray showed that the number of differentially expressed genes (DEGs) increased with time post-nadir (false discovery rate [FDR]-adjusted p < 0.05 and fold change [FC] > 1.5; Benjamini-Hochberg step-up procedure) (Fig. 1c) . Functional annotations of the DEGs revealed that cytokine and IFN-dependent immune responses, inflammatory signaling pathways, neutrophil mediated immunity and necrotic cell death pathways were most significantly reduced with time (Supplemental Figure 2a) . Notably, the most downregulated genes belonging to these pathways were IFIT1, IFI27, IFI44, IFI44L, CXCL10, RSAD2 and OAS3 (Fig. 1d) . In contrast, no consistent pathway was found among the upregulated genes, with the most upregulated genes being involved in metabolism and transport (RPIA, SLC14A1, SELENBP1), cytoskeleton (PLEK2, EPB42) and Notch signaling (TSPAN5) (Fig. 1d) .

Changes in expression of the innate and adaptive immune response genes were validated using the nCounter assay (NanoString Technologies), which allows simultaneous digital quantification of 579 immune-related mRNA transcripts. Mild COVID-19 patients and healthy controls were included in this analysis. Principal component analysis (PCA) confirmed that COVID-19 patients clearly segregated from healthy subjects in the PCA space (Fig. 1e) . Moreover, distinct temporal immune signatures between mild and severe patients could also be largely explained by the variance in principal component (PC) 1 (Fig. 1e) . The top contributors in PC1 were transcripts related to TLR signaling (TLR2, TLR4, TLR5, TLR8) and inflammatory responses (PRKCD, MAPK14) (Supplemental Figure 2b) . Critically, these transcripts were over-expressed in severe COVID-19 patients before or on the day of respiratory function nadir; these genes were universally downregulated after the nadir of respiratory function (Supplemental Fig. 2b ). On the other hand, PC3 segregated the healthy subjects from the infected patients. The top contributors are genes related to cytokine and IFN responses such as CXCL10, MX1, BST2, LAMP3 and IFI35, where higher levels of these transcripts were observed in both mild and severe COVID-19 patients, but not in healthy subjects (Supplemental Figure 2b ). Of note, PC2 did not segregate the healthy from the infected patients (Supplemental Fig. 2c-d) . Consistent with the microarray observations, the most significant immune pathways that showed reduced expression with time were the type-I IFN and TLR signaling pathways (Supplemental Figure 2e-f ). On the other hand, immune pathways that were increased during the recovery phase included pathways related to T cell proliferation and activation (Supplemental Figure 2e, g) . Collectively, these findings highlight that the suppression of the type-I IFN and TLR signaling pathways precedes restoration of the T cell proliferation and activation pathways.

To determine the genes that tracked with disease progression and resolution, we used extraction and analysis of differential gene expression (EDGE) to obtain greater insight on temporally regulated transcripts over the course of illness in severe patients. We identified 4,831 genes (FDR-adjusted p-value < 0.05; q-value < 0.05; Benjamini-Hochberg step-up procedure), that were significantly altered over time in the 6 severe COVID-19 patients (Supplemental Fig. 3a b). Self-organizing map (SOM) analysis of the 4,831 genes further revealed that these genes could be grouped into 8 distinct expression clustersÀ linear plots of expression of genes within each of these clusters against time are shown in Fig. 2a-b . These clusters can be broadly segregated into 2 main modules: Module A highlights the genes with reduced expression after respiratory nadir; module B highlights the genes with increased expression after respiratory nadir ( Fig. 2a-b) . In module A, clusters 1 and 2 showed a decreasing trend in gene expression from day -4 to day 13. Expression of genes in clusters 3 and 4 trended with respiratory dysfunction -expression increased from day -4 to day 0 and declined after day 0 during the recovery phase (Fig. 2a) . In module B, increased expression of the majority of the genes were observed after day 0 (Clusters 5-7), albeit at different rates in the different clusters (Fig. 2b) . The identity of the genes and heatmaps of the top 150 genes in each cluster are displayed in Supplemental Figure 3c -d.

We next examined the molecular pathways that were enriched in module A. Consistent with the observations from the nCounter assay, the top enriched pathway in cluster 1 is the cellular response to type I IFN (Fig. 2c) , where a monotonic decrease in gene expression was observed from day -4 to day 13 (Fig. 2d ). Cluster 1 contains 24 genes involved in the IFN response pathway, the majority of which are IFNstimulated genes (ISGs) (Supplemental Figure 3e) . Similar trends were also observed in mild COVID-19 patients although changes in the expression of these genes were smaller compared to severe cases (Supplemental Fig. 3f ). In cluster 2, the most enriched pathway is the regulation of immune effector process (Fig. 2c) , where the majority of these transcripts encodes for the variable and constant regions of the immunoglobulin heavy and light chains (Supplemental Figure 3g) . These genes were less abundant in whole blood from day -4 to day 0, increased transiently from day 0-4, and then declined thereafter in the severe patients (Fig. 2e) . In contrast, expression of these genes increased over time in mild patients (Supplemental Figure 3h) . Along with the differences in immunoglobulin transcript temporal profiles, severe cases developed higher titers of neutralizing antibodies to SARS-CoV-2 compared to mild cases (Supplemental Table 1 ). Taken collectively, our data suggests that the reduction in immunoglobulin transcripts may be explained by greater proportion of B cell migration to lymphoid organs in severe compared to mild cases.

The expression of genes that tracked most closely in temporal terms relative to disease progression and resolution were in clusters 3 and 4 (Fig. 2a) . Neutrophil mediated immunity and neutrophil activation were the top 2 enriched pathways in cluster 3 (Fig. 2c) . Other enriched pathways in cluster 3 were cellular response to mechanical stimulus, regulation of NFkB transcription factor activity, pattern recognition receptor signaling and IFN-gamma-mediated signaling pathways (Fig. 2c ). Cluster 4 only had neutrophil mediated immunity pathway being significantly modulated. Combining both clusters 3 and 4, we identified 87 genes involved in neutrophil mediated immunity that showed peak expression during the nadir of respiratory function; expression waned with upswing in respiratory function (Fig. 2f , Supplemental Figure 3i) . Critically, the expression of these genes remained unchanged throughout the period of monitoring in mild patients (Supplemental Figure 3j) . Collectively, these data support the notion that neutrophils may play a key role in shaping disease severity.

The abundance of genes identified in the SOM analysis could be influenced by changes in cell counts during acute and convalescent phase of infection. We thus analysed our data using CIBERSORT, which deconvolutes bulk gene expression data to enumerate immune cell type frequencies [27] . The plasma cell frequencies significantly correlated with the transcript levels of immunoglobulins (IGHV4-39, IGKV2D-40) from cluster 2, indicating that the reduced abundance of the immunoglobulin transcripts involved in immune effector processes reflected lower plasma cell frequencies (Supplemental Figure 4a ). In contrast, monocyte, neutrophil and naïve CD4 + T cell frequency levels did not correlate with temporally regulated genes in the IFN response, neutrophil activation and T cell activation pathways (Supplemental Figure 4b-d) , suggesting that transcriptional changes observed in these pathways in severe COVID-19 patients reflected the activation status of immune cells. Increased protein phosphorylation signatures after respiratory function nadir in severe COVID-19 patients Pathway enrichment on module B genes only led to significant enrichment of protein phosphorylation, particularly peptidyl-serine phosphorylation, in cluster 6 (Fig. 3a) . No significant pathways were found in the other clusters (adjusted p-value > 0.05; hypergeometric test). We identified 36 genes that were significantly upregulated in the protein phosphorylation pathway, of which 20 were serine-threonine kinases (Fig. 3b-c) . When these genes were cross-referenced against RegPhos, which allows for exploration of the human protein phosphorylation regulatory network [24] , many of these kinases were identified as interacting partners (Supplemental Figure 4e ) in the phosphatidylinositol 3-kinase (PI3K)/protein kinase B (AKT) pathway. Critically, these kinases are involved in cell growth, cell survival and T cell development (Fig. 3d) , which could be crucial for renewal and regeneration of immune cells after viral infection. We also found that 12 of the 36 genes (33.3%) were differentially phosphorylated during SARS-CoV-2 infection (Supplemental Table 2 ). Of note, genes associated with the protein phosphorylation pathway were similarly increased in mild patients over time, albeit to a lesser extent as compared to severe patients (Supplemental Figure 4f) . Collectively, increased expression of serine-threonine kinases over time suggests a protective role for these kinases in COVID-19 disease resolution.

Our analyses thus far suggest that genes in clusters 3 and 4 that are involved in neutrophil activation, pattern recognition receptor signaling and NFkB-mediated responses are functionally involved in COVID-19 resolution, as these genes are most strongly associated with severe COVID-19 disease progression and resolution (Fig. 2a, c) . To validate this notion, we cross-referenced the transcripts in clusters 3 and 4 with the Drug Gene Interaction Database (DGIdb), which is a curated resource that can be used to search for druggable targets of the human genome [20] . If indeed our analysis was clinically pertinent, we should identify therapies with demonstrated efficacy in reducing COVID-19 severity. One such drug is dexamethasone. This corticosteroid has been found to reduce mortality rates in severe COVID-19 patients [28] . We found that amongst the module A genes, 6 were known targets of dexamethasone. Of the 6 genes targeted by dexamethasone, 2 were in clusters 3 and 4 (Fig. 4) . Specifically, dexamethasone was found to act on JUNB and IRS2 (insulin receptor substrate 2) of clusters 3 and 4, respectively (Supplemental Table 3 ). Both these genes are involved in the cytokine-mediated signaling pathway and are most highly expressed in neutrophils compared to other blood cells [29] . Notably, aspirin was identified as the top drug that inhibits the most genes in clusters 3 and 4, which corroborates a recent cohort study showing that aspirin use is associated with less mechanical ventilation, ICU admissions and mortality in hospitalised COVID-19 patients [30] . In contrast, host targets of a different corticosteroid, methylprednisolone [12, 31] , and other drugs such as azithromycin, lopinavir-ritonavir and hydroxychloroquine [31À34], all of which have shown no efficacy in improving clinical outcomes in COVID-19 patients, were absent from the hits in clusters 3 and 4. Likewise, host-directed anti-IL-6 therapies such as tocilizumab and sarilumab were also not highlighted in this cross referencing. Our findings are thus consistent with previously published efficacy trials for severe COVID-19.

Given how we were able to identify targets of dexamethasone but not other host-directed therapies in our analysis, we next sought to identify other drugs that may be suitable for repurposing as anti-COVID-19 therapies ( Fig. 4 ; Supplemental Table 3 ). Moreover, as the cluster 6 genes from module B consist of serine-threonine kinases that could be important for convalescence, we also highlighted the drug candidates that inhibited module A genes without known activity on cluster 6 genes (Fig. 4) . Many of these drugs are antiinflammatory, such as aspirin, sulfasalazine, celecoxib and indomethacin. Interestingly, several statins such as lovastatin, atorvastatin and simvastatin can specifically target module A genes without affecting the expression of genes from cluster 6. The full list of the potential drug candidates and their corresponding gene targets are detailed in Supplemental Table 3 .

Longitudinal sampling of blood and other samples from COVID-19 patients has enabled the identification of multiple associative factors Drugs that also target cluster 6 genes are indicated in blue, whereas drugs do not interfere with cluster 6 genes are indicated in purple. The genes that the drugs target are shown in Supplemental Table 3. of respiratory dysfunction. However, associative factors that could also drive pulmonary dysfunction in COVID-19 patients have not been systematically defined. A major limitation to identifying drivers of SARS-CoV-2 infection induced pulmonary dysfunction is the lack of an animal model that accurately recapitulates COVID-19 pathogenesis in humans. To overcome this hurdle, we postulated that the host response to SARS-CoV-2 infection around the period of the nadir of respiratory function could provide us with a window into distinguishing the pathogenic from the associative factors of respiratory dysfunction. The principal drivers of COVID-19 severity should resolve before or track tightly with respiratory function.

To test our postulate, we coupled daily sampling with full genome gene expression measurement in combination with EDGE and SOM analyses. We identified 8 distinct gene expression clusters that were altered over the course of disease in severe patients. Cluster 1 and 2 contained transcripts, many of which are IFN-stimulated genes (ISGs) [35] and other immune effectors, all of which showed a sustained decrease in gene expression starting from day -4 relative to the nadir of respiratory function. That these pathways did not track tightly with respiratory function suggest that the pro-inflammatory and other immune responses may initiate the cascade of processes necessary for disease pathogenesis but are not the principal drivers of respiratory dysfunction.

Transcripts that tracked most closely with progression and resolution of respiratory dysfunction were in clusters 3 and 4. These transcripts belong to neutrophil-mediated activation, pattern recognition receptor signaling and NF-kB activation pathways. They were significantly enriched in severe COVID-19 patients but not in patients that did not develop progressive pulmonary deterioration. As the bone marrow is the major site of neutrophil production, the increase in cytokine and chemokine levels before respiratory function nadir (cluster 1 and 2) could have led to an increased neutrophil-to-lymphocyte ratio in severe COVID-19 patients [36] , leading to the accumulation of activated neutrophils in the lung. Our finding is congruent with previous studies showing elevated neutrophil infiltration, activation and neutrophil extracellular trap (NET) markers in the lung of COVID-19 patients at post mortem examination [37, 38] . Similarly, RNA-seq analysis of bronchoalveolar lavage fluid (BALF) from COVID-19 patients identified an upregulation of neutrophil enriched and chemotaxis genes (e.g. CD177, FPR2, CXCR1, CXCR2, CTSZ, S100A8, S100A9) in COVID-19 BALF cells [39, 40] . Of note, S100A8 and S100A9, which were identified in cluster 3 (Supplemental Figure 3i ), are also known to form a stable heterodimer during inflammation known as calprotectin [41] À plasma calprotectin levels have been found to discriminate between mild and severe COVID-19 patients [42] . Besides COVID-19, neutrophils are also enriched in the lungs of patients with ARDS caused by other viral infections [43] . For example, pathogenic influenza viral strains can cause increased neutrophil mobilisation and recruitment to the lungs in mice and humans, that can result in excessive inflammation and fatality [44, 45] . Taken together, neutrophil activation may thus be an important underpinning of respiratory dysfunction in COVID-19 patients.

An additional new finding was the increased transcript levels of serine-threonine kinases (cluster 6) during the recovery from severe COVID-19, raising the possibility that these kinases play an important role in convalescence. Interestingly, many of these are serine-threonine kinases involved in cell cycle progression, which is congruent with the phosphoproteomics analysis conducted by Bouhaddou et al that demonstrated that SARS-CoV-2 infection in Vero E6 cells influences the phosphorylation levels of cell cycle proteins such as PDPK1, PRKACA, RAF1 and RPS6KA3 (Supplemental Table 2 ) [19] . Increased expression and activation of these genes during recovery could facilitate cell cycle progression and aid in restoring the frequencies of immune cell populations, particularly the T cells, that are decreased over the course of the disease [4, 6, 37, 38] . This explanation is further supported by the observed increase in T cell transcripts starting from day~6 after respiratory function nadir. Thus, although the transcripts in cluster 6 did not track with progression of respiratory dysfunction, the coincidence of their expression in respiratory function recovery suggests that T cell proliferation and development may be important for convalescence.

The tracking of genes in clusters 3 and 4 with respiratory dysfunction are consistent with multiple transcriptomics and proteomics studies in COVID-19 patients [46À49] (Supplemental Table 4 ). That we observed these genes to be significantly altered in a high-resolution time series analysis further highlights the potential for these genes to be targeted for treating severe COVID-19. Indeed, the finding of dexamethasone [12] and aspirin [30] when we cross-referenced our list with DGIdb, suggests the clinical validity of this approach. We also identified sulfasalazine, an anti-inflammatory drug, as another potential drug candidate that targets the most number of genes from clusters 3 and 4. Interestingly, the third candidate, methotrexate is an FDA-approved inhibitor of purine synthesis that can inhibit SARS-CoV-2 replication [50] . Clinical trials to test the translatability of these drug candidates could thus be particularly fruitful.

Our present study has some limitations. While this study provides high resolution on the temporal dynamics of the host response involved in disease progression and recovery, gene expression was profiled from peripheral blood rather than infected organs. Direct assessment of gene expression in infected organs would require invasive procedures like bronchoalveolar lavage, which is ethically challenging for mild COVID-19 cases, and risks aerosolization of infectious respiratory fluids that could compromise the safety of healthcare workers. Measuring gene expression in peripheral blood hence provides an opportunity to understand how the host response to SARS-CoV-2 could impact clinical outcome. Furthermore, our findings were based on a small sample size involving 10 COVID-19 patients, of which 8 were male (Supplemental Table 1 ). As males and females may respond to COVID-19 differently [51À53], these temporal host responses may not be fully generalisable to females. As such, increasing the sample size, as well as a deeper temporal characterisation of neutrophil gene expression changes, phenotype and function at the target sites of infection could provide further mechanistic insights into the molecular processes that drive disease pathogenesis and recovery from severe COVID-19.

In conclusion, our results reveal the host responses involved in the progression and resolution of respiratory dysfunction in COVID-19 patients, which could be useful in guiding the choice of drugs for repurposing as COVID-19 therapeutics. 

Data in this study is available upon request from the corresponding author at kuanrong.chan@duke-nus.edu.sg. Individual participant data will be made available, study protocol, statistical analysis plan, the raw and processed transcriptomics data for all participants will be made available for anyone who wishes to access the data. The data can be used for any purpose and the data will be made available indefinitely in the ArrayExpress and Expression website:

http://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-9721 http://wwwdev.ebi.ac.uk/gxa/experiments/E-MTAB-9721?access Key=edad7f29-fb8e-4a51-87d3-661176fa6e5e

Dr Eng Eong Ooi reports grants from National Medical Research Council, Singapore, during the conduct of the study; personal fees from Tychan Private Limited, outside the submitted work; the other authors declare no competing interests.

Intensive care management of coronavirus disease 2019 (COVID-19): challenges and recommendations

Complex immune dysregulation in COVID-19 patients with severe respiratory failure

Clinical features of patients infected with 2019 novel coronavirus in Wuhan

Deep immune profiling of COVID-19 patients reveals distinct immunotypes with therapeutic implications

Heightened Innate Immune Responses in the Respiratory Tract of COVID-19 Patients

A dynamic immune response shapes COVID-19 progression

Interleukin-6 blockade with sarilumab in severe COVID-19 pneumonia with systemic hyperinflammation: an open-label cohort study

Tocilizumab in patients with severe COVID-19: a retrospective cohort study

Retrospective multicenter cohort study shows early interferon therapy is associated with favorable clinical responses in COVID-19 patients

Early IL-1 receptor blockade in severe inflammatory respiratory failure complicating COVID-19

Trials of anti-tumour necrosis factor therapy for COVID-19 are urgently needed

Association between administration of systemic corticosteroids and mortality among critically Ill patients with COVID-19: a meta-analysis

Marked T cell activation, senescence, exhaustion and skewing towards TH17 in patients with COVID-19 pneumonia

Systems biological assessment of immunity to mild versus severe COVID-19 infection in humans

A dynamic COVID-19 immune signature includes associations with poor prognosis

Longitudinal analyses reveal immunological misfiring in severe COVID-19

Genomic signature of early T-cell response is associated with lower antibody titer threshold for sterilizing immunity

Robust enumeration of cell subsets from tissue expression profiles

The global phosphorylation landscape of SARS-CoV-2 infection

DGIdb 3.0: a redesign and expansion of the drug-gene interaction database

Enrichr: a comprehensive gene set enrichment analysis web server 2016 update

REVIGO summarizes and visualizes long lists of gene ontology terms

Significance analysis of time course microarray experiments

RegPhos: a system to explore the protein kinase-substrate phosphorylation network in humans

Lactate dehydrogenase and C-reactive protein as predictors of respiratory failure in CoVID-19 patients

Lymphopenia predicts disease severity of COVID-19: a descriptive and predictive study

Determining cell type abundance and expression from bulk tissues with digital cytometry

Dexamethasone in hospitalized patients with COVID-19 -preliminary report

A genome-wide transcriptomic analysis of protein-coding genes in human blood cells

Aspirin use is associated with decreased mechanical ventilation, icu admission, and inhospital mortality in hospitalized patients with COVID-19

Methylprednisolone as adjunctive therapy for patients hospitalized with COVID-19 (Metcovid): a randomised, double-blind, phase IIb, placebo-controlled trial

Azithromycin in addition to standard of care versus standard of care alone in the treatment of patients admitted to the hospital with severe COVID-19 in Brazil (COALITION II): a randomised clinical trial

Lopinavir-ritonavir in patients admitted to hospital with COVID-19 (RECOVERY): a randomised, controlled, open-label, platform trial

Chloroquine does not inhibit infection of human lung cells with SARS-CoV-2

Interferon-stimulated genes and their antiviral effector functions

Longitudinal characteristics of lymphocyte responses and cytokine profiles in the peripheral blood of SARS-CoV-2 infected patients

Targeting potential drivers of COVID-19: neutrophil extracellular traps

Neutrophil extracellular traps in COVID-19

Excessive neutrophils and neutrophil extracellular traps in COVID-19

Transcriptomic characteristics of bronchoalveolar lavage fluid and peripheral blood mononuclear cells in COVID-19 patients

S100A8/A9 in inflammation

Elevated calprotectin and abnormal myeloid cell subsets discriminate severe from mild COVID-19

A role for neutrophils in viral respiratory disease

Neutrophils-related host factors associated with severe disease and fatality in patients with influenza infection

The influenza virus protein PB1-F2 increases viral pathogenesis through neutrophil recruitment and NK cells inhibition

Immune and metabolic signatures of COVID-19 revealed by transcriptomics data reuse

Transcriptional and proteomic insights into the host response in fatal COVID-19 cases

Transcriptome-based drug repositioning for coronavirus disease 2019 (COVID-19)

The transcriptomic profiling of SARS-CoV-2 compared to SARS, MERS, EBOV, and H1N1

Methotrexate inhibits SARS-CoV-2 virus replication "in vitro

Sex differences in immune responses that underlie COVID-19 disease outcomes

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

Gender differences in patients with COVID-19: focus on severity and mortality

We thank all doctors and nurses in the isolation facility in the Singapore General Hospital. ViREMiCS was established through generous support from the Tanoto Foundation. This study was sponsored in part by a generous gift from The Hour Glass. EEO and JGL are funded by the National Medical Research Council of Singapore, through the Clinician Scientist Awards.

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.ebiom.2021.103262.