key: cord-0966952-shpx8nc5
authors: Yan, Q.; Li, P.; Ye, X.; Huang, X.; Mo, X.; Wang, Q.; Zhang, Y.; Luo, K.; Chen, Z.; Luo, J.; Niu, X.; Feng, Y.; Ji, T.; Feng, B.; Wang, J.; Li, F.; Zhang, F.; Feng, L.; Lei, C.; QU, L.; Chen, L.
title: Longitudinal peripheral blood transcriptional analysis of COVID-19 patients captures disease progression and reveals potential biomarkers
date: 2020-05-08
journal: nan
DOI: 10.1101/2020.05.05.20091355
sha: 2a625f987cb85d04493003f1607ae082ac822850
doc_id: 966952
cord_uid: shpx8nc5

COVID-19, caused by SARS-CoV-2, is an acute self-resolving disease in most of the patients, but some patients can develop a severe illness or even death. To characterize the host responses and identify potential biomarkers during disease progression, we performed a longitudinal transcriptome analysis for peripheral blood mononuclear cells (PBMCs) collected from 4 COVID-19 patients at 4 different time points from symptom onset to recovery. We found that PBMCs at different COVID-19 disease stages exhibited unique transcriptome characteristics. SARS-CoV-2 infection dysregulated innate immunity especially type I interferon response as well as the disturbed release of inflammatory cytokines and lipid mediators, and an aberrant increase of low-density neutrophils may cause tissue damage. Activation of cell death, exhaustion and migratory pathways may lead to the reduction of lymphocytes and dysfunction of adaptive immunity. COVID-19 induced hypoxia may exacerbate disorders in blood coagulation. Based on our analysis, we proposed a set of potential biomarkers for monitoring disease progression and predicting the risk of severity.

The recent global outbreak of COVID-19 is caused by a highly contagious new 40 coronavirus named SARS-CoV-2 (1-3). WHO has declared the outbreak of 41 COVID-19 a Public Health Emergency of International Concern (4). As of April 28, 42 there have been more than 3 million confirmed cases and more than 210,000 deaths 43 worldwide according to the reports of WHO. Most patients with COVID-19 showed 44 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 5 89 Principal component analysis (PCA) of gene expression grouped the patient samples 90 into 4 clusters. Interestingly, these 4 clusters coincided with disease progression in 91 clinical observation (Fig. 1B, Fig. S1A ). We named these four clusters as: 1) Stage 1 92 (S1), representing the early onset; 2) Stage 2 (S2), representing the clinically most 93 severe stage; 3) Stage 3 (S3), representing improving stage; and 4) Stage R, 94 representing the recovering or convalescent stage. Notably, all three patients in S2 95 (PtQ5, PtJ7, PtL9, PtW9) were at the most severe disease state, demonstrated by the 96 highest C-reactive protein in plasma, the lowest lymphocyte counts, and the worst 97 chest radiography (Fig. S1A ). Four samples from the healthy donor formed a distinct 98 cluster, which was named as cluster H. Of note, cluster R is adjacent to cluster H, 99 suggesting the transcriptomes in the recovery stage is approaching the healthy state. 100 This result demonstrates that the gene expression profiles have distinct patterns along 101 with disease progression. These patterns were reproducible in different COVID-9 102 patients, at least for non-ICU patients.

