Single-cell RNA-seq reveals profound monocyte changes in Paediatric Inflammatory Multisystem Syndrome Temporally associated with SARS-CoV-2 infection (PIMS-TS)

Paediatric inflammatory multisystem inflammatory syndrome temporally associated with SARS-CoV-2 infection (PIMS-TS) is a new disease with overlapping features of Kawasaki disease (KD) and toxic shock syndrome. Unbiased single cell RNA sequencing analysis of peripheral blood mononuclear cells from PIMS-TS and KD patients shows monocytes are the main source of pro-inflammatory cytokines and large changes in the frequency of classical, intermediate and non-classical monocytes occur in both diseases.


MAIN
SARS-CoV-2, the viral cause of coronavirus disease 2019 , usually causes mild respiratory infection in children and young adults [1]. Multiple reports show some children subsequently develop a severe hyperinflammatory syndrome comprising fever, inflammation and evidence of single or multi-organ failure that manifests with cardiac dysfunction, hypotension and life-threatening shock [2,3,4]. This new disease, called Paediatric inflammatory multisystem inflammatory syndrome temporally associated with SARS-CoV-2 infection (PIMS-TS, also known as multisystem inflammatory syndrome in children MIS-C) occurs several weeks after SARS-CoV-2 infection which may or may not have been symptomatic. PIMS-TS shares clinical features with Kawasaki Disease (KD) a vasculitis affecting mostly children <5 years old that can present with features of systemic and cardiac inflammation [5,6]. KD can lead to fatal coronary artery aneurysms and is the leading cause of acquired heart disease in children in developed nations. This risk is significantly reduced with early treatment using intravenous immunoglobulin (IVIG) although 20% of patients require additional interventions due to failure to respond to initial treatment [6]. The immunology of KD is still not fully understood but patient data and animal models suggest monocytes are involved [5,7,8]. The CAWS mouse model of KD shows increased recruitment and accumulation of bone marrow derived monocytes into the heart and the inflammatory cytokine IL1-beta are both essential for pathogenesis [8]. Given the recent emergence of PIMS-TS, little is known about the immunological processes driving the disease and whether these are similar or different to KD.
We recruited nine children with signs of PIMS-TS or KD who presented at Birmingham Children's Hospital from April to May 2020. Demographic and clinical features are presented in Figure 1 and SI Table 1  admission whereas 5/7 PIMS-TS patients required PICU for 3 to 8 days ( Figure 1B). A range of treatment intensities were required to achieve clinical response. All PIMS-TS patients were lymphopaenic and four had elevated neutrophils ( Figure 1C). Frequencies of T and B cells were variable but the frequency and absolute counts of NK cells were low in PIMS-TS patients compared to KD ( Figure 1D).
To resolve the complexity and heterogeneity of cellular immune subsets in the blood of patients with PIMS-TS, we performed unbiased single cell RNA sequencing (scRNAseq) on the pre-treatment peripheral blood mononuclear cells (PMBCs) of two representative patients . CC-BY 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted August 7, 2020. . https://doi.org/10.1101/2020 with PIMS-TS, P6 and P7. Both were admitted to PICU but patient P7 responded rapidly to one cycle of IVIG whereas patient P6 required additional IVIG, steroids and Tocilizumab. A convalescent sample from P6 was also analysed. For comparative purposes we analysed pretreatment PBMCs from a KD patient (KD2) who similarly required multiple treatment cycles ( Figure 1B). Unsupervised clustering of 10,031 cells produced 19 different clusters comprising all major lymphocyte subsets ( Figure 1E, SI Fig 1 and 2). Examining cytokine mRNA data showed monocytes were by far the predominant cell type expressing IL1-beta, TNF-a, IL-8, IL-18 and IL1-RA (SI Fig 3). A small number of B-cells positive for IL6, monocytes positive for IL10 and NK cells positive for IFN-g mRNA were present but we did not detect any robust evidence of T-cell activation in our data. Expression of KI-67, a marker of replicating cells, was rarely detected.
Analysing each patient separately, we observed that all possessed lymphocytes assigned to one of the nine B-cell, T-cell or NK cell clusters, although the frequencies varied between patients and for P6 from acute disease to convalescence ( Figure 1F and 1G). In contrast, each patient's monocytes were assigned to only one or two clusters specific to each patient/disease state. To explore this in greater detail we repeated dimensionality reduction analyses on just the subset of data representing monocytes. This yielded nine monocytes clusters (numbered mc0 to mc8, Figure 2A) that we aligned to the three canonical monocyte subsets based on CD14 and CD16 expression (SI FIG 4 and 5). Clusters mc0, mc1, mc3, mc4 and mc5 were designated CD14 + CD16classical monocytes (CM), highly phagocytic scavenger cells that comprise 80-90% of monocytes in blood [9]. Cluster mc6 contained CD14 int CD16 int intermediate monocytes (IM) that play diverse roles in inflammation and angiogenesis [9,10]. Cluster mc2 contained CD14 lo CD16 hi non-classical monocytes (NCM), also called patrolling monocytes, that actively surveil the vasculature via their expression of CX3CR1 and exert pro-or anti-inflammatory roles in different contexts [9,11]. CX3CR1 was strongly expressed by cells in mc2 (SI Fig 5).
The IM and NCM subsets are normally rare in blood [9]. The reclustering analysis ( Figure 2B) showed the acute sample from KD patient KD2 had high frequencies of IM (8% of total monocytes, cluster mc6) and NCM (38% of total monocytes, cluster mc2). The PIMS-TS patient P6, who like KD2 failed to respond to IVIG, also had an abnormally high frequency of NCM (34% of total monocytes) but lacked IM; upon recovery their monocytes returned to the normal state with only CM detected in their convalescent sample. The acute sample from patient P7, who unlike the other patients responded rapidly to IVIG, lacked IM or NCM. The division of CM into five distinct clusters that varied by disease type and patient severity showed this monocyte subtype is also altered during KD and PIMS-TS. Interestingly, clusterin (CLU), a protein with multiple roles in autoimmunity and inflammation [12], was highly . CC-BY 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted August 7, 2020. . https://doi.org/10.1101/2020.08.06.20164848 doi: medRxiv preprint expressed by CMs during the acute phase of KD and PIMS-TS but not in convalescence Figure 2C, SI Fig 4 and 5). Clusterin has not been reported to be expressed by monocytes [9,13,14].
Our study is the first to reveal that KD and PIMS-TS are both characterised by profound changes in monocyte gene expression and frequencies of key monocyte subpopulations.
While supported by a range of known marker genes, our results require confirmation in a larger patient cohort before generalisation. Nevertheless, given the well-known role of monocytes in KD animal models [5,8]   is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted August 7, 2020. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted August 7, 2020. . https://doi.org/10.1101/2020.08.06.20164848 doi: medRxiv preprint is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted August 7, 2020. . https://doi.org/10.1101/2020.08.06.20164848 doi: medRxiv preprint is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted August 7, 2020. . https://doi.org/10.1101/2020.08.06.20164848 doi: medRxiv preprint

