The Notch1/CD22 signaling axis disrupts Treg function in SARS-CoV-2–associated multisystem inflammatory syndrome in children

Multisystem inflammatory syndrome in children (MIS-C) evolves in some pediatric patients following acute infection with SARS-CoV-2 by hitherto unknown mechanisms. Whereas acute-COVID-19 severity and outcomes were previously correlated with Notch4 expression on Tregs, here, we show that Tregs in MIS-C were destabilized through a Notch1-dependent mechanism. Genetic analysis revealed that patients with MIS-C had enrichment of rare deleterious variants affecting inflammation and autoimmunity pathways, including dominant-negative mutations in the Notch1 regulators NUMB and NUMBL leading to Notch1 upregulation. Notch1 signaling in Tregs induced CD22, leading to their destabilization in a mTORC1-dependent manner and to the promotion of systemic inflammation. These results identify a Notch1/CD22 signaling axis that disrupts Treg function in MIS-C and point to distinct immune checkpoints controlled by individual Treg Notch receptors that shape the inflammatory outcome in SARS-CoV-2 infection.


Experimental Procedures
Single-cell RNA sequencing. Cryopreserved PBMCs were thawed in plain RPMI (HyClone) pre-warmed to 37°C, washed in PBS (HyClone) and resuspended in FACS buffer (PBS with 1.5% FBS (Genesee Scientific) and 2.5 mM EDTA (Invitrogen)) for CD4 T cell enrichment through negative selection (Miltenyi Biotec). Samples were studied in 2 independent experiments: experiment 1 included 3 pediatric controls, 1 pre-treatment MIS-C patient, and 4 post-treatment MIS-C patients; experiment 2 included 1 pediatric control, 2 pre-treatment MIS-C patients, and 1 post-treatment MIS-C patient.
CD4 T cells isolated from each sample were stained with Hashing antibodies targeting CD298 and β2 microglobulin (BioLegend TotalSeq™-C anti-human Hashtag, Clones: LNH-94 and 2M2) for subsequent sample identification. To that end, cells were spun at 500 g for 7 min at 4°C, resuspended into 75 μL Fc block (BioLegend, Cat. No. 422302, 1:20 dilution) and incubated for 10 min at 4°C. 75 μL of hashing antibody (6.7 μg/mL working concentration) were then added and samples were incubated for 30 min at 4°C, with gentle resuspension midway. Cells were washed three times in Hash Staining buffer (BioLegend, Cat. No. 420201) and resuspended in PBS with 0.4% BSA (Sigma, Cat. No. A7030) at a concentration of 1,000 cells per μL. Finally, samples were pooled in equal ratio (7,500 cells per sample) for further processing, thus limiting technical batch effects (1). For the first experiment, the total of 60,000 cells pooled from 8 samples were split across two 10x Genomics chip channels, while the second experiment (30,000 cells from 4 samples) was loaded on one single channel. Cells were encapsulated, barcoded and lysed to enable the generation of cDNA libraries for transcriptome and HTO sequencing using the 10x Genomics technology (2). Libraries were sequenced on an Illumina NovaSeq 6000.
Single-cell RNA sequencing clustering analyses. Sequencing data from each 10x run were processed with the CellRanger pipeline (10x Genomics) for demultiplexing and gene alignment (2). The resulting raw count matrices were imported in R (v4.0.2 and above) using Seurat (v4.0.3) (3). Data from all 3 runs were merged into one Seurat object. Genes detected in <1 per 10,000 cells were filtered out, leaving a transcriptomic coverage of 21,675 genes. High quality cells with >1400 unique molecular identifiers (UMIs), >700 genes, a log10(gene) to log10(UMI) ratio >0.84 and mitochondrial to nuclear gene ratio <0.08 were retained for downstream analyses. Quality control revealed no significant batch effect: similar distributions were observed for the metrics mentioned above across different runs and experiments.
HTO data were normalized using centered log ratios before applying Seurat::HTODemux() with the clara method and a positive quantile cutoff of 0.98. Doublets and cells with unclear HTO assignment were excluded (Stoeckius et al., 2018).
Transcriptomic data for the remaining cells were normalized using Seurat:: SCTransform() and regressing out the effects of the mitochondrial gene ratio, number of UMIs and number of genes detected. Principal components (PCs) were calculated from the top 2000 variable features to reduce the data before mapping to a reference PBMC dataset using Azimuth (3). Cells mapped to CD4 T cell subsets were retained while contaminating lymphocytes were excluded, leaving a total of 29,754 Azimuth-annotated CD4 T cells for downstream analyses. SCT normalization and PCs calculation were repeated at this stage, to account for the top 3000 variable features after exclusion of TRAV, TRAJ, TRBV and TRBJ genes, thereby enabling cell clustering by transcriptomic profile independently of clonal identity. To control for inter-sample variability, the data were integrated by source sample using Harmony (4). Uniform manifold approximation and projection (UMAP) coordinates were then computed from the first 50 components of the harmony reduction, and graph-based clustering analysis was run on the first 40 components using Seurat::FindNeighbors() and Seurat::FindClusters(). A resolution of 0.6 was retained to define clusters. Seurat-defined clusters were manually annotated, with an initial coarse characterization based on the abundance of cells classified as naïve versus effector/memory by Azimuth. In parallel, genes significantly upregulated in each cluster were identified with Seurat::FindAllMarkers() using the Wilcoxon rank sum test.
Significance was defined as a p-value of <0.05 and a log2 fold change (LFC) in gene expression >0.25. Heatmaps were generated with Seurat::DoHeatmap() applied to the scaled SCT data on a random subsample of 100 cells per cluster. To confirm upregulation of NF-kB genes, a signature score for the TNFa signaling via NF-kB geneset, sourced from the MSigDB Hallmark collection (5), was calculated at the single-cell level using Seurat::AddModuleScore().
Pseudobulk differential expression analyses (DEA). For pseudobulk differential expression analyses (DEA), gene expression level data was aggregated at the patient level for each subset of interest, namely Tregs and activated Tconv. For this analysis, we considered as Treg any cell assigned to Seurat Cluster 15 (FOXP3-expressing cells) or annotated as Treg by Azimuth, which added up to 1,925 Treg cells across all 12 patients.
Similarly, we considered as activated Tconv any cell assigned to Seurat Clusters 9 to 14 and annotated as CD4 TCM, CD4 TEM, CD4 CTL or CD4 Proliferating by Azimuth (6,674 cells). PC analyses of the aggregated transcriptomic data highlighted healthy control 4 as a strong outlier among both Treg and activated Tconv subsets, leading us to exclude this patient from pseudobulk DEA. Independent pairwise analyses contrasting each of the 3 patient groups (MIS-C pre-treatment, MIS-C post-treatment and control) were run using Heatmaps of gene expression for significant genes (defined as an adjusted p-value <0.05) were generated from the centered rlog-normalized count data using pheatmap (version 1.0.12).