We performed a digital cytometry CIBRSORTx (15) to delineate the transcriptome 105 into abundances of cell subsets (Fig. S1D) in the PBMCs at different stages of disease 106 progression. This analysis showed a dramatic increase of monocytes and pathological 107 low-density neutrophils in peripheral blood during S1 and S2 (Fig. S1D ). In contrast, 108 there was a reduction of T and NK cells in S1 and especially S2. The perturbation in 109 the proportion of cell types in PBMCs, especially lymphocytopenia is one of the 110 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 6 leading clinical manifestations of COVID-19. Our analysis was verified by clinical 111 tests of T lymphocyte and neutrophil counts in these patients (Fig. S1A) . To identify the differentially expressed genes (DEGs) at different disease stages, we 116 first performed paired comparison between S1 and R, S2 and R, S3 and R for each 117 person (Fig. 1C) . We also combined all samples from different patients at the same 118 disease stage for group comparison (Fig. 1D ). Compared to the convalescent state, the 119 number of DEGs was 2661, 4811, and 834 in S1, S2, and S3 respectively, 120 representing 8.7%, 15.8%, and 2.7% genes in the transcriptomes (Fig. 1D ). In contrast, 121 there were few changes in the number of DEGs before and after vaccination in a 122 healthy person. We found that there were many identical DEGs among different 123 COVID-19 patients, especially when these samples belonged to the same disease 124 stage. The 4 patients shared 155, 341 and 23 DEGs at S1, S2 and S3 respectively (Fig. 125 S1E-F). Therefore, SARS-CoV-2 infection resulted in common changes in gene 126 expression profile in different regular COVID-19 patients. Since the samples in R 127 may not fully represent the healthy state, we also performed grouped comparison 128 between S1 and H, S2 and H, S3 and H, R and H (Fig. S1G) . Indeed, there were more 129 DEGs when compared with gene expression in a healthy state. An average of 19.3%, 130 17.1%, 11.0%, and 3.9% genes in the transcriptome have undergone significant 131 expression changes in S1, S2, S3, and R, respectively. 132 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. To identify genes closely related to disease progression, we performed a time series 134 analysis for all DEGs. The healthy samples were assumed to be a stage representing 135 pre-infection to capture the time-dependent changes. All DEGs could be divided into 136 6 clusters based on their expression patterns (Fig. 1E ). The pattern of gene expression 137 in each cluster somewhat appeared to resemble an alphabet, so we named each cluster 138 as "n", "A", "M", "U", "V", and "W". Each cluster contains 1659, 1668, 1717, 1780, 139 1703, and 1599 genes respectively (Fig. 1F ). To investigate the biological functions 140 associated with each cluster, we performed gene ontology analysis. The top 5 biology 141 process (BP) terms of each cluster were listed (Fig. 1G) . Cluster "n", "A", and "M" 142 contained 1043，609, and 339 BP terms respectively (q-value <0.05). Genes in cluster 143 "U", "V" and "W" showed few discernible biological themes; each contained only 23, 144 11, and 0 BP terms (Fig. 1G, Fig. S1H ). To verify the contribution of the genes in 145 cluster "n", "A" and "M" to disease progression, we took the top 2000 genes with the 146 largest contribution to PC1 and PC2 (termed gene set PC1 and PC2 hereafter) and 147 compared with genes in the cluster "n", "A" and "M", respectively (Fig. 1B, Fig. S1I ).

The result revealed that genes in cluster "n" highly overlapped with gene set PC2, 149 which reflects the difference between illness and healthy state. Genes in cluster "M" 150 highly overlapped with gene set PC1, which reflects disease severity. Genes in cluster 151 "A" partially overlapped with both PC1 and PC2. Therefore, we mostly focused on 152 these three clusters (n, A and M) that were likely associated with disease progression 153 in the subsequent analysis.

154 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 8 Inflammatory cytokines and lipid mediators are disturbed in COVID-19 patients 155 We noticed that a set of DEGs in cluster "n" were enriched in the biological processes: 156 immune effector process, response to cytokine, inflammatory response, cytokine 157 production, and cell chemotaxis ( Fig. 2A-B) . Their expression increased in S1 and S2, 158 declined in S3 and R (Fig. 2C, Fig. S2A ). These "n" type genes include: IL1B, IL6, 159 IL10, IL12A, IL18, IL19, IL27, and chemokine CXCL10 (IP-10), are known to 160 promote the inflammatory response. Persistent high-level expression of chemokines 161 CXCL1, CXCL2, CXCL3, and CXCL8 promote neutrophil infiltration into the lung 162 and/or other tissues, while CCL2 promotes monocyte migration to the site of infection 163 (Fig. 2C, Fig. S2A ). Interestingly, we found another set of genes in cluster "A", 164 including cytokines such as IL2, IL7, IL15, and IL21, their expression sharply 165 elevated in S2 and returned to near-normal levels, as clinical symptoms resolved in S3 166 and R (Fig. 2C, Fig. S2A ). These cytokines are associated with the proliferation and 167 activation of T, B, and NK cells. To verify our analysis, we used Luminex to confirm 168 the protein levels of some cytokines in the plasma of COVID-19 patients in this study.

