Multi-layered transcriptomic analyses reveal an immunological overlap between COVID-19 and hemophagocytic lymphohistiocytosis associated with disease severity

TL;DR: In this article, the authors show that COVID-19 and HLH have an overlap of signaling pathways and gene signatures commonly dysregulated, which were defined by investigating the transcriptomes of 1253 subjects (controls, COVID19, HLH patients) using microarray, bulk RNA-sequencing (RNAseq), and single-cell RNAseq (scRNAseq).

Abstract: Clinical and hyperinflammatory overlap between COVID-19 and hemophagocytic lymphohistiocytosis (HLH) has been reported. However, the underlying mechanisms are unclear. Here we show that COVID-19 and HLH have an overlap of signaling pathways and gene signatures commonly dysregulated, which were defined by investigating the transcriptomes of 1253 subjects (controls, COVID-19, and HLH patients) using microarray, bulk RNA-sequencing (RNAseq), and single-cell RNAseq (scRNAseq). COVID-19 and HLH share pathways involved in cytokine and chemokine signaling as well as neutrophil-mediated immune responses that associate with COVID-19 severity. These genes are dysregulated at protein level across several COVID-19 studies and form an interconnected network with differentially expressed plasma proteins which converge to neutrophil hyperactivation in COVID-19 patients admitted to the intensive care unit. scRNAseq analysis indicated that these genes are specifically upregulated across different leukocyte populations, including lymphocyte subsets and immature neutrophils. Artificial intelligence modeling confirmed the strong association of these genes with COVID-19 severity. Thus, our work indicates putative therapeutic pathways for intervention.

Summary

INTRODUCTION

• More than one year of Coronavirus disease 2019 (COVID-19) pandemic caused by the severe acute respiratory syndrome Coronavirus (SARS-CoV)-2, more than 197 million cases and 4,2 million deaths have been reported worldwide (July 30th 2021, WHO COVID-19 Dashboard).
• The clinical presentation ranges from asymptomatic to severe disease manifesting as pneumonia, acute respiratory distress syndrome (ARDS), and a life-threatening hyperinflammatory syndrome associated with excessive cytokine release 1–3.
• Thus, the need for a better understanding of the immunopathologic mechanisms associated with severe SARS-CoV-2 infection.
• Thus, HLH has been proposed as an underlying etiologic factor of severe COVID191,3,20.
• The authors found shared gene signatures and cellular signaling pathways involved in cytokine and chemokine signaling as well as neutrophil- mediated immune responses that associate with COVID-19 severity.

RESULTS

• The transcriptional overlap between COVID-19 and HLH Considering the fact that COVID-19 and HLH14,21 share clinical and biological features the authors suspected that these diseases have a transcriptional overlap underlying their phenotypic convergence.
• Among them are important genes encoding molecules involved in activation of inflammatory immune responses (e.g., PGLYRP1, OLR1, FFAR2), cytokine and chemokine mediated immune pathways (e.g., IL1R2, CXCR2, CXCR8, CCL4, CCL2), and neutrophil activation (e.g., CD177, MPO, ELANE).
• Bivariate correlation analysis showed a similar phenomenon .
• Altogether, these multi-layered transcriptomic results associate COVID-19 and HLH common genes with disease severity.
• In agreement with the recent observation that neutrophil hyperactivation plays a key role in the severity of COVID-1981–84 and HLH18,19, their approach indicates that COVID-19 and HLH have a common transcriptional profile formed predominantly by a group of regulatory molecules related to cytokine/chemotaxis and by a group of effector molecules that are linked to neutrophil hyperactivation and disease severity.

ACKNOWLEDGMENTS