Study design, participants and ethical approval
All samples were obtained from children at Birmingham Children's Hospital as part of a Health Research Authority approved study (TrICICL) in accordance with the Declaration of Helsinski.
The study was reviewed and approved by South of Birmingham Research Ethics Committee

ScRNAseq Data Pre-processing and Quality Control
Processing of raw reads including 10X barcode-aware demultiplexing from BCL to FASTQ files, transcriptome alignment to human genome assembly GRCh38 and unique molecular . CC-BY 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted August 7, 2020. . https://doi.org/10.1101/2020.08.06.20164848 doi: medRxiv preprint identifier (UMI) counting were performed using the 10X Cell Ranger pipeline version 3.1.0 with GRCh38-version 3.0.0 as the reference. Seurat v3.1.5 [21] was used for sample merging, quality control (QC), clustering and reporting. Both the gene expression (RNA) and antibodyderived tag (ADT) assays were loaded from the CellRanger count platform into a Seurat object for each sample excluding cells with less than 200 genes and features detected in less than 3 cells. The four samples were then merged to create a single aggregated object. QC was conducted on both assays separately. For the RNA assay, to mitigate the influence of sex on clustering and differential expression, the XIST and RSP4Y1 genes were removed. Cells with a number of features between 200 and 600, counts of greater than 1000 and mitochondrial percentage less than 15% were kept for further processing. The gene expression data was normalised using the 'LogNormalize' method with the scale.factor set to the default 10000.