IP-10 and IL-6 were indeed significantly elevated in S1 and S2, but returned to near 170 normal range in S3 and R as the disease starting to resolve (Fig. 2D ). These findings 171 further indicated that dysregulated cytokines release strongly associated with the 172 progression of the disease.

We found that lipid mediators also play an important role in the inflammatory 175 response of COVID-19 patients. The expression pattern of enzymes involved in lipid 176 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 mediator synthesis also belonged to cluster "n", including cyclooxygenases and 177 lipoxygenases which convert arachidonic acid (AA) to lipid mediator prostanoids.

The increased expression of TLR7 and NFKB2 promoted the expression of 179 PLA2G4A and PTGS2 (COX-2) in S1 (Fig. 2E) . The key enzymes of lipid mediator 180 synthesis, including PLA2G7, PLA2G15 and PTGES3 (PGE2 synthase), also 181 increased in S1 and S2 ( Fig. 2E-F, Fig. S2B ). Elevation of these enzymes can promote 182 the production of inflammatory prostanoids and inflammatory response. Interestingly, 183 lipoxygenase ALOX15B, which produces anti-inflammatory mediator LXA4, had a 184 surge in S1 as an initial response, but declined in S2 and resumed to near healthy state 185 in R as the disease resolved ( Fig. 2E-F, Fig. S2B ). We measured the plasma 186 concentration of LXA4 in these patients and found that LXA4 also increased 187 significantly in S1 and then declined in S2 and thereafter (Fig. 2G) . Therefore, lipid 188 mediators are likely to play an important role in inflammation-associated disease 189 severity.

Elevated pathological low-density neutrophils may contribute to tissue damage 192 We found that neutrophil-associated genes were highly enriched in cluster "n" (Fig. 193 3A-B). During viral infections or autoimmune diseases, neutrophils may abnormally 194 differentiate to pathological low-density neutrophils (LDNs) that remain in PBMCs 195 after Ficoll-gradient centrifugation (16, 17). Our digital cytometry showed an increase 196 of low-density neutrophils during S1 and S2 (Fig. S1D, Fig. 3D ). Besides, complete 197 blood count verified that there was an elevated neutrophils level in S2 in these 198 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020. granulocytes (LRG1). The expression of these genes significantly increased in S1 and 202 S2, and gradually returned to near normal in S3 and R (Fig. 3C, Fig. S3A ). These of NETs formation in COVID-19 increased in S1, S2, and S3. As clinical symptoms 211 alleviated in R, these genes declined to near-normal levels (Fig. 3C, Fig. S3A ). The "n" 212 cluster also contains genes related to neutrophil chemotaxis and migration. CXCR1, 213 CXCR2, FPR1, FPR2, CD177 and AQP9 are known to regulate neutrophil 214 chemotaxis, which also increased in S1 and S2 (Fig. 3C, Fig. S3A ). Furthermore, the 215 transcripts of inflammatory-related receptors (GPR84, IL1R2) also elevated in S1, S2, 216 and S3 (Fig. 3C, Fig. S3A ). CXCL8 (IL-8) is secreted primarily by neutrophils. Its 217 elevation in S1 and S2 (Fig. 3F ) may acts as a chemotactic factor by further recruiting 218 the neutrophils to the site of infection. We measured the plasma concentration of IL-8 219 in these patients and found that plasma IL-8 started to increase in S1, peaked in S2, 220 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 decreased in S3, and returned to normal level in R (Fig. 3G) . Therefore, the increase 221 of LDNs may contribute to pneumonia and tissue damage.

