key: cord-0978235-35w3m49i
authors: Gisby, J. S.; Buang, N. B.; Papadaki, A.; Clarke, C. L.; Malik, T. H.; Medjeral-Thomas, N.; Pinheiro, D.; Mortimer, P. M.; Lewis, S.; Sandhu, E.; McAdoo, S. P.; Prendecki, M. F.; Willicombe, M.; Pickering, M. C.; Botto, M.; Thomas, D. C.; Peters, J. E.
title: Multi-omics identify LRRC15 as a COVID-19 severity predictor and persistent pro-thrombotic signals in convalescence
date: 2022-05-01
journal: nan
DOI: 10.1101/2022.04.29.22274267
sha: cf1399ea49cfa3dd25c3f2be0d0c19ef11ca035f
doc_id: 978235
cord_uid: 35w3m49i

Patients with end-stage kidney disease (ESKD) are at high risk of severe COVID-19. Here, we performed longitudinal blood sampling of ESKD haemodialysis patients with COVID-19, collecting samples pre-infection, serially during infection, and after clinical recovery. Using plasma proteomics, and RNA-sequencing and flow cytometry of immune cells, we identified transcriptional and proteomic signatures of COVID-19 severity, and found distinct temporal molecular profiles in patients with severe disease. Supervised learning revealed that the plasma proteome is a superior indicator of clinical severity than the PBMC transcriptome. Notably, we show that both the levels and trajectory of plasma LRRC15, a proposed co-receptor for SARS-CoV-2, are the strongest predictors of clinical outcome. Strikingly, we observed that two months after the acute infection, patients still display dysregulated gene expression related to vascular, platelet and coagulation pathways, including PF4, providing a potential explanation for the prolonged elevation of thrombotic risk following COVID-19.

induction of interferon-alpha/beta", " Wilk et al., 2021 IFN module" [26] , "Host-pathogen 166 interaction of human coronaviruses interferon induction" and "SARS-CoV-2 innate immunity 167 evasion and cell-specific immune response", reflecting host anti-viral responses and providing 168 validation of our analysis ( Figure 2C, Supplementary File 1E) . Highly up-regulated proteins 169 within these pathways included STAT1; DDX58 and ISG15, both crucial to the IFN-mediated 170 antiviral response in ; IFITM3, which is up-regulated in lung epithelial cells 171 during early SARS-CoV-2 infection [28] ; and the chemokines CXCL11, CXCL1, CXCL6, 172 CXCL5 and CXCL10. Another significantly up-regulated pathway was "Senescence-173 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 associated secretory phenotype", which included up-regulated ubiquitin-conjugating enzymes 174 (UBE2S, UBE2E1), histones (H2BC21, H2BU1) and STAT3 ( Figure 2D ). Down-regulated pathways included "Integrin cell surface interactions" and "Collagen biosynthesis and 176 modifying enzymes" which contained collagen proteins (e.g. COL11A2, COL13A1, COL15A1) 177 and related enzymes (e.g. P4HB, PCOLCE) ( Figure 2D) . 178 179 Transcriptomic and proteomic changes associated with COVID-19 severity 180 181 In both cohorts, the PCA of the PBMC transcriptomics revealed differences according to both 182 severity at time of sampling and overall clinical course (defined by peak severity score) ( Figure  183   3A ). There was a gradient of severity reflected in the molecular phenotype. We next assessed 184 molecular features associated with severity at time of blood sampling, encoded as an ordinal 185 variable. We identified 3,522 genes that were significantly associated with contemporaneous 186 severity in the Wave 1 cohort and 657 genes in the Wave 2 cohort (LMM, 1% FDR, 187 Supplementary File 1F, Supplementary Figure 6 ). We then applied GSVA to identify 188 pathways and used RRA to combine results from each cohort (Supplementary File 1G) . 189