Monocyte Clustering
To examine the monocyte populations in finer detail, we extracted clusters 4, 6, 9, 10 and 13 from Figure 1E using the 'Subset' function. As monocytes were abundant within each sample, all monocyte clusters were included (2248 monocytes extracted). Both assays were normalized separately again (RNA: 'LogNormalize' & ADT: 'CLR'). The 'vst' method was used to find the 2000 most variable features within the RNA assay and the data was scaled with the variables 'number of counts' and 'percent mitochondria' regressed. The ADT assay was scaled with default parameters. In this instance, 100 principal components were used in the 'RunPCA' function with the top 50 of those components used to find cell neighbours and clusters. The is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted August 7, 2020. . https://doi.org/10.1101/2020.08.06.20164848 doi: medRxiv preprint UMAP reduction was selected here as the clusters visually displayed a better path from classical to intermediate to non-classical monocytes over the tSNE reduction. Nine clusters were originally generated with further assessment resulting in the identification of seven monocyte clusters. CD14 and CD16 markers from the ADT assay were used to classify classical (CD14+), non-classical (CD16+) and intermediate (CD14+CD16+)

Reporting Summary
Further information on research design is available in the research reporting summary linked to this article.

Acknowledgements
We thank the children and their families for consenting to join this research study and Birmingham Women's and Children's Hospital Charity for funding the single cell RNA sequencing analysis. ES is supported by a CRUK Birmingham Centre clinical PhD fellowship.
BS is funded by a National Institute for Health Research (NIHR) Clinician Scientist fellowship.
The views expressed are those of the author(s) and not necessarily those of the NHS, the NIHR or the Department of Health and Social Care. We thank Genomics Birmingham (University of Birmingham) for generating the single cell sequencing data. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted August 7, 2020. . https://doi.org/10.1101/2020 We would like to acknowledge staff at the Birmingham Women's and Children's Hospital NHS is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted August 7, 2020. . https://doi.org/10.1101/2020 SI is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted August 7, 2020. . https://doi.org/10.1101/2020 SI Table 2 is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted August 7, 2020. . https://doi.org/10.1101/2020.08.06.20164848 doi: medRxiv preprint SI Figure 1. Gene expression profile of 19 clusters generated from patient PBMC samples.
. CC-BY 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted August 7, 2020. . https://doi.org/10.1101/2020.08.06.20164848 doi: medRxiv preprint SI Figure 3. RNA expression and protein levels in patient PBMCs. tSNE representation of combined PBMCs from all four patient samples coloured by the level of RNA per cell or level of protein per cell measured using antibody derived tags.
. CC-BY 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted August 7, 2020. . https://doi.org/10.1101/2020.08.06.20164848 doi: medRxiv preprint SI Figure 4. Heatmap showing scaled expression of discriminative genes sets for each of the new clusters generated by reclustering of cells from the originally identified monocyte clusters. Only the top five discriminative genes per cluster are shown.  RETN  CD163  NKG7  RNASE2  CLU  IER3  ID2  HLA−DQB1  LGALS2  HBA2  HBB  HBA1  S100A12  HP  C1QB  C1QA  IFITM3  CDKN1C  FCGR3A  CRIP1  TMEM170B  IRS2  AREG  AC020656.1  NFKBIA  RNASE1  KLF10  LYZ   . CC-BY 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted August 7, 2020. . https://doi.org/10.1101/2020.08.06.20164848 doi: medRxiv preprint SI Figure 5. RNA expression and protein levels in patient monocytes.
UMAP representation of monocytes from all four patient samples coloured by the level of RNA per cell or level of protein per cell measured using antibody derived tags. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted August 7, 2020. . https://doi.org/10.1101/2020.08.06.20164848 doi: medRxiv preprint . CC-BY 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted August 7, 2020. . https://doi.org/10.1101/2020.08.06.20164848 doi: medRxiv preprint