The interferon responses are imbalanced in COVID-19 patients 224 We found that the genes involved in innate immunity process were enriched in cluster IFNB (IFN-β), IFNE (IFN-ε) and IFNW (IFN-ω) was almost undetectable from S1 to 230 R (Fig. 3K ). In contrast, another subtype of type I IFNs, IFNK (IFN-κ), and type II 231 IFNs, IFNG (IFN-γ), started to increase in S1, peaked in S2, decreased in S3 and 232 reached the normal level in R (Fig. 3K ). While type III IFNs, IFNL (IFN-λ) was 233 unchanged from S1 to R (Fig. 3K ). These results indicate that the expression of 234 interferons is dysregulated in COVID-19 patients.

Nevertheless, there was a significant increase in type I interferon-stimulated genes 237 (ISGs) in S1, including ISG15, ISG20, TRIM5, TRIM25, and APOBEC3A, 238 IFN-induced GTP-binding protein (MX1 and MX2), 2'-5'-oligoadenylate synthase 239 (OAS1, OAS2, OAS3, and OASL), and interferon-inducible protein (IFI6, IFI27, 240 IFI35, and IFITM3). These ISGs decreased in S2, and gradually reduced to normal 241 level in S3 and R (Fig. 3J, Fig. S3B ). These results suggested that the antiviral 242 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020. . https://doi.org/10.1101/2020.05.05.20091355 doi: medRxiv preprint response was significantly enhanced in early infection. Interestingly, the IFN-γ 243 inducible genes CCL2, CXCL10 showed similar expression patterns as IFNG, which 244 increased in S1-S2, decreased in S3-R (Fig. 2C , 3J-K). These results demonstrated 245 that the significant elevation of type II interferon and type II interferon inducible 246 cytokines might be a hallmark that the disease is in the most severe stage. in S1, declined to the lowest level in S2, then elevated in S3 and slightly decreased 254 again in R (Fig. 4C, Fig. S4A ). Expression of TCR was at the lowest level in S2, and 255 rebounded in S3-R (Fig. 4D) . Interestingly, the decrease of TCR diversity occurred in 256 S1, remain low in S2, increased in S3 and R, but was still much lower than the healthy 257 state ( Fig. 4E ), suggesting T-cells started to return to circulation in S3 and R. Of note, 258 GZMB (granzyme B) increased in S1, but decreased in S2, then increased again in antibody-secreting cells (ASCs), CD38 and CD27 elevated in S1, reduced in S2, 264 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020. . https://doi.org/10.1101/2020.05.05.20091355 doi: medRxiv preprint 13 increased in S3, and then decreased in R (Fig. 4F, Fig. S4B ). BCR expression 265 increased in S1, decreased in S2, increased in S3, then decreased again in R (Fig. 4G ).