The up-regulated transcriptomic pathways in more severe disease included those involved in 191 oxidative stress ("Glutathione metabolism", "Detoxification of reactive oxygen species"), 192 "Transcriptional regulation of granulopoiesis", pathways containing numerous histone-193 encoding genes ("HDACs deacetylate histones", "Diseases of programmed cell death", " RHO 194 GTPases activate PKNs") and "Complement and coagulation cascades" (Figure 3B -C, 195

Supplementary File 1G). Down-regulated pathway terms included "TCRA pathway ", 196 "Pathogenesis of SARS-CoV-2 mediated by nsp9-nsp10 complex", "TP53 activity", and "PD1 197 signaling", suggesting T cell activation in more severe COVID-19 ( Figure 3B pathways in more severe disease were "HDACs deacetylate histones", pathways related to 207 transcriptional regulation (e.g. "mRNA splicing minor pathway", "Spliceosome", "RNA 208 polymerase II transcription termination", "Processing of capped intron-containing pre mRNA") 209 and "RUNX1 regulates genes involved in megakaryocyte differentiation and platelet function", 210 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 while the most down-regulated pathways included "PD-1 signaling" and "T-cell receptor and 211 costimulatory signaling". Example proteins from these pathways are shown in Supplementary 212 Figure 7C . 213 214 Severe COVID-19 is associated with dynamic multi-omic modular trajectories 215 216 We next examined the temporal trajectories of the transcriptome and the proteome during 217 COVID-19 by explicitly modelling molecular profiles with respect to time following symptom 218 onset (Methods). To aid biological interpretation, we first applied a dimension reduction 219 strategy using weighted gene correlation network analysis (WGCNA) [29] . WGCNA identified 220 23 modules of co-expressed genes (which we denote with the prefix 't') ( Supplementary File 1N). The severity-associated modules tB and tJ were both strongly 245 positively associated with myeloid cell proportions, particularly neutrophils, and negatively 246 associated with lymphocyte subsets (Supplementary Figure 11) . The presence of a 247 . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted May 1, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 neutrophilic gene signature in the PBMC preparations may indicate the presence of low-248 density granulocytes. Consistent with this, hub genes in Module tB (including TECPR2, 249 CSF3R, STX3; Figure 4A ) are associated with granulocytes and autophagy, and pathway 250 analysis of the module genes revealed enrichment for pathways including "Neutrophil 251 degranulation" and "ROS and RNS production in phagocytes" (including genes encoding the 252 key cytosolic components of the phagocyte NADPH oxidase such as NCF1, NCF2 and NCF4). 253

Module tB also contains genes encoding calcium-binding proteins (e.g. S100A6, S100A9, 254 S100A11, S100A12) that play important roles in regulating inflammatory pathways [31], as 255 well as integrins (e.g. ITGA1, ITGAM, ITGB4, ITGAX, ITGAD), adhesion molecules (e.g. 256 CEACAM1, CEACAM3, CEACAM4, ICAM3), OSM (encoding Oncostatin M) and CSF1 257 (encoding M-CSF). The tL module, which also displayed a rising trajectory in worse disease, 258 was strongly positively associated with imputed plasma cell proportion (Supplementary 259 Figure 11 ) and many of its members encoded immunoglobulins. The severity-associated 260 proteomic modules that strongly correlated with transcriptomic modules tB, tJ and tL were p8 261 and p9 (both enriched for pathways related to RNA splicing), and p12 (significantly enriched 262 for the pathway "HDACs deacetylate histones") (Supplementary Table 4 ). The latter is 263 consistent with our earlier observations that a histone pathway signature was prominently 264 associated with COVID-19 severity in both the RNA-seq ( Figure proportions, we performed flow cytometry on a subset of PBMC samples from the Wave 2 281 cohort. We found no major difference in the overall proportions of myeloid or lymphoid cells 282 within the PBMC fraction between pre-infection and COVID-19 positive samples, except for a 283 reduction in the proportion of type 2 dendritic cells (Supplementary Figure 12) . Similarly, 284 there was little difference in the distribution of cells between mild/moderate and severe/critical 285 patients. We observed some severity-related differences within cell subsets. Within lymphoid 286 cells, we noted higher expression of the activation marker CD69 on CD4+ T cells at day 7 in 287 severe/critical disease compared to either pre-infection or mild/moderate disease 288 (Supplementary Figure 13A) . At day 14, there was an increase in CD38 hi plasmablasts in 289 severe/critical disease compared to pre-infection or mild/moderate samples (Supplementary 290 Figure 13B ). We also found that in severe/critical patients, there was a progressive drop in 291 the proportion of non-classical monocytes over the first 14 days of the illness that was more 292 marked than in mild/moderate patients (Supplementary Figure 14A) Figure 14B) , likely reflecting enhanced activation [30] . In classical 296 monocytes there was a similar, but non-significant, trend. We found higher expression of 297 proliferation-associated Ki67 on classical monocytes in COVID-19 versus pre-infection 298 samples in both mild/moderate and severe/critical patients (Supplementary Figure 14C) . In 299 our transcriptomic data we identified increased SIGLEC1 gene expression in COVID-19 300 ( Figure 2B ). SIGLEC-1 is exclusively expressed by CD14+ monocytes at the protein level. 301 SIGLEC-1 expression measured by flow cytometry correlated with GSVA enrichment score of 302 type I IFN signatures (Supplementary Figure 14D) . We observed SIGLEC1 expression 303 increased at greater intensity as early as day 0-3 post infection in severe/critical versus 304 mild/moderate patients, suggesting stronger and a more immediate type I IFN response in 305 severe COVID-19 (Supplementary Figure 14E) . 306 307 Longitudinal cytokine/chemokine analysis reveals distinct temporal profiles that distinguish 308 disease severity 309 310 Many plasma proteins associated with severe COVID-19 are canonically intra-cellular 311 proteins. Their elevation in severe COVID-19 may therefore be a readout of increased cell 312 turnover, death, stress, and viral hijacking of host cellular machinery. Consequently, we 313 performed a more focussed analysis examining proteins whose primary biological role is to 314 act extra-cellularly (e.g. cytokines, chemokines, growth factors and their receptors). These 315 classes of proteins are important therapeutic targets in inflammatory diseases [32] . 316 Accordingly, we modelled the temporal profiles of 232 proteins that fell within the KEGG 317 pathway "Cytokine-cytokine receptor interaction". Fifty proteins had significantly different 318 profiles in patients with a severe/critical clinical course versus those with mild/moderate ones 319 (TxCC interaction effect, 5% FDR; Supplementary File 1O). Proteins exhibited distinct 320 patterns of divergence between severe/critical and mild/moderate disease over time (Figure 321 . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted May 1, 2022. ; https://doi.org/10.1101 doi: medRxiv preprint 5A). Some (e.g. IL1b, IL6, IL15RA, CCL2) showed a relatively stable temporal profile in 322 mild/moderate patients but rising trajectories in severe/critical patients ( Figure 5B ). Others 323 (e.g. CCL15, TNFSF13B (BAFF), PDGFRB, EDAR, IFNA10, IFNA13, IFNA16, IFNE and  324 IFNL3) were elevated early in the disease course and decreased over time, but displayed 325 more marked initial elevations in severe/critical patients ( Figure 5C ). Yet other proteins 326 displayed temporal profiles in mild/moderate patients that were inverted compared to 327 severe/critical. For example, CD40LG, TNFSF10 (TRAIL) and IL11 were reduced in the 328 severe/critical versus the mild/moderate group at early timepoints but increased in 329 severe/critical patients later ( Figure 5D) . Conversely, leptin, INHBA (inhibin A), and CCL22 330 were initially higher in severe/critical than mild/moderate patients but with the reverse pattern 331 later on ( Figure 5E ). These data illustrate the dynamic nature of the soluble protein response 332 and how this varies according to disease severity, highlighting the limitation of studies that use 333 a single snapshot. 334 335 Plasma LRRC15 as a predictor of COVID-19 severity 336 337 We next investigated whether clinical severity could be inferred from the transcriptomic and/or 338 proteomic data and which had the better predictive performance. For this analysis, we 339 combined the COVID-19 cases from both cohorts. For each COVID-19 patient, we selected 340 the first sample at the patient's peak severity score so that there was one sample per patient. 341

To predict COVID-19 severity at time of sampling, we employed two supervised learning 342 methods, lasso and random forests. We applied these separately on i) the plasma proteomic 343 data; ii) the PBMC transcriptomic data; and iii) the combination of both (the multi-omic data). CC-BY-ND 4.0 International license It is made available under a 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 1, 2022.

