Downregulation of exhausted cytotoxic T cells in gene expression networks of multisystem inflammatory syndrome in children

Abstract

Multisystem inflammatory syndrome in children (MIS-C) presents with fever, inflammation and pathology of multiple organs in individuals under 21 years of age in the weeks following severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection. Although an autoimmune pathogenesis has been proposed, the genes, pathways and cell types causal to this new disease remain unknown. Here we perform RNA sequencing of blood from patients with MIS-C and controls to find disease-associated genes clustered in a co-expression module annotated to CD56^dimCD57⁺ natural killer (NK) cells and exhausted CD8⁺ T cells. A similar transcriptome signature is replicated in an independent cohort of Kawasaki disease (KD), the related condition after which MIS-C was initially named. Probing a probabilistic causal network previously constructed from over 1,000 blood transcriptomes both validates the structure of this module and reveals nine key regulators, including TBX21, a central coordinator of exhausted CD8⁺ T cell differentiation. Together, this unbiased, transcriptome-wide survey implicates downregulation of NK cells and cytotoxic T cell exhaustion in the pathogenesis of MIS-C.

Fig. 1. Study overview.

Workflow of the analyses presented in the paper. A RNA-seq was generated on whole blood samples from MIS-C cases, pediatric COVID-19 cases, and healthy controls. B Deconvolution estimated the relative cell-type proportions in each transcriptome, which were compared between cases and controls to identify the immune cell types involved in the pathogenesis of MIS-C. C Expression of each gene in the transcriptome was tested for association with disease, resulting in a MIS-C signature that was queried to resolve the dysregulated molecular pathways. D Co-expression network construction organized genes into coherent units called modules. E and F Modules loaded with genes of the MIS-C signature were empirically identified (E), validated using DE signatures from a large Kawasaki disease (KD) gene-expression dataset (F), and functionally annotated to pathway, cell type, and other disease signatures (F). G The module with the strongest enrichment for MIS-C that also enriched for KD signatures was further annotated to pinpoint cell subtypes, and key regulators of the processes captured by this module were identified in a regulatory network built from whole blood gene expression in an independent cohort. This figure was created with BioRender.com.

Fig. 2. Differential expression analyses identify the transcriptional signature of MIS-C.

A Differential expression (DE) analysis for MIS-C patients versus HCs. The x-axis is the mean normalized count for each gene and the y-axis is the log₂(fold-change) for the differential expression. Positive and negative log₂(fold change) represent genes upregulated and downregulated in MIS-C, respectively, and the significance of association between gene expression and MIS-C status is indicated by the color of the dots as defined in the legend. B Overlap of MIS-C and pediatric COVID-19 transcriptional signatures: Venn diagram of the overlap of genes across DE signatures. Each comparison is labeled on the plot. C GO terms for MIS-C signature: GO term enrichment results for the top 10 upregulated and downregulated processes in MIS-C compared to HCs. Two-sided Fisher tests were used and p-values were adjusted for multi-testing as described in the “Methods” section. Full DE results and pathway enrichments for all comparisons in B can be found in Supplementary Data 7 and 8.

Fig. 3. Co-expression network analysis identifies modules of genes dysregulated in MIS-C.

A Module GO term enrichments. The y-axis is the most significant GO term associated with each module, the x-axis is the −log₁₀(adjusted p-value) for the enrichment (performed using Fisher’s test). Bars are colored by module names, which are also specified for clarity. Only modules with an enrichment p-value < 1 are shown. The threshold for significance, −log₁₀(0.05), is indicated by the red dashed line. In bold lettering are the modules enriched for the MIS-C signature. B Module cell type signature enrichments. The x-axis and y-axis are the names of the modules and of the cell type signatures, respectively. All signatures shown here are derived from the LM22 reference (Newman et al. Nat. Methods, 2015). The color and size of the circles represents the log₂(odds ratio) of the enrichments as defined in the legend. In bold are the modules enriched for the MIS-C signature. Only enriched cell types are shown. C Modules enrichment for MIS-C signatures. The x-axis is the module names and the y-axis the odds ratio of the enrichment of the modules for the genes upregulated and downregulated in MIS-C. Only modules significantly enriched for MIS-C DEGs are shown. The color of the bars represent direction and the opacity represent the significance as defined in the legend. All p-values were adjusted for multi-testing as described in the “Methods” section. All module enrichments can be found in Supplementary Data 7, 9–12, 21.

Fig. 4. Disease enrichments in MIS-C modules and cyan.

Slices are modules and each circular row represents the corresponding disease signature defined in the legend on the right. The purple outer rim of the plot represents the sum of the OR for all enrichments in that slice (performed using Fisher’s test). Each slice is divided into two components that show the OR for the enrichment of genes with ±log(fold change) in the corresponding disease signature (red and blue, respectively, as defined in legend). Numbers in the category slice map the circular rows to disease signatures. ORs are shown for disease signature enrichment adjusted p-values < 0.05. Diseases were grouped in biologically meaningful clusters as defined in the legend. Parenthesis in the legend refers to the source of the signature and the tissue or cell type assessed. References for disease signatures are defined in the “Methods” section and all enrichments can be found in Supplementary Data 15.

Fig. 5. Further dissection of skyblue implicates cytotoxic lymphocyte subtypes in MIS-C.

A Defining NK cell and CD8⁺ T cell subtypes in skyblue (Supplementary Data 19). The x-axis represents the OR for the enrichment of skyblue for specific cytotoxic cell subtypes and the y-axis the signatures used. Cell one and cell two refer to the cell subtypes being compared to generate the DE signatures projected onto skyblue and the direction of the bar is shown towards the upregulated cell type in the comparison. The color and length of bars represent the OR of the enrichment test as defined by the x-axis and the legends. For NK cells, we used signatures from Yang et al., Nat. Commun., 2019 (labeled NK cells), and from Collins et al. Cell 2019 (labeled CD56^dim NK subtypes). For CD8⁺ T cells, we used signatures from Wherry, et al. Immunity, 2007 (labeled CD8⁺ T cells), and from Beltra et al. Immunity, 2020 (labeled CD8⁺ Tex subtypes). Detailed descriptions of the subtypes can be found in the Methods. B–E Skyblue key drivers measures. The x-axes of the panels represent the key driver genes. B Key driver analysis results (hypergeometric test).The y-axis of this panel is the −log₁₀(Bonferroni adjusted p-value) of key driver analysis (Supplementary Data 21). The red dashed line is the significance threshold at −log₁₀(0.05). C Related disease differential expression. In this panel, the y-axis is the disease name corresponding to Fig. 4. The color represents the direction in the comparison as defined in the legend (defined per the associated reference). D Cell-type-specific expression of key driver genes (Supplementary Data 22). The y-axis is the cell type signatures from the following references: 1 = Newman et al. Nature Methods, 2015; 2 = Park et al. Science, 2020; 3 = Wilk et al. Nature Medicine, 2020; 4 = Liao et al. Nature Medicine, 2020; and 6 = Szabo et al. Nature Communications, 2019. E CD8⁺ T cell and NK cell subtype-specific signatures. The y-axis is the cell-type name corresponding to A and the color represents the direction in the comparison as defined in the legend (defined per the associated reference).