BCR diversity was at the lowest level in S2, elevated in S3-R, and peaked in R ( Activation of cell death, exhaustion, and migration contribute to a reduction of 279 lymphocytes and NK cells 280 We found that genes involved in cell death pathway were enriched in several clusters 281 (Fig. 5A, B) . Notably, many pro-apoptotic molecules were enriched in cluster "M", 282 which was similar to T and B cell-related genes. TP53 (p53), BAX, BAK1, BAD, 283 BIK increased in S1, decreased in S2, elevated in S3, and then decreased to normal 284 level in R. The extrinsic apoptosis, necrotosis, and pyroptosis related genes 285 TNFRSF12A (TWEAK), TNFSF10 (TRAIL), FADD, CASP3, CASP7, CASP9, 286 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020. . https://doi.org/10. 1101 14 CYCS (cytochrome c), RIPK1, MLKL, and GSDME, also elevated in S1 or S2 ( Fig.   287 5C, Fig. S5A ). In contrast, the anti-apoptotic gene BCL2 was low in S1-S2, but 288 increased in S3-R. Therefore, the induction of cell death molecules during early onset 289 might be associated with lymphocyte and NK cells death in S2. We found that ARG1 290 (arginase) had an increase in S2 (Fig. 5D ). Arginase is constitutively expressed by 291 neutrophils and is an enzyme that metabolizes arginine into ornithine and urea. The Interestingly, the S1PR1, S1PR2 and S1PR4 and S1PR5 which are associated with 298 lymphocyte and NK cells egression from peripheral lymphoid organs decreased to the 299 lowest level in S2, while the S1PR3 that mostly expresses in monocytes increased in 

308 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020. . https://doi.org/10.1101/2020.05.05.20091355 doi: medRxiv preprint by hypoxia 309 We observed that blood coagulation related genes were significantly enriched in 310 cluster "n", including blood coagulation, platelet activation, plasminogen activation, 311 blood vessel development, fibrin clot formation (Fig. 6A-B) . Coagulation factors (F3, 312 F5, F12, and F13A1) and platelet-activating genes (GP1BA, GP1BB, and GP9) were 313 elevated in S1 and maintained at a high level until S3. Several transcripts related to 314 the process of fibrin dissolution, including PLAU (urokinase plasminogen activator, 315 uPA) and PLAUR (uPA receptor), increased during S1-S3 ( Fig. 6C-D, Fig. S6 ).

PLAU and PLAUR acts to convert plasminogen to plasmin, which promotes 317 fibrinolysis, leading to generation of D-dimers (fibrin degradation products). D-dimer 318 is a known hemostasis marker that reflects ongoing fibrin formation and degradation.

In clinical tests, the 4 patients in the study had average D-dimer level at 725 ± 371 320 μg/L in S1, which was much higher than that in healthy people (45 ± 23 μg/L) (Fig.   321 6E). Paradoxically, SERPINE1 which encodes plasminogen activator inhibitor type 322 1(PAI1), also elevated during S1-S3. PAI1 inhibit tissue-type plasminogen activator 323 (tPA) and uPA. Moreover, additional serpin family members (SERPINA1, SERPINB2 324 and SERPING1) also showed elevated expression during S1-S3 (Fig. 6C, Fig. S6 ). 325 We speculate that PAI1 may inhibit uPA and tPA mediated fibrinolysis and stimulate 326 thrombus formation in COVID-19 patients. These data suggest that coagulation 327 disorder occurs early and lasts most period of COVID-19 disease progression.

We found a subset of genes involved in hypoxia response was enriched in cluster "A" 330 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020. . https://doi.org/10.1101/2020.05.05.20091355 doi: medRxiv preprint ( Fig. 6A-B ). HIF1A and EPAS1 (HIF2A) are the key genes responding to low oxygen, 331 were significantly elevated in S2, although the hypoxia was likely occurred in S1 (Fig.   332 6C, Fig. S6 ). Its downstream genes, including VEGFA (vascular endothelial growth 333 factor A), EGF (epidermal growth factor), ANGPT1 (angiopoietin1), TIMP1 (TIMP 334 metallopeptidase inhibitor1), EDN1 (Endothelin1), and TFRC (transferrin receptor) 335 increased in S2 and/or S3 (Fig. 6C, Fig. S6 ). The upregulation of these genes 336 contributes to angiogenesis, accelerated iron metabolism, and vascular contraction, 337 which can help to improve oxygen delivery. Indeed, clinical observation of the 338 patients in this study showed their arterial oxygen partial pressures started to decline 339 in S1 and were the lowest in S2 (Fig. 6F) . 

350 progression and prognosis 351 We applied multi-category classification (MCC) based on logistic model to identify 352 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020. analysis using the combination of these 25 genes was sufficient to define different 361 disease stages without using the whole transcriptomes (Fig. 1B, Fig. S7B ).

Finally, we attempted to further identify biomarkers that may predict the risk of severe 364 disease during the early onset (S1), preferably before the change occurs for traditional 365 biomarker CRP (which was elevated during S2 in our patients). We performed another 366 RNA-seq analysis for PBMCs collected from 3 severe cases during their early stages 367 (day 4, 6, and 8 after symptom onset, respectively) before these patients were 368 transferred to ICU. 12 genes showed good predictive power of disease severity prior 369 to or parallel to clinically severe manifestation were revealed, including IL6, IL10, 370 CXCL2, ALOX15B, IL1R2, MMP9, GPR84, F3, SERPINE1, CD8B, CD3G, and 371 S1PR1 (Fig. 7A ). AUC analysis based on the expression of these 12 genes showed an 372 excellent prediction of the severity (Fig. 7B) . This result demonstrated that the 373 combinatorial analysis of 12 genes could help to predict if a patient has a high risk of 374 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020. . https://doi.org/10.1101/2020.05.05.20091355 doi: medRxiv preprint developing severe disease even when the disease is in its early course. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020. . https://doi.org/10.1101/2020.05.05.20091355 doi: medRxiv preprint 19 response in patients with COVID-19 may be related to disease progression and severe 397 outcomes (3, 23). We and others have demonstrated that plasma levels of IL-6 and 398 IP-10 (CXCL10) in COVID-19 patients are significantly increased (3, 12). In our 399 study, the mRNA expression of cytokines IL1B and IL6 in PBMCs is also increased. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 

The IFNs play an essential role in the antiviral response, but viruses can also escape with SARS-COV, MERS-COV as well as H7N9 infection (11, 43, 44) . We speculated 455 that IFN-γ might induce large amounts of cytokines secretion to aggravate the illness.

Therefore, the drugs that lower the expression of IFNγ may help to alleviate the 457 disease severity. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020. to COVID-19 patients that may help to recover from lymphopenia. Finally, we found 475 that the four S1PR receptors decreased to the lowest level in S2, which may reduce 476 the egression of T and NK cells from peripheral lymphoid organs to the blood (50, 477 51).

In contrast, S1PR3, which is mainly expressed by monocytes, had increased 480 expression in S2. The inhibitors of S1PR1 and S1PR5, Fingolimod (FTY720, trade (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020. . https://doi.org/10.1101/2020.05.05.20091355 doi: medRxiv preprint 23 Some COVID-19 patients developed disordered coagulation, including disseminated 485 intravascular coagulation (DIC) (7, 52, 53) . We found that COVID-19 patients had 486 increased expression of many coagulation associated genes started from early-onset 487 and lasted two-third of the disease course. The patients in this study also had 488 increased D-dimer in the serum, indicating fibrin deposition. As a consequence, 489 PLAU and PLAUR were elevated for fibrin degradation. However, it was reported (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020. We identified 25 gene expression markers that can be used to identify the disease 509 progression. Importantly, combinatorial analysis of expression of 12 genes may be 510 able to predict if a patient is at high risk of progressing to severe disease, even when (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 Peripheral blood mononuclear cells (PBMCs) from 4 female COVID-19 patients were 529 obtained from Guangzhou Eighth People's Hospital of Guangzhou Medical University.

All patients signed informed consent for this study. This study was approved by (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 26 forward, 5'-TGCCTCAGTATGCTGGCTCT, reverse, 5'-GAGACCTTTGCCTCCTT 551 GTTC. Primers for CD8A: forward, 5'-TCCTCCTATACCTCTCCCAAAAC, reverse, 552 5'-GGAAGACCGGCACGAAGTG. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 using the (QIA seq FastSlect-rRNA HRM KIT, QIAGEN) according to the 573 manufacturer's instructions, Ribosomal RNA depletion was confirmed by using 574 Agilent Bioanalyzer analysis and noting the absence of ribosomal peaks. Next, Prime, 575 and Fragment Mix from the (NEB Next RNA library prep Kit) were added to each Raw RNA-seq reads were filtered according to their base qualities, read sequences 586 were trimmed at 3'end after reaching a 2-base sliding window with PHRED quality 587 score lower than 20. Following filtering, Illumina adapter sequences at 3'end were 588 removed using Trimmomatic v0.36 (58). After low-quality filtering and adapter 589 trimming, reads more than 50 nt in length were retained for further analysis. Next, the 590 trimmed reads were mapped to the human (hg38) and SARS-COV-2 viral 591 (Wuhan-Hu-1) reference genomes (3) using HISAT v2.1 (59) with corresponding 592 gene annotations (Gencode GRCh37/V32 for the human genome) with default 593 settings RF, respectively. Total counts per mapped gene were determined using 594 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. 

Raw counts matrix obtained from featureCounts was used as input for differentially 598 expression gene analysis with the bioconductor package edgeR v3.28 (61) or DESeq2 599 v1.26 (62) in R v3.6. Gene counts more than 5 reads in a single sample or more than To determine the similarity between the PBMCs sample, the normalized counts matrix 616 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020. . https://doi.org/10. 1101 29 was used to calculated pairwise person coefficients using build-in function cor in R 617 software. Hierarchical clustering across all samples was based on pairwise Pearson 618 correlation coefficients matrix. And the build-in function prcomp in R software was 619 enlisted for principal component analysis. We inferred the immune cell quantities in 620 each blood sample using the CIBERSORT server (https://cibersortx.stanford.edu/). cluster number and all the DEGs were divided into 6 clusters according to their 630 expression pattern. subsequently, gene ontology analysis was performed to assess 631 their biological relevance using R package clusterprofiler v 3.14.3 (64). 634 We determine the expression of rearranged TCR and BCR in RNA-seq samples using 635 MiXCR v3.0.3 (65). MiXCR is a universal tool for fast and accurate analysis of T-636 and B-cell receptor repertoire sequencing data, which also provides parameter -p 637 rna-seq for processing RNA-seq data. Here, we aligned the filtered reads (see in 638 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 Pre-processing of the raw RNA-seq data) against reference V(D)J genes that 639 download from IMGT (http://www.imgt.org/). And the total matched counts of TCR 640 or BCR of each sample were normalized according to their sample size factors. To 641 analysis the diversity of TCR and BCR, we defined a clone as group of sequences 642 with the same VH gene, JH gene and identical CDR3 amino acid sequence. The 643 clonal diversity of TCR and BCR were analyzed using the R package Alakazam (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020. All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 

Clinical Characteristics of 138 Hospitalized Patients With

Clinical course and outcomes of critically ill patients with SARS-CoV-2 674 pneumonia in Wuhan, China: a single-centered, retrospective, observational study

Clinical features of patients infected with 2019 novel coronavirus in Wuhan

Clinical characteristics and imaging manifestations of the 2019 novel 681 coronavirus disease (COVID-19):A multi-center study in Wenzhou city

Clinical Characteristics of Children with Coronavirus Disease

Clinical Characteristics of Coronavirus Disease 2019 in China

Clinical characteristics of fatal and recovered cases of coronavirus disease 2019 688 (COVID-19) in Wuhan, China: a retrospective study

The cytokine storm of severe influenza and development of 690 immunomodulatory therapy

An interferon-gamma-related cytokine storm in SARS patients

Plasma inflammatory cytokines and chemokines in severe acute respiratory 694 syndrome

Pathogenic T cells and inflammatory monocytes incite inflammatory storm in 696 severe COVID-19 patients

Elevated exhaustion levels and reduced functional diversity of T cells in 698 peripheral blood may predict severe progression in COVID-19 patients

Functional exhaustion of antiviral lymphocytes in COVID-19 patients

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

An emerging role of neutrophils and 705

NETosis in chronic inflammation and fibrosis in systemic lupus erythematosus (SLE) and 706

ANCA-associated vasculitides (AAV): Implications for the pathogenesis and treatment

The role of neutrophils in the pathogenesis of systemic lupus 33

inhibitor leads to a pro-thrombotic potential in endothelial cells

HIF transcription factors, inflammation, 716 and immunity

pROC: an open-source package for R and S+ to analyze and compare ROC 718 curves

Clinical progression and viral load in a community outbreak of 720 coronavirus-associated SARS pneumonia: a prospective study

COVID-19: combining antiviral and anti-inflammatory treatments

Reducing mortality from 2019-nCoV: 724 host-directed therapies should be an option

IL-15 is required for sustained 726 lymphopenia-driven proliferation and accumulation of CD8 T cells

Naive T cell homeostasis: from awareness of space to a sense of 729 place

Normal T cell homeostasis: the conversion of naive cells into 731 memory-phenotype cells

COX-2-derived prostacyclin confers atheroprotection on female mice

Targeted prostaglandin E2 inhibition enhances antiviral immunity through 735 induction of type I interferon and apoptosis in macrophages

Prostaglandin E2 inhibits IFN-alpha secretion and Th1 costimulation by 737 human plasmacytoid dendritic cells via E-prostanoid 2 and E-prostanoid 4 receptor 738 engagement

Prostaglandin E2 740 suppresses lipopolysaccharide-stimulated IFN-beta production

PGE2 production at sites of tissue injury promotes an anti-inflammatory 743 neutrophil phenotype and determines the outcome of inflammation resolution in vivo

Coronavirus disease 19 (Covid-19) and 746 non-steroidal anti-inflammatory drugs (NSAID)

Neutrophil extracellular traps: double-edged swords of innate 748 immunity

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

Multi-platform 'Omics Analysis of Human Ebola Virus Disease Pathogenesis

Neutrophil Activation and Early Features of NET Formation Are 754

Associated With Dengue Virus Infection in Human

Dysregulation of immune response in patients with COVID-19 in Wuhan

SARS: prognosis, outcome and sequelae

A major outbreak of severe acute respiratory syndrome in Hong Kong

Expression profile of immune response genes in patients with Severe 761

Acute Respiratory Syndrome

Immune evasion strategies of flaviviruses

Early hypercytokinemia is associated with interferon-induced transmembrane 765 protein-3 dysfunction and predictive of fatal H7N9 infection

MERS-CoV 768 infection in humans is associated with a pro-inflammatory Th1 and Th17 cytokine profile

Hypothesis for potential pathogenesis of SARS-CoV-2 infection--a 771 review of immune changes in patients with viral pneumonia

Defining 'T cell exhaustion

Transcriptional regulation and T cell exhaustion

T-Cell Exhaustion in Chronic Infections: Reversing the State of Exhaustion

Reinvigorating Optimal Protective Immune Responses

L-Arginine modulates CD3zeta expression and T cell function in activated 779 human T lymphocytes

S1PR5 is essential for human natural killer cell migration toward 781 sphingosine-1 phosphate

Lymphocyte egress from thymus and peripheral lymphoid organs is 783 dependent on S1P receptor 1

Abnormal coagulation parameters are associated with poor 785 prognosis in patients with novel coronavirus pneumonia

Attention should be paid to venous thromboembolism prophylaxis in the 789 management of COVID-19

Incidence of thrombotic complications in critically ill ICU patients with 793 COVID-19

Mechanisms of severe acute respiratory syndrome coronavirus-induced 795 acute lung injury

Neutrophil extracellular trap (NET) impact on deep vein 797 thrombosis

Trimmomatic: a flexible trimmer for Illumina sequence 799 data

HISAT: a fast spliced aligner with low memory

All rights reserved. No reuse allowed without permission.(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 35 requirements. Nat Methods 12, 357-360 (2015 818 819 820 All rights reserved. No reuse allowed without permission.(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 37 stages (S1, S2, S3, and R) and healthy control (H). Principal components 1 (PC1) and 841 2 (PC2), which represent the greatest variation in gene expression, are shown. (C) 842 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 38 Paired comparison between S1 and R, S2 and R, S3 and R for each COVID-19 patient. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020. . https://doi.org/10.1101/2020.05.05.20091355 doi: medRxiv preprint 41 876 All rights reserved. No reuse allowed without permission.(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020. . https://doi.org/10.1101/2020.05.05.20091355 doi: medRxiv preprint 1 All rights reserved. No reuse allowed without permission.(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020. arcs (orange, red, and blue) are used to present S1, S2, and S3 for each patient, curved 10 bands connecting each pair of arcs indicates the overlaps between each pair of stages, 11 All rights reserved. No reuse allowed without permission.(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 the width of the curved bands represent the number of overlapped DEGs at the two (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101 (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020. foldchange (middle) and adjusted p values (right) for S1PR1-S1PR5. (C) Heatmaps 47 All rights reserved. No reuse allowed without permission.(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020. . https://doi.org/10. 1101 showing relative expression level for S1PR1-S1PR5 in individual patients in four (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.The copyright holder for this preprint this version posted May 8, 2020 . . https://doi.org/10.1101