important predictors (10/15 for lasso and 9/15 for random forests). This, and our finding that 359 the plasma proteome was a superior predictor of severity than the PBMC transcriptome, 360 highlights that plasma proteins provide a valuable read-out of the pathophysiological 361 processes in severe 363 Importantly, both lasso and random forests identified plasma LRRC15 protein levels as the 364 most important predictor of COVID-19 severity. Interestingly, this protein was recently 365 identified by two pre-prints as a receptor for SARS-COV-2 [33, 34] . We next examined 366 Up-regulated clotting-related genes included PF4 (encoding platelet factor 4) and the related 384 gene PF4V1 (platelet factor 4 variant 1). Of note, these genes are located in the same genomic 385 region on chromosome 4, along with the chemokine CXCL5, which was also significantly up-386 regulated. Another nearby gene, PPBP (Pro-Platelet Basic Protein, aka CXCL7), was also up-387 regulated in convalescent samples, although it did not quite reach significance at 1% FDR 388 (nominal P = 3.33x10 -5 , adjusted P = 0.0159). The upregulation of these neighbouring genes 389 suggests they are influenced by a shared genomic regulatory element. Overrepresentation 390 analysis of the 25 differentially expressed genes revealed significant enrichment of terms 391 including "Platelet activation, signaling and aggregation", "Formation of fibrin clot/clotting 392 cascade", "Chemokine signaling pathway", "SARS-CoV-2 innate immunity evasion and cell-393 specific immune response" and "Smooth muscle contraction" (Figure 7C, Supplementary  394 File 1T). These data suggest persistent activation of abnormal processes for a considerable 395 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. genes reflecting leukocyte-vascular interactions. Other up-regulated pathways included "Polo-428 like kinase mediated events" and "Golgi-cisternae peri-centriolar stack re-organisation". Both 429 are likely to reflect the extensive cell division of immunocytes that occurs in COVID-19. For 430 instance, the pericentriolar stacks of Golgi cisternae undergo extensive fragmentation and 431 reorganization in mitosis. Similarly, polo-like kinase is crucial for facilitating the G2/M 432 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10. 1101 transition. These findings are consistent with the up-regulation of APC-Cdc20 mediated 433 degradation of Nek2A and other APC-Cdc20 related processes that we observed in the 434 proteomic data; Cdc20 is a protein that is key to the process of cell division. cytokines/chemokines revealed divergence temporal trajectories between disease severity 469 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101/2022.04.29.22274267 doi: medRxiv preprint strata, manifesting in several patterns ( Figure 5 ). For example, in patients with a 470 severe/critical disease course, IL11 was reduced early on but increased later relative to more 471 indolent disease ( Figure 5D ). IL11 is known to cause progressive fibrosis [44, 45] , and the 472 marked increases late in severe/critical disease may have implications for the development of 473 pulmonary sequelae. Leptin, INHBA (inhibin A), and CCL22 showed the opposite pattern 474 ( Figure 5E ). Leptin has roles in both cell metabolism and immunity with many immune cells 475 responding to leptin directly via the leptin receptor, resulting in a pro-inflammatory phenotype 476

[46]. It is produced by adipocytes, so its elevation early in severe/critical disease may be a 477 read-out of higher body mass index, which is a risk factor for severe or 50] , our data suggest that the picture may be more complex. 492

Thus, IFNs may act as a double-edged sword, with harm to the host from both insufficient 493 responses (leading to failure to control the virus) and from excessive responses (resulting in 494 immunopathology). While we cannot exclude the possibility that increased IFNs is a 495 consequence rather than a cause of severe disease, their elevation very early in disease 496 suggests this is less likely. Another consideration is that the greater IFN response in severe 497 disease might reflect higher viral burden. 498 499 Using two distinct supervised learning methods, we observed that the plasma proteome better 500 captures disease severity than the PBMC transcriptome. When supervised learning algorithms 501 were trained on both the proteomic and transcriptomic data simultaneously, plasma proteins 502 dominate the list of important biomarkers. There are several reasons why this might be the 503 case. Plasma is under strong homeostasis: derangement is a marker of loss of physiological 504 control. Plasma proteins may provide important read-outs of both pathogenesis and tissue 505 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101/2022.04.29.22274267 doi: medRxiv preprint injury by reflecting the activity of cell types other than PBMCs, such as neutrophils, 506 endothelium and hepatocytes (a major source of coagulation and complement proteins). 507 508 A striking finding was the predictive value of plasma levels of LRRC15 in indicating COVID-509 19 severity. Longitudinal profiling revealed that LRRC15 levels remain stable in those with a 510 mild/moderate clinical course but decrease over time in severe/critical illness ( Figure 6C ). Our 511 findings are particularly intriguing as two recent pre-prints have identified LRRC15 as an 512 accessory factor for SARS-CoV-2 entry to cells. Using arrayed transmembrane protein and 513 pooled genome-wide CRISPR activation screens, Shilts and colleagues demonstrated that 514 the SARS-CoV-2 spike protein interacts with LRRC15 [33] . Both screens identified the 515 interaction and the CRISPRa screen identified LRRC15 and the established SARS-CoV-2 516 binding partner, ACE2, as the two most prominent interactors. This work also showed that 517 ACE2 and LRRC15 bind the C-terminal domain of the spike protein, which contains the 518 receptor binding domain, suggesting that the two proteins may compete for spike protein 519 binding. Song and colleagues also used a CRISPRa approach to identify proteins that could 520 bind the SARS-CoV-2 spike protein to the A375 melanoma cell line [34] . The screen identified 521 ACE2 and LRRC15, and further showed that the interaction took place with the receptor 522 binding domain of the spike protein. Expression of LRRC15 on a HeLa cell line that expresses 523 ACE2 inhibited the entry of a SARS-CoV-2 spike pseudovirus. This paper notes, however, 524 that LRRC15 is expressed on different cells from those that express ACE2 and proposes that 525 LRRC15 inhibits virally entry in trans, acting as a decoy and binding virions that cannot then 526 enter cells via ACE2. Our data are consistent with a model in which a failure to up-regulate 527 LRRC15 increases risk of severe COVID-19 disease because of the lack of a receptor that 528 inhibits its entry to cells. Thus our study is the first human in vivo study to highlight the 529 importance of LRRC15 in the response to SARS-CoV-2. 530 531 Another unique strength of our study was the availability of baseline pre-infection samples for 532 the Wave 2 cohort, as well as samples taken two months after the acute COVID-19 episode. 533

