The raw scRNA-seq FASTQ files of PBMCs from 13 patients and 5 healthy controls were downloaded from the Genome Sequence Archive of the Beijing Institute of Genomics (BIG) Data Center (http://bigd.big.ac.cn/gsa-human, accession HRA000150). The 13 patients with COVID-19 were classified into three clinical conditions: moderate (n = 7), severe (n = 4) and convalescent (conv; n = 6, of whom 4 were paired with moderate cases). The downloaded reads were then processed individually using the Cell Ranger (v.4.0.0, 10xgenomics, https://www.10xgenomics.com/) count pipeline with the GRCh38 human reference genome to generate gene expression matrices. The subsequent analyses were performed by R (v.4.0.2) scripts with the Seurat (v.3.2.2) package as described in our previous study (13). Briefly, the Cell Ranger output filtered gene expression matrices were further filtered according to the parameters in our previous study (13). The genes expressed at a percentage of 0.1% cells or more were kept for every sample. The cells were filtered as following criteria: (1) the number of genes is equal to or larger than 500; (2) the number of unique molecular identifiers (UMIs) is equal to or larger than 800; (3) the percentage of the mitochondrial genes is no more than 10%. Those cells that did not satisfy the above criteria were removed. Then, the datasets from different samples from four conditions were integrated into an integrated and unbatched dataset using the “standard workflow”, as described at https://satijalab.org/seurat/v3.2/integration.html. The integrated dataset was scaled and principal components analysis (PCA) was calculated. The top 20 principal components (PCs) were selected to construct a shared nearest neighbor (SNN) network and an unsupervised graph-based clustering approach, the Louvain algorithm, was applied to cluster cells with a parameter resolution = 1.5. Clusters were then classified and annotated based on the selected classic markers (13). Specifically, MAIT cells were located using marker SLC4A10 (encoding solute carrier family 4 member 10) and TRAV1-2 (encoding T cell receptor alpha variable 1-2) and were extracted to another dataset to perform downstream analyses. Finally, uniform manifold approximation and projection (UMAP) was applied to visualize the clustering result in a two-dimensional space (27).

To identify the DEGs across different clusters and/or conditions, the “FindMarker” function in the Seurat package was performed with multiple threshold parameters, including an average log2 (fold change) ≥ 0.5, a Benjamini–Hochberg-corrected P value ≤ 0.01, and detection in ≥ 10% of cells in at least one condition. The obtained DEGs were uploaded to the Metascape webtool (www.metascape.org) (28) and the gene sets derived from gene ontology (GO) Biological Process Ontology were selected to obtain their function profiles.

We used cell scores to evaluate the degree to which individual cells express a certain pre-defined gene set (13, 29, 30). For a given gene set (G_j ) reflecting a specific cell state or biological function, the score for every cell i, SC_j(i), quantifying the relative expression of G_j in cell i as the average relative expression (Er) of genes in G_j compared to the average relative expression of a control gene set (G_j ^cont): SC_j(i) = average[Er(G_j, i)] - average[Er(G_j ^cont , i)]. The control gene set was defined by first binning all the analyzed genes into 25 bins of aggregate expression levels. For each gene in the given gene set, we randomly chose 100 genes from the same expression bin. The “AddModuleScore” function in Seurat package was used to implement this approach with “nbin” set to 25. We respectively used CELL KILLING (GO:0001906), CELL CHEMOTAXIS (GO:0060326), APOPTOTIC SIGNALING PATHWAY (GO:0097190), 9 activation-related genes (CD38, CD69, CD25(IL2RA), CD95(FAS), CD134(TNFRSF4), CD137(TNFRSF9), CD154(CD40LG), MKI67, KLRG1) to define the cell killing, cell chemotaxis, apoptosis, and the activation score.

The Cell Ranger (v.4.0.0) vdj pipeline with GRCh38 as the reference was used to perform the gene quantification and TCR clonotype assignment. Then, the files with contig annotations and clonotype frequencies were obtained. According to our previous study, only cells with at least one productive TCR α-chain (TRA) and one productive TCR β-chain (TRB) were kept for further analysis (13). The cell barcodes and sample identifiers were concatenated to define the cell identifiers and to associate the gene expression data with TCR information for each MAIT cell. Here, a clonotype refers a unique TRA(s)-TRB(s) pair (complementary determine region 3 (CDR3) amino acid sequences included). If the cell numbers of one clonotype were greater than 1, this clonotype was considered to be clonal, and the number indicated the degree of clonality of the clonotype.

MAIT cells were identified as CD161^hiTCR Vα7.2⁺ cells among CD3⁺ T cells. To detect surface markers, the following antibodies were used: anti-CD3-APC-Cy7 (BD Biosciences, San Diego, California, USA. Clone: OKT3), anti-CD161-PE (BD Biosciences. Clone: HP-3G10), and anti-TCR Vα7.2-BV421 (BD Biosciences. Clone: 3C10). For fluorochrome-labeled inhibitors of caspases (FLICA) caspase-1 detection, FLICA staining was conducted in accordance with the manufacturer’s instructions (Bio-Rad, Hercules, CA, USA). Cells were incubated with the FLICA caspase-1 reagent for 1 hour at 37 C and then cells were washed for downstream surface markers staining and intracellular anti-capase-3-Alexa 647 (BD Biosciences. Clone: C92-605). Data were acquired on a BD FACSCanto II flow cytometer (BD Biosciences), and further analyzed using FlowJo (Ashland, OR, USA) Tree Star software.

All the boxplots in this study were plotted using “geom_boxplot” in ggplot2 (ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016) R package. Every dot represented a sample. The horizontal line within each box acted as the median, and the bottom and top of each box indicated the 25-th and 75-th percentile. The Wilcoxon rank-sum test was applied to test the significance of the difference between conditions using the “geom_signif” function in ggsignif (Significance Brackets for ‘ggplot2’) R package. The statistical tools and methods for each analysis are explicitly described in the figure legends.

To characterize the transcriptional profiling of MAIT cells in patients with COVID-19, we analyzed the published scRNA-seq datasets of PBMCs from thirteen patients and five healthy donors (HDs) (13) ( Figures 1A, C ). The 13 patients with COVID−19 consisted of moderate cases (n = 7), severe cases (n = 4), and convalescent (n=6; of whom 4 were paired with the moderate cases) ( Figure 1C ). Fourteen cell types comprising a total of 121464 cells were selected and annotated after quality control and doublets removal. Specifically, according to the expression level of canonical markers SLC4A10 and TRAV1-2, a total of 2502 MAIT cells were identified and extracted for subsequent detailed analyses ( Figures 1B and S1A ). Among these MAIT cells, 417 cells (16.7%) were from HDs, 1019 cells (40.7%) were from moderate cases, 130 cells (5.2%) were from severe cases, and 936 cells (37.4%) were from convalescent samples. Meanwhile, there were an average of 1.97%, 2.58%, 0.53%, and 2.53% MAIT cells among PBMCs per sample in the HD, moderate, severe, and conv groups, suggesting profound depletion of MAIT cells in patients with severe disease. The changed MAIT proportions in the different conditions agreed with the previous findings (13, 25, 26). Visualization of MAIT cells via UMAP clearly demonstrated their batch-less nature and comparability ( Figure 1C ). As such, we clearly defined the MAIT cells in the peripheral blood of individuals with COVID-19.

To investigate the transcriptomic features of MAIT cells comprehensively, we first compared the expression patterns of MAIT cells with 13 other immune cell types, including naïve T cell, activated T cells, γδ T cells, proliferative T cells, natural killer cells, B cells, plasma B cells, CD14+ monocytes, CD16+ monocytes, monocyte−derived dendritic cells, plasmacytoid dendritic cells, platelets, and hemopoietic stem cells, in PBMCs and identified a total of 230 differentially expressed genes (DEGs). Further analysis found that these DEGs were mainly involved in myeloid leukocyte activation, lymphocyte activation and migration, antigen processing and presentation, and cytokine production ( Figure 1D ), suggesting that MAIT cells may engage in the immune response against SARS-CoV-2. Second, we compared the expression patterns of MAIT cells from patients with moderate or severe disease with the counterparts from HDs. We found that DEGs were mostly enriched in the type I interferon signaling pathway, response to interferon (gamma and/or beta), regulation of innate immune response, positive regulation of cytokine production, and NF-kappaB transcription factor activity (TFA) processes ( Figures 2A, B ). These results agreed with the previously observed results in 13 other immune cell types from patients with COVID-19 (13, 31). Third, we compared the transcriptional profiles of MAIT cells from patients with severe disease with the counterparts from patients with moderate disease to explore the correlation of disease severity and gene expression in MAIT cells. A total of twelve DEGs were recognized. Ten of them (TRBV9, TRAV8-2, S100A8, GZMH, S100A9, KLF6, CD8B, KLRD1, IGLV3-19, and JCHAIN) were upregulated, whereas the only two were downregulated genes (SLC4A10 and TRAV1-2), which are canonical markers of MAIT cells ( Figure 2C ). Further analysis found that the DEGs were mainly enriched in adaptive immune response and humoral immune response ( Figure 2D ), suggesting that the expression of these DEGs might correlate with the disease severity. Interestingly, we found that 18 DEGs were in at least two conditions ( Figure 2E ). Specifically, twelve of them (IFITM1, IFI44L, IFI6, ISG15, XAF1, LY6E, MX1, IRF7, OAS1, EIF2AK2, TRIM22, and TXNIP) showed increased expression during COVID−19 and later declined in the convalescent condition; however, six DEGs (RGCC, LMNA, ZFP36, MT-ND6, JUN, and FOS) showed the opposite; i.e., the expression levels of these genes were reduced during infection with COVID-19 and recovered after being restored to health ( Figure 2F and Table S1 ). These findings suggested that the encoded proteins might be involved in the virus-host interaction and host immune response, thus reflecting disease severity.

To further investigate the features of MAIT cells in patients with COVID-19, we used a scoring system to evaluate certain key biological processes, such as activation, cell killing, cell chemotaxis, and apoptosis. Optimal cell activation is a prerequisite for cell function, and we found that MAIT cells in patients with COVID-19 showed trends if higher activation levels than in HDs, suggesting a consistent response by MAIT cells to SARS-CoV-2 infection. Notably, the activation levels of MAIT cells in severe cases were less than those in moderate and conv cases. For functional evaluation, we found that the cell killing level in MAIT cells increased consistently with disease severity and displayed the highest levels in severe cases. Importantly, cell chemotaxis and apoptosis levels in MAIT cells were also highest in severe cases, suggesting that the increased migration and cell death pathways in the MAIT cells of patients with severe COVID-19 may be associated with their MAIT-cell loss in peripheral blood ( Figure 3A ).

To investigate the association between pyroptosis, a newly identified programmed cell death pathway, and disease severity during SARS-CoV-2 infection, we first compared CASP1 gene expression and found the CASP1 gene expression in patients with COVID-19 showed higher levels than those in in HDs ( Figures 3B and S2 ). To confirm this finding, we used FLICA to detect the percentages of the active form of caspase-1⁺ cells. We observed that the frequencies of FLICA caspase-1⁺ MAIT cells were significantly increased in the peripheral blood of patients with COVID-19. By contrast, the frequencies of active caspase-3⁺ MAIT cells was not significantly increased in patients with COVID-19 ( Figures 3C and S3 ). Interestingly, CD8 T cells, but not CD4 T cells, also showed an increased pyroptotic phenotype in patients with COVID-19 ( Figures 3C and S3 ). Collectively, these data, at least in part, confirm that MAIT cell pyroptosis occurs in SARS-CoV-2 infection and pyroptosis as possible cause of MAIT cell loss.

To better determine the clonal relationship among individual MAIT cells and their usage of V(D)J genes across the four conditions, we reconstructed TCR sequences of MAIT cells from the TCR sequencing of 20 cases. In this analysis, 1687 cells matching TCR information were detected from 2318 MAIT cells and the percentage of cells with TCR information is more than 70% in all conditions, except for the HD ( Figures 4A, B ). We then performed statistical processing of all clonotypes ( Figures 4C, D ). Overall, the number of clonotypes was negatively logarithmically correlated with the number of cells per clonotype ( Figure 4C ). It should be noted that the number of high frequencies of clonotypes of MAIT cells was relatively high. Compared with the HDs, clonal expansion was obvious in patients with COVID-19 both in disease progression and in convalescence ( Figure 4D ). Interestingly, there were more large clonal expansions (clonal size >10) in the severe cases than in the other conditions, indicating that some clonotypes might be abnormally over expanded in the MAIT cells of severe cases. To study the abnormality and gene preference of TCRs in MAIT cells of severe cases, we compared the usage of V(D)J genes across the four conditions ( Figures 4E–G ). The top 10 complementarity determining region 3 (CDR3) sequences in TRB were quite different across the four conditions ( Figure 4E ) and the moderate and conv conditions shared three CDR3 sequences because four samples from these conditions were paired. However, for all the top 10 CDR3 sequences across the four conditions, there were only 21 unique CDR3 sequences in TRA. In total, the usage percentage of the top 10 CDR3 sequences in the HD condition was lower and more balanced compared with those of the other three conditions. Notably, we found that the percentage of cells using TRAV1−2 or TRBV33 was much higher than that of other genes in TRA and the percentage in severe cases was significantly decreased compared with that in the moderate and conv cases ( Figure 4E ). TRAV1-2 is a marker gene of MAIT cells and it was also expressed at a lower level inside MAIT cells of severe cases compared with that in the other conditions. We also found that the usage levels of the gene pairs TRAV1-2 and TRAJ33, and TRAV1-2 and TRAJ20 were high across all conditions ( Figure 4G ). Besides, the usage level of gene pair TRAV1-2 and TRAJ12 was also high across all conditions except for that in the severe cases. Thus, it seems that in TRA, these three gene pairs were mainly used only in moderate and conv cases. The lost gene pair TRAV1-2 and TRAJ12 might influence MAIT cell antigen recognition of patients with severe disease.