• The authors acknowledge the Latin American Society of Immunodeficiencies for providing the research funding of LFS ( Fellowship award 2020), and the São Paulo Research Foundation (FAPESP grants.
• This study was financed in part by the coordination for the improvement of higher education personnel – Brazil – finance code 001.
• The authors thank Prof. Luis Carlos de Souza Ferreira for discussion and suggestions to develop this project.

FIGURE LEGENDS:

• Schematic view of datasets and results obtained indicating new candidate biomarkers and therapeutic targets for COVID-19 and HLH, also known as Graphic Abstract.
• The label (gene name) colours represent transcripts from Overlap 1 , Overlap 2 (red), and Overlap 3 (blue).
• Correlation coefficient and significance level (p-value) for each correlation are shown within each graph.
• Multi-layered transcriptomic analysis associates COVID-19 and HLH common genes with disease severity, also known as Figure 5.
• (A) Schematic overview of sample cohort and classification of scRNAseq dataset obtained by Schulte-Schrepping et al.65 and used for the following analysis.

Data and code availability

• This paper analyzes existing, publicly available data.
• The accession numbers for the datasets are listed in the key resources table.
• All original codes used for data analysis have been deposited at github (https://github.com/lschimke/COVID19-and-HLH-paper) and are publicly available as of the date of publication.
• R packages are listed in the key resources table.
• Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Data Curation

• The authors searched in public functional genomics data repositories (Gene Expression Omnibus111 and Array Express112 for human transcriptome data from patients with HLH and COVID-19 published until February 2021.
• This resulted in 7 studies, 6 from COVID-19 patients and one from HLH patients (containing 4 patients with germline mutations and 7 without identified mutation) with transcriptome data generated from different platforms (Table 1).
• Combined, the authors investigated a total of 1253 transcriptome samples.
• Statistical cut-offs to define DEGs are described below in the section statistical analysis.
• Regression for the number of UMIs and scaling were performed as previously described65.

Interactome analysis

• For more comprehensive Protein–Protein Interaction (PPI) analyses, the authors used NAViGaTOR 3.0.14120 to visualize genes commonly dysregulated in COVID-19 and HLH datasets, highlighting the biological processes enriched by each gene.
• Then, the resultant network was annotated, analysed, and visualized using NAViGaTOR 3.0.14120.
• Enrichment analysis and data visualization Elsevier Pathway Collection analysis for selected gene lists (7 genes underlying fHLH/IEI and 11 genes associated with severe COVID-19) was carried out using Enrichr webtool123–125.
• Circular heatmaps were generated using R version 4.0.5 (The R Project for Statistical Computing.

Correlation Analysis

• Principal Component Analysis (PCA) of genes associated with COVID-19 severity (25 transcripts) was performed using the R functions prcomp and princomp, through factoextra package127.
• Canonical Correlation Analysis (CCA)128 of genes associated with cytokines/chemokines and neutrophil-mediated immune responses was performed using the packages CCA and whitening128.
• In addition, multilinear regression analysis for combinations of different variables was performed using the R package ggpubr, ggplot2 and ggExtra.

Proteome Data Analysis

• The authors also evaluated the proteomics data obtained from plasma samples of COVID-19 patients previously reported by Overmyer et al.58.
• Briefly, raw LFQ abundance values were quantified, normalized and log2 transformed as previously described58.
• Differences in protein expression between COVID-19_ICU and COVID-19_nonICU were calculated as described below in the section statistical analysis.
• Follow-up analysis used the Gini decrease, number of nodes, and mean minimum depth as criteria to determine variable importance.
• The adequacy of the Random Forest model as a classifier was assessed through out of bags error rate and ROC curve.

STATISTICAL ANALYSIS

• The authors used the Fisher´s method to combine p values from multiple studies for information integration114.
• Differences in protein expression between COVID-19_ICU and COVID-19_nonICU was calculated using the nonparametric MANOVA (multivariate analysis of variance) test129 followed by analysis of nonparametric Inference for Multivariate Data130 using the R packages npmv, nparcomp, and ggplot2.
• All supplemental figures, titles, and legends are provided in separate document file.

Full text

Multi-layered transcriptomic analyses reveal an immunological overlap between COVID-19
and hemophagocytic lymphohistiocytosis associated with disease severity
Lena F. Schimke
a,5
, Alexandre H.C. Marques
a
, Gabriela Crispim Baiocchi
a
, Caroline Aliane de Souza Prado
b
,
Dennyson Leandro M. Fonseca
b
, Paula Paccielli Freire
a
, Desirée Rodrigues Plaça
b
, Igor Salerno Filgueiras
a
,
Ranieri Coelho Salgado
a
, Gabriel Jansen-Marques
c
, Antonio Edson Rocha Oliveira
b
, Jean Pierre
Schatzmann Peron
a
, José Alexandre Marzagão Barbuto
a,d
, Niels Olsen Saraiva Camara
a
, Vera Lúcia Garcia
Calich
a
, Hans D. Ochs
e
, Antonio Condino-Neto
a
, Katherine A. Overmyer
f,g
, Joshua J. Coon
h,i
, Joseph
Balnis
j,k
, Ariel Jaitovich
j,k
, Jonas Schulte-Schrepping
l
, Thomas Ulas
m
, Joachim L. Schultze
l,m
, Helder I.
Nakaya
b
, Igor Jurisica
n,o,p
, Otavio Cabral-Marques
a,b,q
a
Departamento de Imunologia, Instituto de Ciencias Biomedicas, Universidade de São Paulo, São Paulo,
SP, Brazil. ID Ringgold USP + ICB: 54544, ISNI: 0000000406355304
b
Department of Clinical and Toxicological Analyses, School of Pharmaceutical Sciences, University of São
Paulo, São Paulo, SP, Brazil.
c
Information Systems, School of Arts, Sciences and Humanities, University of Sao Paulo, São Paulo, SP,
Brazil.
d
Laboratory of Medical Investigation in Pathogenesis and targeted therapy in Onco-Immuno-Hematology
(LIM-31), Department of Hematology, Hospital das Clínicas HCFMUSP, Faculdade de Medicina,
Universidade de Sao Paulo, Sao Paulo, Brazil.
e
Department of Pediatrics, University of Washington School of Medicine, and Seattle Children's Research
Institute, Seattle, WA, USA.
f
National Center for Quantitative Biology of Complex Systems, Madison, WI 53562, USA
g
Morgridge Institute for Research, Madison, WI 53562, USA
h
Department of Biomolecular Chemistry, University of Wisconsin, Madison, WI 53506, USA Department
of
i
Chemistry, University of Wisconsin, Madison, WI 53506, USA.
j
Division of Pulmonary and Critical Care Medicine, Albany Medical Center, Albany, NY 12208, USA
k
Department of Molecular and Cellular Physiology, Albany Medical College, Albany, NY 12208, USA
l
Life and Medical Sciences (LIMES) Institute, University of Bonn, Germany.
m
German Center for Neurodegenerative Diseases (DZNE), PRECISE Platform for Genomics and
Epigenomics at DZNE, and University of Bonn, Bonn, Germany.
n
Osteoarthritis Research Program, Division of Orthopedic Surgery, Schroeder Arthritis Institute and Data
Science Discovery Centre for Chronic Diseases, Krembil Research Institute, University Health Network,
Toronto, Canada.
o
Departments of Medical Biophysics and Computer Science, University of Toronto, Toronto, Canada.
p
Institute of Neuroimmunology, Slovak Academy of Sciences, Bratislava, Slovakia.
q
Network of Immunity in Infection, Malignancy, and Autoimmunity (NIIMA), Universal Scientific
Education and Research Network (USERN), Sao Paulo, Brazil.
5
Lead contact
Corresponding authors:
Lena F. Schimke, MD
Department of Immunology
Institute of Biomedical Sciences - University of
São Paulo
Lineu Prestes Avenue, 1730,
São Paulo, SP, CEP 05508-900, Brazil
e-mail: lenaschimke@hotmail.com
phone.: +55 (11) 3091-7387
Otavio Cabral-Marques, MSc, PhD
Department of Immunology
Institute of Biomedical Sciences - University of
São Paulo
Lineu Prestes Avenue, 1730,
São Paulo, SP, CEP 05508-900, Brazil
e-mail: otavio.cmarques@usp.br
phone: +55 11 974642022
.CC-BY-NC-ND 4.0 International licenseavailable under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprintthis version posted August 1, 2021. ; https://doi.org/10.1101/2021.07.30.454529doi: bioRxiv preprint

ABSTRACT
Clinical and hyperinflammatory overlap between COVID-19 and hemophagocytic
lymphohistiocytosis (HLH) has been reported. However, the underlying mechanisms are unclear.
Here we show that COVID-19 and HLH have an overlap of signaling pathways and gene
signatures commonly dysregulated, which were defined by investigating the transcriptomes of
1253 subjects (controls, COVID-19, and HLH patients) using microarray, bulk RNA-sequencing
(RNAseq), and single-cell RNAseq (scRNAseq). COVID-19 and HLH share pathways involved in
cytokine and chemokine signaling as well as neutrophil-mediated immune responses that
associate with COVID-19 severity. These genes are dysregulated at protein level across several
COVID-19 studies and form an interconnected network with differentially expressed plasma
proteins which converge to neutrophil hyperactivation in COVID-19 patients admitted to the
intensive care unit. scRNAseq analysis indicated that these genes are specifically upregulated
across different leukocyte populations, including lymphocyte subsets and immature neutrophils.
Artificial intelligence modeling confirmed the strong association of these genes with COVID-19
severity. Thus, our work indicates putative therapeutic pathways for intervention.
Keywords:
COVID-19, Hemophagocytic lymphohistiocytosis, common signaling pathways, gene signatures
.CC-BY-NC-ND 4.0 International licenseavailable under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprintthis version posted August 1, 2021. ; https://doi.org/10.1101/2021.07.30.454529doi: bioRxiv preprint

INTRODUCTION
More than one year of Coronavirus disease 2019 (COVID-19) pandemic caused by the severe
acute respiratory syndrome Coronavirus (SARS-CoV)-2, more than 197 million cases and 4,2
million deaths have been reported worldwide (July 30
th
2021, WHO COVID-19 Dashboard). The
clinical presentation ranges from asymptomatic to severe disease manifesting as pneumonia,
acute respiratory distress syndrome (ARDS), and a life-threatening hyperinflammatory
syndrome associated with excessive cytokine release (hypercytokinaemia)
1–3
. Risk factors for
severe manifestation and higher mortality include old age as well as hypertension, obesity, and
diabetes
4
. Currently, COVID-19 continues to spread, new variants of SARS-CoV-2 have been
reported and the number of infections resulting in death of young individuals with no
comorbidities has increased the mortality rates among the young population
5,6
. In addition,
some novel SARS-CoV-2 variants of concern appear to escape neutralization by vaccine-induced
humoral immunity
7
. Thus, the need for a better understanding of the immunopathologic
mechanisms associated with severe SARS-CoV-2 infection.
Patients with severe COVID-19 have systemically dysregulated innate and adaptive
immune responses, which are reflected in elevated plasma levels of numerous cytokines and
chemokines including granulocyte colony-stimulating factor (GM-CSF), tumor necrosis factor
(TNF), interleukin (IL)-6, IL-6R, IL18, CC chemokine ligand 2 (CCL2) and CXC chemokine ligand 10
(CXCL10)
8–10
, and hyperactivation of lymphoid and myeloid cells
11
. Notably, the
hyperinflammation in COVID-19 shares similarities with cytokine storm syndromes such as those
triggered by sepsis, autoinflammatory disorders, metabolic conditions and malignancies
12–14
,
often resembling a hematopathologic condition called hemophagocytic lymphohistiocytosis
(HLH)
15
. HLH is a life-threatening progressive systemic hyperinflammatory disorder
characterized by multi-organ involvement, fever flares, hepatosplenomegaly, and cytopenia due
to hemophagocytic activity in the bone marrow
15–17
or within peripheral lymphoid organs such
as pulmonary lymph nodes and spleen. HLH is marked by aberrant activation of B and T
lymphocytes and monocytes/macrophages, coagulopathy, hypotension, and ARDS. Recently,
neutrophil hyperactivation has been shown to also play a critical role in HLH development
18,19
.
This is in agreement with the observation that the HLH-like phenotype observed in severe
COVID-19 patients is due to an innate neutrophilic hyperinflammatory response associated with
.CC-BY-NC-ND 4.0 International licenseavailable under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprintthis version posted August 1, 2021. ; https://doi.org/10.1101/2021.07.30.454529doi: bioRxiv preprint

virus-induced hypercytokinaemia which is dominant in patients with an unfavorable clinical
course
17
. Thus, HLH has been proposed as an underlying etiologic factor of severe COVID-
19
1,3,20
. HLH usually develops during the acute phase of COVID-19
1,20–27
. However, a case of HLH
that occurred two weeks after recovery from COVID-19 has recently been reported as the cause
of death during post-acute COVID-19 syndrome
28
.
The familial form of HLH (fHLH) is caused by inborn errors of immunity (IEI) in different
genes encoding proteins involved in granule-dependent cytotoxic activity of leukocytes such as
AP3B1, LYST, PRF1, RAB27A, STXBP2, STX11, UNC13D
29–31
. In contrast, the secondary form
(sHLH) usually manifests in adults following a viral infection (e.g., adenovirus, EBV, enterovirus,
hepatitis viruses, parvovirus B19, and HIV)
32,33
, or in association with
autoimmune/rheumatologic, malignant, or metabolic conditions that lead to defects in T/NK cell
functions and excessive inflammation
16,31
. fHLH and sHLH affect both children and adults,
however, the clinical and genetic distinction of HLH forms is not clear since immunocompetent
children can develop sHLH
34,35
, while adult patients with sHLH may also have germline
mutations in HLH genes
36
. Of note, germline variants in UNC13D and AP3B1 have also been
identified in some COVID-19 patients with HLH phenotype
37
, thus, indicating that both HLH
forms may be associated with COVID-19.
Here, we characterized the signaling pathways and gene signatures commonly
dysregulated in both COVID-19 and HLH patients by investigating the transcriptomes of 1253
subjects (controls, COVID-19, and HLH patients) assessed by microarray, bulk RNA-sequencing
(RNAseq), and single-cell RNAseq (scRNAseq) (Table 1). We found shared gene signatures and
cellular signaling pathways involved in cytokine and chemokine signaling as well as neutrophil-
mediated immune responses that associate with COVID-19 severity.
.CC-BY-NC-ND 4.0 International licenseavailable under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprintthis version posted August 1, 2021. ; https://doi.org/10.1101/2021.07.30.454529doi: bioRxiv preprint

RESULTS
The transcriptional overlap between COVID-19 and HLH
Considering the fact that COVID-19 and HLH
14,21
share clinical and biological features (Figure 1A)
we suspected that these diseases have a transcriptional overlap underlying their phenotypic
convergence. First, we obtained differentially expressed genes (DEGs) for each dataset from
different high-throughput transcriptome technologies present in peripheral blood lymphocytes
(PBLs), peripheral blood mononuclear cells (PBMCs), and nasopharyngeal swabs from COVID-19
patients, HLH patients and controls
(Figure 1B and Supp. Table S1). Then, we defined the
transcriptome overlap between DEGs from COVID-19 cohorts and HLH by cross-technology
comparisons, performed enrichment analysis and association studies between specific DEGs and
severity status of disease (Figure 1C). To identify the common DEGs we divided the datasets into
three subgroups based on type of samples and RNAseq platforms: Overlap 1 (HLH and COVID-19
blood transcriptomes), Overlap 2 (HLH and COVID-19 nasopharyngeal swab transcriptomes),
and Overlap 3 (HLH and COVID-19 scRNAseq transcriptomes) (Supp. Figure 1A and 1B, and
Supp. Table S2 and S3). We found a total of 239 unique common DEGs between HLH and all
COVID-19 datasets, most of them (237 DEGs) up-regulated (Figure 1D). Hereafter, we focused
on the implications of the up-regulated genes, since the 2 common down-regulated genes
(granulysin or GNLY; myomesin 2 or MYOM2) alone did not enrich any significant pathway.
However, this might also indicate a defect in cytotoxicity activity, typical of HLH
31
, that will
require future investigation. The 237 common up-regulated DEGs encode proteins mainly
involved in immune system, metabolic and signaling processes, forming a highly connected
biological network based on physical protein-protein interactions (PPI, Figure 1E). Among them
are important genes encoding molecules involved in activation of inflammatory immune
responses (e.g., PGLYRP1, OLR1, FFAR2), cytokine and chemokine mediated immune pathways
(e.g., IL1R2, CXCR2, CXCR8, CCL4, CCL2), and neutrophil activation (e.g., CD177, MPO, ELANE). Of
note, the transcriptional overlap between HLH and COVID-19 contains several molecules
interacting with 7 genes causing fHLH due to IEI which itself were not among our DEGs (Figure
1E).
.CC-BY-NC-ND 4.0 International licenseavailable under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprintthis version posted August 1, 2021. ; https://doi.org/10.1101/2021.07.30.454529doi: bioRxiv preprint