Leveraging this, we demonstrate that there is chronic activation of vascular, platelet and 534 coagulation pathways for a prolonged period after clinical resolution of disease. The elevated 535 risk of thrombotic events during acute COVID-19 is well-documented. In a large study 536 encompassing both hospitalised and non-hospitalised patients [51], the risk of pulmonary 537 embolism (PE) and deep vein thrombosis (DVT) were 27-fold and 17-fold increased, 538 respectively, in the seven days following diagnosis. These risk ratios are much higher than 539 those previously associated with upper respiratory tract infections, suggesting unique features 540 specific to SARS-CoV-2 infection. The risk of arterial thrombosis was also significantly 541 increased, although smaller in magnitude than the risk of venous thromboembolism (VTE). 542 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10. 1101 The pathophysiology underlying COVID-19 associated coagulopathy is complex and involves 543 the convergence of several pathways [52] . Invasion of ACE2-expressing epithelial cells by 544 SARS-CoV-2 results in down-regulation of ACE2 and increased angiotensin II levels. This in 545 turn leads to increased expression of PAI1 which impairs breakdown of fibrin and promotes 546 increased vascular tone, via smooth muscle contraction. Endothelial cell activation, 547 complement activation, NETosis, hypoxia and cytokine/chemokine secretion all promote 548 coagulopathy through increases in tissue factor and concomitant fibrin formation. Remarkably, 549 our data suggest that these pathways remain dysregulated months after acute infection has 550 resolved (Figure 7, Table 1 ). This is particularly important given emerging evidence indicating 551 that the risk of thrombo-embolism extends beyond the acute phase. Ho et al showed that risk 552 of a PE was 3.5-fold higher even in the time window 28 to 56 days after diagnosis of . A recent population-wide registry study revealed that following COVID-19 the risk of 554 DVT and PE was significantly elevated for 70 and 110 days, respectively [53] . Although VTE 555 risk was greatest for those with severe disease, even patients with mild disease had elevated 556 VTE risk. Our data provide a molecular basis that begins to explain this risk. Intriguingly, 557 among the genes up-regulated in convalescent samples compared to pre-infection was 558 platelet factor 4 (PF4). PF4 is expressed in platelets and leucocytes. It is released from the 559 alpha granules of activated platelets, contributing to platelet aggregation. The prolonged up-560 regulation of PF4 after COVID-19 is therefore likely to contribute to a prothrombotic state. Of 561 note, autoantibodies to PF4 are the pathogenic entity in both vaccine-induced thrombotic 562 thrombocytopenia (VITT) [54,55] and heparin-induced thrombocytopenia (HIT). PF4 becomes 563 an autoantigen when it forms complexes with adenoviral vaccine components or heparin 564 respectively, unmasking epitopes to which autoantibodies bind [56] . It will therefore be 565 interesting for future studies to investigate whether autoantibodies to PF4 might contribute to 566 post-COVID-19 thrombosis in some patients. Whether the molecular abnormalities found in 567 our study also apply to more general patient populations without background ESKD needs to 568 be determined. Ongoing studies focusing on the sequelae of COVID-19 are well placed to 569 address this. 570 571 Our study has several limitations. ESKD patients have considerable multi-morbidity and 572 deranged physiology, and our findings may not all be generalisable to other patient 573 populations. We lacked a comparator group of ESKD patients with another viral infection to 574 delineate COVID-19 specific features. We studied peripheral blood; while this can provide 575 valuable information, it does not always reflect processes at the site of tissue injury. We 576 performed bulk RNA-seq on PBMCs. Thus, transcriptomic signatures may reflect both 577 changes in gene expression and also alteration in the distribution of cell subtypes within 578

PBMCs. We mitigated this issue through use of deconvolution methods and flow cytometry, 579 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101/2022.04.29.22274267 doi: medRxiv preprint but future studies using single cell RNA-seq and CITE-seq will provide further granularity. We 580 did not have measurements of viral load which would have aided interpretation of the 581 magnitude of host responses (e.g. interferon signaling). Finally, the convalescent samples 582 were taken relatively soon after clinical recovery: it will be important for future studies to 583 establish how long molecular abnormalities persist. 584

In summary, we demonstrate dynamic transcriptomic, proteomic and cellular signatures that 586 vary both with time and COVID-19 severity. We show that in patients with a severe clinical 587 course there is increased type 1 interferon signaling early in the illness, with increases in pro-588 inflammatory cytokines later in disease. We identify plasma levels of the proposed alternative 589 SARS-CoV-2 receptor, LRRC15, as the strongest predictor of COVID-19 severity. Finally, we 590

show that immune cells display dysregulated gene expression two months following COVID-591 We recruited two cohorts of ESKD patients with COVID-19 ( Figure 1A ). All patients were 605

receiving haemodialysis prior to acquiring COVID-19. The first cohort ('Wave 1') were recruited 606 during the initial phase of the COVID-19 pandemic (April-May 2020). Blood samples were 607 taken from 53 patients with COVID-19 (Supplementary Table 1) . Serial blood sampling was 608 carried out where feasible ( Figure 1B) , given the pressure on hospital services and the effects 609 of national lockdown. We also contemporaneously recruited 59 non-infected haemodialysis 610 patients to provide a control group, selected to mirror the age, sex and ethnicity distribution of 611 the COVID-19 cases (Supplementary Figure 1A-C) . 612 613