Gene pathway analysis using the Fischer and Monte-Carlo tests. To identify if a
pathway is relevant to MIS-C or acute-COVID-19 (mild and severe pediatrics patients), a comparison between MIS-C or acute-COV19 and the eight databanks described above was performed using the following steps below. To minimize false positive and artifactual results, all samples were processed using the same pipeline, Variant Explorer (VExP) (6), starting with their raw data (Fastq files).
Step 2 (Variant filtering): Variant analysis was performed in each family based on three filtering criteria: first, include variants predicted by ANNOVAR to have a potential functional coding consequence, including stop gain or loss, splice site disruption, indel, and nonsynonymous. Second, variants are filtered based on allele frequency in control populations (gnomAD, ExAC, EVS, 1000GP, and internal data from 8114 unaffected individuals from BCH). Heterozygous/hemizygous variants were included if minor allele frequency (MAF) was <0.0005 (0.05%) in any database. In comparison, homozygous variants were included only if MAF was ≤0.00005 (0.005%) and for compound heterozygous models the MAF cutoff was ≤0.01 (1%) with no homozygous variant reported in any database. The variants were further prioritized to include those with read depth ≥10X, alternative depth ≥5X, allele balance ≥0.20, and deleterious prediction (4 or more of 23 softwares, including PolyPhen, SIFT, FATHMM, and CADD).
Step 3  Heterozygous variants and c) Homozygous and/or 2 or more heterozygous variants in the same gene with a minimum distance between them of 100 base pairs (compound heterozygous filters). P values were calculated using 2 methods: traditional Fisher test (two sided) and Monte-Carlo method. The expectation for one event (pathway) using Monte Carlo method is described by the following formula: Where "F" represents the number of families with rare pathogenic variants in the "k" pathway (k=1:8,626 pathways) and "X" is a random control group with the same number of samples of the comparison group, for MIS-C, 39 samples and for acute-19, 24 samples.
Independent samples were taken random using a uniform distribution and 4682 samples described above. "N" is the total number of independent simulations (10,000 in total). The use of independent samples was very important to establish fairness in our tests, so then we use only one sample per family (probands).

Fig. S4. Characterization of circulating CD4 + Treg and Tconv cells in MIS-C. A,B.
Flow cytometric analysis and cell frequencies of T cell activation state markers (CD45RA,

Flow cytometric analysis and frequencies of IFNg and IL-17-expressing Treg (C) and
Tconv (D) cells of the respective subject groups. Each symbol represents one subject.        ▲All children with MIS-C and KD had an echocardiogram preformed.