The Wave 2 cohort consisted of 17 ESKD patients with COVID-19 infected during the 614 resurgence of cases in January-March 2021 (Supplementary Table 2 ). These 17 individuals 615 had all been recruited as part of the COVID-19 negative control group during Wave 1, and so 616 a pre-infection sample collected in April/May 2020 (8-9 months preceding infection) was also 617 available. For the Wave 2 cohort, we systematically acquired serial samples for all patients at 618 regular intervals (every 2-3 days over the course of the acute illness) ( Figure 1C) . Additionally, 619

for 12 of these 17 patients, we acquired convalescent samples at approximately 2 months 620 post the acute COVID-19 episode (range 41-55 days from the initial sample). Convalescent 621 samples were unavailable for four patients who died and for one patient due to logistical 622 difficulties in sample collection. 623 624

To minimise variation related to the timing of dialysis, blood samples were taken prior to 625 commencing a haemodialysis session. 626 627 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10. 1101 For the Wave 1 cohort, following quality control (QC), transcriptomic data were available for 683 179 samples from 51 COVID-19 positive ESKD patients (median 3 samples per patient, range 684 1-8) (Supplementary Figure 1D) proteins and 1,160 as secreted (Supplementary Figure 16A) . Many proteins were labelled 707 as both intracellular and as membrane or secreted, reflecting the biology of protein storage 708 and extra-cellular secretion/excretion (Supplementary Figure 16B) . 709 710

We report proteins by their corresponding HGNC gene ID [66] , which provides a more 711

standardised nomenclature compared to protein names and allows direct comparison with the 712 transcriptomic data. 713 714

Where multiple Somamers related to the same protein, we retained these Somamers for 715 univariate analyses such as differential abundance analyses. However, for analyses that 716 considered multiple proteins simultaneously (PCA, WGCNA, MEFISTO, supervised learning), 717

we selected one Somamer at random to represent each protein. One COVID-19 positive 718 sample in the wave 2 cohort failed QC and was excluded from the analyses. The raw 719

SomaScan data was separated by cohort. individuals in the study and thus account for repeated measures. We performed differential 733 expression analyses for the transcriptomic data and the proteomic data. The regression model 734 for these analyses in Wilkinson-style notation was: 735 736 E ~ covid_status + sex + age + ethnicity + (1 | individual) 737 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. When reporting the number of differentially expressed proteins in the text we refer to the 756 number of unique proteins rather than the number of significant Somamers. 757 758

Testing transcriptomic and proteomic features for association with COVID-19 severity 759 We performed a within-cases analysis, testing for the association of gene expression with 760 COVID-19 severity at time of sampling. We used the four-level WHO severity rating (mild, 761 moderate, severe, critical), which could vary between samples from the same individual 762 reflecting the clinical status at the time the same was taken. We again used a linear mixed 763 model to account for samples from the same individual. The regression model was: 764 765 E ~ covid_severity + sex + age + ethnicity + (1 | individual) 766 767

The "covid_severity" variable represents severity at the time of the sample and was encoded 768 using orthogonal polynomial contrasts to account for ordinal nature of severity levels. 769 770 COVID-19 positive samples from the Wave 1 cohort were analysed separately to those from 771 the Wave 2 cohort. 772 773

The same approach was used for the proteomics data. 774 775

Gene set variation analysis 776

To identify pathways that were up-or down-regulated in COVID-19 positive versus negative 777 samples, we applied gene set variation analysis ( To dissect out the key molecules underpinning enriched pathways, we examined the genes 787 that comprise these pathway terms and identified which of these featured most prominently in 788 the differential gene expression analysis. 789 790

We repeated this procedure for testing of association of pathways with severity at the time of 791 sample using the 4-level ordinal score. 792 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101/2022.04.29.22274267 doi: medRxiv preprint

We then applied the same approach to the proteomics data for the COVID-19 positive versus 794 negative analysis, and for testing associations with COVID-19 severity at the time of sample. 795 796

Robust rank aggregation 797

The Wave 1 and Wave 2 cohorts were analysed separately for both the differential expression 798 analyses between COVID-19 positive and negative samples and for the within-cases severity 799

analyses. To identify the associations that were most consistent between the Wave 1 and 800

Wave 2 cohorts, for each analysis, we integrated the P-values for each cohort using robust 801 rank aggregation (RRA) [75] . This method identifies features that are ranked higher than 802 expected across multiple lists. RRA generates a significance score analogous to a P-value; 803

we -log10 transform these values such that a larger score indicates more consistent 804 associations between the Wave 1 cohort and the Wave 2 cohort. RRA was applied to the 805 results of the transcriptomic, proteomic and GSVA analyses comparing COVID-19 positive 806 versus negative samples from Wave 1 and Wave 2. Similarly, it was applied to the analyses 807 testing for association of molecular features with COVID-19 severity at the time of sampling. 808 809

Modelling modular longitudinal trajectories 810

We examined the temporal trajectories of the transcriptome following infection, by explicitly 811 modelling molecular markers with respect to time following COVID-19 symptom onset. We 812 used a two-step approach. 813 814

Step 1. To aid biological interpretation, we first applied a dimension reduction strategy using 815 weighted gene correlation network analysis (WGCNA) [29] to identify modules of correlated 816 molecular features. For this analysis, we combined samples from the Wave 1 and Wave 2 817 cohorts. Additionally, since our goal was to perform longitudinal analysis, we only selected 818 patients who had been sampled at least three times prior to 21 days following COVID-19 819 symptom onset. The default implementation of WGCNA is not designed for use with non-820 independent samples [76], so we modified the analysis pipeline by generating a correlation 821 matrix using a repeated measures correlation metric (rmcorr) that is appropriate for repeated 822 measures [77] . We used WGCNA's pickSoftThreshold.fromSimilarity function to pick the 823 minimum soft-thresholding power that satisfied the minimum scale free topology fitting index 824 (R 2 >0.85) and maximum mean connectivity (100). We subsequently defined signed adjacency 825 and topological overlap matrices before applying average-linkage hierarchical clustering. We 826 cut this tree with a hybrid dynamic tree cutting algorithm, with the parameters deepSplit = 4 827 and minClusterSize = 30 [78] . Finally, we defined eigengenes for each module and merged 828 those with a distance less than 0.25. The eigen-genes provide a numerical representation for 829 each module of co-expressed genes. 830 831

We used the same approach to analyse the proteomic data. 832 833

Step 2 The regression model used is displayed using Wilkinson-style notation below. 847 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101/2022.04.29.22274267 doi: medRxiv preprint eigen-expression ~ clinical_course * time + sex + age + ethnicity + wave + (1 | individual) 849 850 We extracted the P-values for the TxCC term in this model and applied Benjamini-Hochberg 851 adjustment, using 5% FDR as the significance threshold. A significant interaction effect for the 852 TxCC term indicates that the module has a different temporal profile in mild/moderate versus 853 severe/critical disease. 854 855

Additional WGCNA module annotation and association testing 856

To better understand the biological information reflected in the transcriptomic and proteomic 857 modules, we further characterised them through a multi-pronged analytical strategy. We 858 tested association of eigen-genes and eigen-proteins with other variables. First, we tested for 859 the association of the modules with WHO severity at the time of the sample using the LMM 860 approach described above in subsection 'Testing transcriptomic and proteomic features for 861 association with COVID-19 severity' i.e.: 862 863 E ~ covid_severity + sex + age + ethnicity + wave + (1 | individual) 864 865

Second, since PBMCs represent a mixed population of immune cells, we investigated whether 866 disease trajectory-associated transcriptomic modules might reflect shift in cell type 867

proportions. To this end, we applied CIBERSORTx, a computational algorithm to impute 868 immune cell fractions from RNA-seq data (see subsection 'Cell fraction imputation' below). 869 We then tested for correlations between these imputed immune cell proportions and module 870 eigengenes using LMM: 871 872 eigen-expression ~ cell_fraction + sex + age + ethnicity + wave + (1 | individual) 873 874

Both these models included an additional fixed effect ("wave") to reflect the cohort. 875 876

Third, we performed pathway enrichment analysis on the modules using the R package 877 clusterProfiler's "enricher" We used this method to find joint factors of variation in the transcriptomic and proteomic 895 datasets. For the MEFISTO analysis, we used the same set of samples as in the network 896 analysis and applied the same pre-processing steps to the data (see Methods -network 897 analysis). Additionally, we removed genes with the lowest maximum absolute deviation [67] 898 such that the number of genes retained were equal to the number of unique proteins measured 899 (6,323) to avoid imbalance numbers of features between the transcriptomic and proteomic 900 data which can impact the MEFISTO algorithm. Using the "slow" convergence criterion, 901 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 MEFISTO identified 8 factors that had a minimal variance explained of 1% in at least one data 902 modality. 903 904 We then applied the longitudinal model described earlier to test for an interaction effect 905 between time from first symptoms and clinical course, with a latent factor identified by 906 MEFISTO as the dependent variable. The regression model used is displayed using 907

Wilkinson-style notation below: 908 909

Latent_factor ~ clinical_course * time + sex + age + ethnicity + wave + (1 | individual) 910 911

Longitudinal modelling of cytokines and cytokine receptors 912

We modelled the temporal profiles of 232 plasma proteins that fell within the KEGG pathway 913

"Cytokine-cytokine receptor interaction". The goal of this analysis was to predict clinical severity from the molecular features 923 (transcriptomic, proteomic or both). We performed supervised learning using the R caret 924 framework [84] ; caret uses the randomForest package to fit random forest models and glmnet 925

[85] to fit lasso models. For this analysis, we only included samples on which both 926 transcriptomics and proteomics had been performed. We then selected the earliest sample for 927 each individual at which they had reached their peak COVID-19 WHO severity score, so that 928 there was one sample per patient. We then categorised the clinical severity score 929 corresponding to each sample into a binary variable such that patients with a WHO severity 930 score of mild or moderate were considered "mild/moderate" and those with a WHO score of 931 severe or critical were considered "severe/critical". This resulted in n=37 mild/moderate 932 samples and n=14 severe/critical samples. 933 934

We trained models using Monte Carlo cross-validation for: i) the plasma proteomic data alone 935 (6,323 features); ii) the PBMC RNA-seq data alone (12,225 features); and iii) the combined 936 proteomic and RNA-seq datasets. The first step in this training process was to create 200 937 random partitions of the data, such that 80% of the data was used to train the model in each 938 resample and 20% was retained as a validation set. In each resample, we calculated the area 939 under the curve (AUC) of the receiver operating characteristic (ROC) curve. We then 940 calculated confidence intervals for the 200 AUC-ROC values generated for each model and 941 feature type. 942 943

The random forest model's parameters were kept constant at 500 trees and the mtry value 944 (number of proteins randomly sampled as candidates at each node) was calculated as the 945 square root of the number of features. After cross-validation, we fitted a final random forest 946 model using the entirety of the dataset. We extracted important features from this model using 947 the R randomForestExplainer package, based on the accuracy decrease metric (the average 948 decrease in prediction accuracy upon swapping out a feature). For the lasso model, the 949 lambda value that maximised the mean AUC-ROC during cross-validation was selected. We 950 recorded the features selected by the lasso model in each data resample; feature importance 951 was subsequently defined as the number of models in which each feature had a non-zero 952 coefficient. The feature importance metrics from both models were scaled by dividing their 953 values by the maximum value, such that the most important feature has an importance metric 954 of 1. 955 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. Data and code availability 1008 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101/2022.04.29.22274267 doi: medRxiv preprint . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. 

. CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. 

. CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. 

. CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. 

. CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101/2022.04.29.22274267 doi: medRxiv preprint 1100 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101/2022.04.29.22274267 doi: medRxiv preprint CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101/2022.04.29.22274267 doi: medRxiv preprint 

. CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101/2022.04.29.22274267 doi: medRxiv preprint 1119 1120 

. CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 1 and Wave 2 infected vs non-infected differential expression analyses for the transcriptome. 1151

Each point is a gene, coloured according to its significance in the Wave 1 and 2 analyses. r= 1152

Pearson's correlation coefficient. B) As A, for the proteome. 1153

. CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 cases and controls, according to robust rank aggregation (RRA). These genes were significant 1157 in both cohorts (1% FDR). Columns are ordered by COVID-19 status and severity. For the 1158

Wave 1 heatmap, genes are ordered by hierarchical clustering; the Wave 2 heatmap is 1159 ordered to match this. Each feature was scaled and centred separately in each dataset. 1160 1161 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 heatmap, proteins are ordered by hierarchical clustering; the Wave 2 heatmap is ordered to 1177 match this. Each feature was scaled and centred separately in each dataset. 1178 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 contemporaneous severity, according to robust rank aggregation (RRA). These genes were 1183 significant in both cohorts (1% FDR). Columns are ordered by COVID-19 status and severity. 1184

For the Wave 1 heatmap, genes are ordered by hierarchical clustering; the Wave 2 heatmap 1185 is ordered to match this. Each feature was scaled and centred separately in each dataset. 1186 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 heatmap, proteins are ordered by hierarchical clustering; the Wave 2 heatmap is ordered to 1206 match this. Each feature was scaled and centred separately in each dataset. 1207 1208 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. represent 95% confidence intervals. 1215 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. effects. The bottom row shows the association with severity score at time of sample. Asterisks 1230 represent significant (5% FDR) associations. 1231 1232 . CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. CC-BY-ND 4.0 International license It is made available under a 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 1, 2022. COVID-19 severity with 95% confidence intervals using random forests. "Both" = supervised 1298 learning on the combined proteomic and transcriptomic data. B) Important proteins (left) and 1299 genes (right) for the random forests model. Feature importance is scaled between 0 and 1, 1300

where 1 represents the most important feature. 1301 1302 . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted May 1, 2022. annotations. 1307 1308 . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted May 1, 2022. . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted May 1, 2022. ; https://doi.org/10.1101/2022.04.29.22274267 doi: medRxiv preprint 6 (66.7%) 2 (22.2%) 0 (0.0%) 1 (11.1%) 5 (62.5%) 2 (25.0%) 1 (12.5%) 0 (0.0%) Diabetes 11 (64.7%) 6 (66.7%) 5 (62.5%) Current smoker 0 (0.0%) 0 (0.0%) 0 (0.0%) ESKD cause DN HTN/vascular GN/autoimmune Genetic Other/unknown 10 (58.8%) 0 (0.0%) 1 (5.9%) 0 (0.0%) 6 (35.3%) 6 (66.7%) 0 (0.0%) 0 (0.0%) 0 (0.0%) 3 (33.3%) 4 (50.0%) 0 (0.0%) 1 (12.5%) 0 (0.0%) 3 (37.5%) Hospitalisation due to COVID-19 † 9 (52.9%) 1 (11.1%) 8 (100%)

Fatal COVID-19 4 (23.5%) 0 (0.0%) 4 (50.0%) DN = diabetic nephropathy. GN = glomerulonephritis. HTN = hypertension. IQR = inter-1320 quartile range. Subsets defined according to peak WHO severity over the course of the illness. 1321 1322 1323

. CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted May 1, 2022. ; https://doi.org/10.1101/2022.04.29.22274267 doi: medRxiv preprint . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted May 1, 2022. ; https://doi.org/10.1101/2022.04.29.22274267 doi: medRxiv preprint . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted May 1, 2022. ; https://doi.org/10. 1101 

• Individual-level data for transcriptomics, proteomics and flow cytometry are available 1009 without restriction from Zenodo

1012 Acknowledgements 1013 The authors thank the patients who volunteered for this study and the staff at Imperial College 1014

Healthcare NHS Trust (the Imperial College Healthcare NHS Trust renal COVID-19 group and 1015 dialysis staff

1021 1022 We also acknowledge the efforts of renal specialist doctors in training for assistance with 1023 recruiting patients to this study

Funding statement 1029 1030 This research was partly funded by Community Jameel and the Imperial President's 1031 Excellence Fund and by a UKRI-DHSC COVID-19 Rapid Response Rolling Call 1032 (MR/V027638/1) (to JEP), and by funding from UKRI/NIHR through the UK Coronavirus 1033

We also acknowledge the National Institute for 1034

A) and the Sidharth Burman endowment. 1040 MCP is a Wellcome Trust Senior Fellow in Clinical Science (212252/Z/18/Z). NM-T and ES 1041 are supported by Wellcome Trust and Imperial College London Research Fellowships, and 1042 CLC by an Auchi Clinical Research Fellowship. 1043 1044 The funders had no role in study design, data collection and interpretation, or the decision to 1045 submit the work for publication

Novartis and Gyroscope; DCT reports speaker and consultancy fees from Astra-1053

Zeneca and Novartis; JEP has received travel and accommodation expenses and hospitality 1054 from Olink to speak at Olink-sponsored academic meetings (none within the past 5 years). 1055 None of the other authors have any interests to declare

Longitudinal immune profiling reveals key myeloid signatures 1428 associated with COVID-19

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

Single-cell multi-omics analysis of the immune response in 1433 COVID-19

Longitudinal Multi-omics Analyses Identify Responses of 1435

Longitudinal profiling of respiratory and systemic immune responses 1438 reveals myeloid cell-driven lung inflammation in severe COVID-19

Longitudinal analysis reveals that delayed bystander CD8+ T cell 1441 activation and early immune pathology distinguish severe COVID-19 from mild disease

A blood atlas of COVID-19 defines hallmarks of disease severity and 1444 specificity

Longitudinal proteomic analysis of severe COVID-19 reveals survival-1446 associated signatures, tissue-specific cell death, and cell-cell interactions

Longitudinal proteomic profiling of dialysis patients with COVID-19 1449 reveals markers of severity and predictors of death

Systems-Level Immunomonitoring from Acute to Recovery Phase 1452 of Severe COVID-19

SARS-CoV-2 RNAemia and proteomic trajectories inform 1457 prognostication in COVID-19 patients admitted to intensive care

Seroconversion stages COVID19 into distinct pathophysiological 1460 states

Proteomic Characterization of Acute Kidney Injury in Patients 1462 Hospitalized with SARS-CoV2 Infection

Circulating proteins to predict adverse COVID-19 outcomes

The COVIDome Explorer researcher portal

Dexamethasone in Hospitalized Patients with Covid-19

Interleukin-6 Receptor Antagonists in Critically Ill Patients with 1471

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

Baricitinib plus Remdesivir for Hospitalized Adults with Covid-19

COVID-19-related mortality in kidney transplant and haemodialysis 1480 patients: a comparative, prospective registry-based study

Immunogenicity Rates After SARS-CoV-2 Vaccination in People With 1484 End-stage Kidney Disease: A Systematic Review and Meta-analysis

Antibody Response to COVID-19 Vaccination in Patients Receiving 1487

GSVA: gene set variation analysis for 1489 microarray and RNA-Seq data

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

ISG15-dependent activation of the sensor MDA5 is antagonized by the 1494 SARS-CoV-2 papain-like protease to evade host innate immunity

Upregulated Explicitly in SARS-CoV-2 Infected Lung Epithelial Cells

WGCNA: An R package for weighted correlation network 1500 analysis

Identifying temporal 1502 and spatial patterns of variation from multi-modal data using MEFISTO

The 1505 prognostic value of S100A calcium binding protein family members in predicting severe 1506 forms of COVID-19

Reframing Immune-Mediated Inflammatory 1509 Diseases through Signature Cytokine Hubs

LRRC15 mediates an accessory interaction with the SARS-CoV-2 spike 1512 protein

LRRC15 is an inhibitory receptor blocking SARS-CoV-2 spike-mediated 1514 entry in trans

Outcomes of patients with end-stage kidney disease hospitalized with 1517 COVID-19

Regulation of DNA replication-coupled histone gene expression

Histone levels are regulated by 1521 phosphorylation and ubiquitylation-dependent proteolysis

The role of extracellular histone in organ injury

Proteomics of SARS-CoV-2-infected host cells reveals therapy 1526 targets

SARS-CoV-2 uses a multipronged strategy to impede host protein 1528 synthesis

Inhibition of Nox2 Oxidase Activity Ameliorates Influenza A Virus-1530 Induced Lung Inflammation

Reactive oxygen species delay control of lymphocytic choriomeningitis 1532 virus

EROS-mediated control of NOX2 and

Interleukin-11 signaling underlies fibrosis, 1536 parenchymal dysfunction, and chronic inflammation of the airway

IL-11 is a crucial determinant of cardiovascular fibrosis

The Role of the Adipokine Leptin in Immune Cell Function 1541 in Health and Disease

The glycoprotein-hormones activin A and inhibin A interfere with 1544 dendritic cell maturation

CCL22 controls immunity by promoting regulatory T cell communication 1547 with dendritic cells in lymph nodes

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

Autoantibodies against type I IFNs in patients with life-threatening 1552 COVID-19

Thromboembolic Risk in Hospitalized and Nonhospitalized COVID-19 1554 Patients: A Self-Controlled Case Series Analysis of a Nationwide Cohort

Current and novel biomarkers of thrombotic risk in COVID-19: a 1557 Consensus Statement from the International

Risks of deep vein thrombosis, pulmonary embolism, and bleeding 1560 after covid-19: nationwide self-controlled cases series and matched cohort study

Thrombotic Thrombocytopenia after

Thrombosis and Thrombocytopenia after

Insights in ChAdOx1 nCoV-19 vaccine-induced immune thrombotic 1567 thrombocytopenia

The nf-core framework for community-curated bioinformatics 1569 pipelines

Nextflow enables reproducible computational workflows

FastQC: a quality control tool for high throughput sequence data

Cutadapt removes adapter sequences from high-throughput sequencing 1575 reads

STAR: ultrafast universal RNA-seq aligner

HTSeq--a Python framework to work with high-1580 throughput sequencing data

edgeR: a Bioconductor package for 1583 differential expression analysis of digital gene expression data

A scaling normalization method for differential 1586 expression analysis of RNA-seq data

Genenames.org: the HGNC and VGNC resources in 2021. Nucleic 1591 Acids Res

M3C: Monte Carlo reference-based consensus clustering

Tissue-based map of the human proteome. Science (80-. )

Fitting Linear Mixed-Effects Models 1597 Using lme4

Tests in 1599 Linear Mixed Effects Models

Dream: powerful differential expression analysis for 1601 repeated measures designs

variancePartition: interpreting drivers of variation in 1604 complex gene expression studies

Molecular signatures database (MSigDB) 3.0

Type I interferons affect the metabolic fitness of CD8+ T cells from 1609 patients with systemic lupus erythematosus

Robust rank aggregation for gene list integration 1612 and meta-analysis

Application of Weighted Gene Co-expression Network Analysis for Data from 1615 Paired Design

Repeated measures correlation

Defining clusters from a hierarchical cluster 1619 tree: The Dynamic Tree Cut package for R

A review of spline 1622 function procedures in R

clusterProfiler 4.0: A universal enrichment tool for interpreting omics data

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

Robust enumeration of cell subsets from tissue expression 1630 profiles

Identifying temporal and spatial patterns of variation from multimodal 1632 data using MEFISTO

Building Predictive Models in R Using the caret Package

Regularization Paths for Generalized Linear 1636 Models via Coordinate Descent

the Wave 1 and Wave 2 cohorts. The column "Aggregated Score" represents the RRA score 1424 for the P-values from both cohorts. 1425