Main

Although SARS-CoV-2 mainly replicates in the respiratory tract, and lower respiratory tract complications are major drivers of morbidity and mortality, it is unclear to what extent immunological dynamics observed in blood can be translated to the respiratory tract. Nasopharyngeal and peripheral blood cytokines and antibodies do not correlate during COVID-19 (ref. 7). While the mucosal immune system remains understudied, several studies demonstrated increased neutrophil levels and activated alveolar macrophages/monocytes in bronchoalveolar lavage (BAL)8 or nasopharyngeal and oropharyngeal swabs9 from individuals hospitalized with COVID-19. T cell recruitment to the respiratory tract might be beneficial, as increased BAL CD8+ T cell clones8 or increased tracheal aspirate CD4+ T cells associate with less severe disease or survival10. Nasopharyngeal swabs have been mostly used to analyze upper respiratory tract responses9,11,12 but predominantly sample superficial cells, such as epithelial cells and neutrophils/monocytes, while they incompletely capture other immune cells, such as T cells. Although BAL and tracheal aspirates provide good insight into lower airways8,10,13,14, they are difficult to collect longitudinally, after recovery or from healthy individuals and individuals that do not require intubation. Therefore, we still have a limited understanding of how COVID-19 affects mucosal immunity.

Here, we performed a prospective observational cohort study to characterize mucosal immune cell dynamics in the upper respiratory tract during acute, early recovery and convalescence stages of COVID-19. Individuals with PCR-confirmed SARS-CoV-2 infection were recruited after hospital admission (Fig. 1a). Longitudinal nasal curettage samples from 20 hospitalized individuals were collected, with up to 4 samples analyzed per individual. Curettage also samples immune cells from the nasal mucosa that are absent from the lumen15. Samples were collected from 2 to 61 d after hospital admission (11–82 d after symptom onset). We stratified hospitalized individuals into those with acute infection (n = 9, 2–11 d since hospital admission) or in early recovery stage (ERS), which is defined by having moved from the intensive care unit (ICU) to ward (n = 11, 15–61 d since hospital admission and with an ICU stay period of 4–55 d). Individual characteristics, comorbidities, outcome and treatment are shown in Supplementary Table 1. Sixteen individuals were also sampled 5–6 weeks after hospital discharge (convalescent individuals with a median of 77 d after symptom onset). Twelve sex- and age-matched healthy individuals with negative SARS-CoV-2 IgG and without symptoms of a respiratory tract infection were included.

Fig. 1: Cytometry by time of flight (CyTOF) analysis of nasal immune cells during and after COVID-19 infection.
figure 1

a, Participant timelines. For each included individual, symptom onset (purple cross), hospitalization (black bar) and ICU stay (red bar) are indicated, aligned to the day of hospital admission. The orange x indicates one individual who was included but later transferred to a different hospital. One individual was discharged and then readmitted 1 d later. Blue diamonds indicate nasal curettage samples. The asterisk indicates one individual who was hospitalized for unrelated reasons at the time of positive test and symptom onset, and hospital admission was set at day of symptom onset. b, Hierarchical stochastic neighbor embedding (SNE) was used to cluster cellular landmarks on 37 markers into 12 populations. Some of these populations were then further divided into subpopulations in a second t-distributed stochastic neighbor embedding (t-SNE) plot embedding at a data level, as indicated in the rectangles. c, Heat map of marker expression per population. Median intensity per population is shown after arcsin transformation; ILCs, innate lymphoid cells; mDCs, myeloid dendritic cells; pDCs, plasmacytoid dendritic cells; TEM, effector memory T cells; DP, double-positive T cells (CD4+CD8+); DN, double-negative T cells (CD4CD8); EMRA, effector memory reexpressing CD45RA; Treg, regulatory T cells; gran, granulocytes; mono, monocytes.

In total, 875,564 nasal CD45+ immune cells were analyzed from 56 samples (44 samples from 29 individuals with COVID-19 and 12 samples from healthy donors) using a 39-marker mass cytometry panel (Supplementary Table 2). Nasal CD45+ immune cells were divided into eight main lineages and further subclustered into 28 populations (Fig. 1b,c). At a lineage level, granulocytes dominated the nasal immune profile of individuals with acute COVID-19 (Fig. 2a). This progressively decreased as granulocyte frequencies were slightly lower in ERS and further reduced during convalescence, more similar to healthy age-matched donors. All other lineages, apart from monocytes, appeared decreased during acute infection compared to healthy donors. To understand whether this apparent depletion was related to increased granulocytes and monocytes or to true lymphopenia, we normalized immune cell numbers to epithelial cell numbers for each sample, permitting an independent assessment of immune cell populations while correcting for variable sample yield. In healthy donors and recovered individuals, there was a strong correlation between epithelial and immune cell yields, as expected (Fig. 2b). This association was absent in hospitalized individuals, suggesting that the vastly increased granulocyte numbers are caused by a strong influx into the nasal mucosa. Indeed, when normalizing to epithelial cells, granulocytes and monocytes were highly increased during acute infection and ERS, while a non-significant trend for increased granulocytes was shown during convalescence compared to healthy donors (Fig. 2c and Supplementary Table 3). Other main cell lineages (B cells, NK cells, monocytes, pDCs, mDCs, CD4+ T cells and CD8+ T cells) were not statistically different between individuals with acute infection and healthy individuals. The observation that nasal lymphocyte numbers were unchanged is in contrast with the general lymphopenia detected in peripheral blood2,6.

Fig. 2: Nasal cell lineage abundance during and after COVID-19.
figure 2

a, Stacked bar charts showing the composition of the nasal immune system in acute COVID-19 (red), during ERS (after ICU but still in the hospital; orange) or in COVID-19 convalescence (5–6 weeks after hospital discharge; pink) or in healthy donors (HD; blue). b, Correlation analyses between nasal immune and epithelial cells for hospitalized individuals (left) and convalescent individuals and healthy donors (right) are shown. Ranks of individuals are shown with color corresponding to group as well as a regression line (blue), 95% confidence intervals (shaded area) and results from Spearman correlation analysis. c, Ratio of nasal immune cell types normalized to the number of epithelial cells from the same sample. Individuals and box plots are shown, and paired samples between are indicated by gray lines. If a cell type was not detected in at least one sample, half the value of the lowest recorded number was added to each sample before log transformation. d, Heat map showing log10 relative abundance (RA) of nasal cell clusters scaled per cluster. e, Box plots of nasal clusters during acute COVID-19. Box plots depict median and interquartile ranges, with whiskers extending to 1.5× interquartile range or maximum value; *P < 0.05 and ***P < 0.001 by linear mixed model with group as fixed effect and individuals as random effect with post hoc testing and Tukey multiple testing correction, followed by Benjamini–Hochberg correction for comparing multiple lineages or subsets. See Supplementary Tables 3 and 4 for exact test results. Samples from acute individuals (n = 9 individuals, 10 independent samples), ERS individuals (n = 11 individuals, 18 independent samples), convalescent individuals (n = 16 individuals) and healthy donors (n = 12 individuals) are used; Neutro, neutrophils. Only the first sample per donor in a time point (acute or ERS) is shown for a, c and e, but all are included in statistical modeling.

On the subclustering level (Fig. 1b,c), 8 of the 28 defined cell clusters were significantly elevated during COVID-19 (Fig. 2d,e and Supplementary Tables 3 and 4). These showed clear associations with time from hospital admission (Extended Data Fig. 1). Three monocyte/macrophage populations, defined by CD163/CD206 expression, were increased during acute infection and ERS compared to convalescence and healthy donors. The granulocyte subsets CD16hi, CD16dim and CD16 neutrophils were all elevated to a variable degree during acute stage and in ERS compared to healthy donors. More detailed flow cytometric analysis (Extended Data Fig. 2) showed that nasal CD16 granulocytes expressed CXCR4 but not Siglec8 and were therefore not eosinophils but neutrophils that had shed CD16. Expressing CXCR4 and losing CD16 is typical for aged neutrophils that are (pre)apoptotic16,17. Furthermore, the expression of CD10 indicated that they were fully maturated18. Although there was no overall change in CD4+ T cell numbers, effector (CCR7CD45RO+) CD4+ T cells were 18-fold increased during acute stage compared to healthy donors. CD8+ EMRA T cells (CCR7CD45RA+) showed a non-significant trend toward increased numbers during acute infection. This agrees with earlier reports on peripheral blood T cells showing higher induction of antigen-specific CD4+ T cells than of CD8+ T cells, while the majority of SARS-CoV-2-specific CD8+ T cells are EMRA19,20. These short-lived effector cells returned during ERS and convalescence to similar numbers as observed in healthy donors. Finally, CD11c+ NK cells were increased in hospitalized individuals (acute and ERS) compared to convalescent individuals and healthy donors. CD11c marks NK cells with increased interferon-producing capacity and effector function21. Thus, a dynamic recruitment of various adaptive and innate populations mediating inflammation and antiviral function to the upper respiratory tract was observed during hospitalization that normalized in convalescence to levels resembling those in healthy donors. Of note, nasal B cells were not increased in hospitalized individuals nor did we detect plasmablasts in nasal mucosa, corroborating observations that mucosal antibody levels are reduced compared to systemic titers in hospitalized individuals7. Alternatively, B cells might not migrate past the submucosa, and further studies would be required to address mucosal B cell responses early in infection and/or in individuals with mild infection. To understand whether factors like sex, comorbidities and medication were drivers of nasal immune profiles, we performed multidimensional scaling (MDS) using all cell subsets (Extended Data Fig. 3). Acute individuals clustered separately from healthy donors and convalescent individuals, with ERS individuals intermediate. There was no clustering based on any of these covariates, although larger sample sizes are needed to conclusively exclude such effects.

Next, we analyzed how nasal immune populations associated with systemic immune cells. We measured paired cryopreserved whole-blood samples using CyTOF (n = 45) and obtained absolute cell counts from fresh whole blood using flow cytometry (n = 43) (Extended Data Fig. 4). With CyTOF, we identified 133 clusters and 15 subsets in blood (Extended Data Fig. 5), showing good concordance for most subsets with the freshly acquired flow cytometry data (Extended Data Fig. 6). Of all CyTOF clusters and subsets, 41 were significantly different between the groups after correction for multiple testing (Fig. 3a and Supplementary Tables 5 and 6). Acute COVID-19 was associated with a general lymphopenia, with reduced numbers of naive and central memory CD4+ T (TCM) and CD8+ T cells, mucosal-associated invariant T (MAIT) cells, certain CD4+ Teff cell clusters, pDCs, T cell antigen receptor-γδ (TCRγδ) T cells, NK cells and ILCs. By contrast, HLA-DRlow and HLA-DRmed monocytes, CD163+ monocytes, CD11c+ neutrophils and two clusters of CD38+PD1+CD4+ Teff cells were increased during acute COVID-19. In convalescence, these perturbations were partially reverted, as only pDCs and NK cells remained significantly lower than healthy donors after multiple testing correction (Supplementary Table 6 and Extended Data Fig. 7). Flow cytometry confirmed in absolute counts the depletion of pDCs, NK cells and CD4+ TCM cells (Extended Data Fig. 7). In addition, it showed that eosinophils and basophils were significantly depleted from blood during acute stage, while plasmablasts were increased (Fig. 3b). Our findings in blood thus agree with previous reports1,2,3,4,5.

Fig. 3: Integration of blood and nasal immune responses during and following COVID-19.
figure 3

a, Heat map showing significant CyTOF whole-blood (WB) clusters. Individuals are shown in columns (n = 34 individuals, 45 samples) ordered by days since hospital admission per group. b, Absolute counts of cellular subsets in blood. Acute (n = 4 individuals, 5 samples), ERS (n = 7 individuals, 10 samples), convalescent (n = 16 individuals) and healthy donors (n = 12 individuals). Individuals and box plots are shown, and paired samples are indicated by gray lines. Box plots depict median and interquartile ranges, with whiskers extending to 1.5× interquartile range or maximum value. Only the first sample per donor per time point is shown. Results were compared against healthy donors from testing using linear mixed models with group as fixed effect and individuals as random effect with post hoc testing and Tukey multiple testing correction. c, MDS of significant nasal mucosa (purple) or whole-blood (red) CyTOF clusters, whole-blood flow subsets (green) and clinical parameters (cyan). Only samples from hospitalized individuals with complete data were included (n = 3 acute, n = 10 ERS with 15 samples). d, Correlation plots between nasal and whole-blood CyTOF subsets. Individuals with paired data (n = 34 individuals, 45 samples) are shown, with colors indicating disease status. Estimates and 95% confidence intervals from linear regression analyses are shown as black lines with shaded areas. Pearson correlation analysis results are depicted. e, Multiomics factor analysis (MOFA) of all three datasets. Explained variation of all data explained per dataset or per latent factor (LF) is shown in the bar chart and heat map, respectively. f, Factor scores of LF1 and LF3 plotted against day of hospital admission. Healthy donors are separated by a dashed line. Loess regression estimate is indicated by the black line, and shaded areas represent 95% confidence intervals. g,h, The top five loadings per dominant dataset are indicated for LF1 (g) and LF3 (h). td, terminally differentiated. i, Association between LF3 and oxygen saturation (spO2) in blood. Symbols indicate hospitalized individuals (n = 20 individuals, 28 samples), and the dark line depicts the estimate from a linear mixed model with viral load, days since admission and spO2 included as fixed effects and individuals as random effects. The model estimate (est) and P value are depicted.

To understand how the perturbations in blood and nasal mucosa were linked to each other and clinical variables, we performed MDS using all significant populations (Fig. 3c). Proinflammatory nasal subsets generally clustered together with viral load. Indeed, viral load positively correlated with nasal monocyte and neutrophil levels (Extended Data Fig. 8). Nasal proinflammatory populations were also associated with blood lymphopenia (Fig. 3c and Extended Data Fig. 8). We next evaluated whether mucosal immune cells correlated with their systemic counterparts to assess whether peripheral responses can serve as correlates for mucosal events. Both MDS and hierarchical clustering of the correlation matrix suggested that increased blood monocyte clusters (HLA-DRlow and CD163+) associate with nasal monocytes/neutrophils and viral load. Indeed, HLA-DRlow monocyte numbers in blood were significantly correlated with nasal monocyte numbers. Moreover, CD38+PD1+CD4+ Teff cell clusters 14 (HLA-DR+) and 2 (CD27+) were strongly linked to nasal CD4+ Teff cells (Fig. 3d). Further characterization of nasal CD4+ TEM cells revealed two global subsets, with a group of activated cells expressing CD38, PD1 and CTLA-4 (Extended Data Fig. 9a–c). These activated cells were mostly Tbet+, suggesting a type 1 helper T cell (TH1) profile. Thus, we found that blood lymphopenia did not associate with the concurrent recruitment of lymphocytes to nasal mucosa, while it was associated with increased proinflammatory populations in nasal mucosa. In addition, we also found populations in blood, that is, CD38+PD1+CD4+ Teff cells and HLA-DRlow monocytes, that correlated with their nasal counterparts.

Nasal monocyte and granulocyte subsets, blood lymphopenia and perturbed clusters, viral load and time since hospital admission were all associated with each other. To understand which factors associated with clinical features and to disentangle these relationships, we integrated the immunological datasets using MOFA (Fig. 3e). This tool detects hidden LFs explaining variation both within and across different datasets22. Here, six LFs were found to which both nasal data and flow cytometric absolute blood counts contributed most. LF1 was strongly associated with time since hospital admission, increasing progressively during recovery and convalescence to levels similar to healthy donors (Fig. 3f). LF1 also was negatively correlated with viral load using marginal correlation, but this was not significant when correcting for time since admission using linear mixed modeling. The nasal monocyte subsets (CD163+ and classical monocytes), blood plasmablasts and CD4+:CD8+ T cell ratio were the most important immune populations negatively driving LF1 (and thus associated with acute infection) (Fig. 3g). Indeed, blood plasmablasts and nasal CD163+ monocytes were positively correlated (Fig. 3d). This links to previous findings showing that antibody development associates with increased nasal and blood inflammatory cytokines during COVID-19 (ref. 7). Immune populations positively contributing to LF1 (thus associated with longer time since infection, resolving of infection and immune homeostasis in healthy donors) were nasal CD8+ T cells and blood lymphocytes (TCRγδ T cells, CD8+ T terminally differentiated effector cells and NK cells). LF3 was strongly driven by nasal neutrophils, with CD16hi neutrophils as top loading factor (Fig. 3h). LF3 was negatively associated with blood oxygen saturation, considering time since hospital admission and viral load in linear mixed modeling (Fig. 3i). A similar linear model including CD16hi neutrophils instead of LF3 confirmed that these cells were independently and inversely associated with oxygen saturation, suggesting that independent of time since admission and viral load, enhanced nasal neutrophils are associated with hypoxia.

Subsequently, we looked more closely at phenotypic expression profiles of abundant nasal cell clusters and to what extent these profiles normalized after hospital discharge. Among all increased monocyte subsets, acute individuals had relatively more CD163+ and fewer CD163+CD206+ monocytes/macrophages than healthy donors, which normalized during recovery (Fig. 4a). These CD206+ cells are likely fully differentiated tissue-resident macrophages, while CD163+ monocytes could represent recently recruited monocytes in individuals with COVID-19 (refs. 10,23). This hypothesis was supported by trajectory analysis (Fig. 4b). HLA-DR expression was reduced on nasal monocytes/macrophages during hospitalization, which normalized during convalescence (Fig. 4c). This agrees with our finding that HLA-DRlow monocytes, but not classical HLA-DRhigh monocytes, in blood strongly correlate with nasal monocyte numbers. Low HLA-DR expression is typical for myeloid-derived suppressor cells (MDSC), and MDSC-like cell expansion in peripheral blood during severe COVID-19 was previously reported24. Our results suggest that these cells may rapidly seed the upper airway mucosa where they might further differentiate into macrophages. CD163+CD206+ monocytes/macrophages also expressed elevated interleukin-3 (IL-3) receptor (CD123) levels during hospitalization, while CD163+ and CD163 monocyte subsets generally lacked CD123 (Fig. 1c). We then investigated more closely CD16 neutrophils by further grouping them into seven subclusters (Fig. 4d). Subcluster 1, characterized by increased IL-7 receptor (CD127) expression, was significantly increased during hospitalization and convalescence compared to healthy donors (Fig. 4e,f). Thus, while total neutrophil numbers normalized during convalescence compared to healthy donors, alterations in their phenotype remained visible. The half-life of granulocytes is relatively low (hours to days)25, suggesting either ongoing recruitment of altered cells or continued local perturbation in the respiratory mucosa.

Fig. 4: Phenotypic changes of innate and adaptive cell subsets during and after COVID-19.
figure 4

a, Percentages of monocytes expressing CD163, CD163 and CD206 or neither (classical). Subset frequencies are shown for individuals with acute COVID-19 (red), during ERS (after ICU but still in the hospital; orange), in convalescent individuals (Conv; pink) and in healthy donors (blue). Box plots and individual data points are depicted, with paired individual samples connected by a gray line. b, Diffusion map showing trajectory analysis of monocyte/macrophage subsets. All cells are colored according to their clustered phenotype. c, Violin plots showing expression of HLA-DR and CD123 on CD163+CD206+ monocytes/macrophages. All cells are shown; ***P = 4.9 × 10–5, **P = 0.002 and *P = 0.032 by statistical modeling as described below on the mean signal intensity per sample. d, Clustering of CD16 granulocytes using t-SNE and Gaussian mean shift. e, Proportion of cluster 1 within CD16 granulocytes. Box plots and individual data points are depicted, with paired individual samples connected by a gray line. f, Violin plots showing expression of CD127 per cluster. Box plots depict median and interquartile ranges, with whiskers extending to 1.5× interquartile range or maximum value. Samples from acute individuals (n = 9 individuals, 10 independent samples), ERS (n = 11 individuals, 18 independent samples), convalescence (n = 16 individuals) and healthy donors (n = 12 individuals) are used. Only the first sample per donor in a time point (acute or ERS) is shown for a and e, but all are included in the statistical modeling. Statistical results from a linear mixed model with group as fixed effect and individuals as random effect with post hoc testing and Tukey multiple testing correction are shown.

In mice, influenza-specific CD8+ TRM cells persist in the nasal mucosa following infection and efficiently control secondary infections26. Therefore, we assessed whether long-term protective CD8+ T cell immunity develops in the nasal mucosa of individuals with COVID-19 to act as ‘gatekeepers’ and protect against reinfection. The majority of nasal CD8+ T cells highly expressed CD69 and very little KLRG1, defining them as TRM cells27 (Fig. 1b,c). Subclustering of CD8+ TRM cells showed variable expression of activation markers CD38, HLA-DR and Tbet (Fig. 5a). CD8+ TRM cells from acute phase, ERS, convalescent and healthy individuals clustered differentially, indicating altered phenotypes during and following COVID-19 (Fig. 5b). Indeed subcluster frequencies significantly differed between groups (Fig. 5c,d). Subcluster 5 (HLA-DR/Tbet/CD38 expression) was increased in hospitalized individuals, while subcluster 3 (HLA-DR/CD38 expression) was higher during convalescence than in healthy donors as were all CD38+ TRM cells. Thus CD8+ TRM cells have an increased activation profile, which persisted at least 5–6 weeks after hospital discharge.

Fig. 5: SARS-CoV-2-specific T cells are present in the nasal mucosa in convalescent individuals with COVID-19.
figure 5

a, t-SNE analysis of all 7,826 CD8+ TRM cells. Expression of CD38, HLA-DR or Tbet overlaid onto a t-SNE plot. Bottom right, clustering of t-SNE using Gaussian mean shift. b, Two-dimensional kernel density estimation from t-SNE plot of all cells divided per group. c, Heat map showing median expression for all markers per TRM cluster. d, Percentage of CD8+ TRM cells belonging to either cluster 3 or cluster 5 or all clusters expressing CD38. Box plots and individual data points are depicted, with paired individual samples connected by a gray line. Box plots depict median and interquartile ranges, with whiskers extending to 1.5× interquartile range or maximum value. Only samples with >10 CD8+ TRM cells are included, and samples from acute individuals (red; n = 4 individuals, 5 samples), ERS individuals (orange; n = 7 individuals, 13 samples), convalescent individuals (pink; n = 14 individuals) and healthy donors (blue; n = 12 individuals) are used. Only the first sample per donor per time point (acute or ERS) is shown, but all are included in the statistical modeling. P values by linear mixed model with group as fixed effect and individual as random effect followed by post hoc testing and Tukey multiple testing correction are shown. e, Bar plots showing the frequency of TCR clonotypes in nasal samples collected from four donors. T cell clones with non-paired TCRα or TCRβ sequences are depicted on the left and right side of the y axis, respectively, and are ranked by their frequency. Only clones with a frequency >1% among reads are shown for TCRβ clones, while all TCRα clones are shown. Red- and yellow-colored bars indicate clones also present in SARS-CoV-2-specific CD8+ T cells sorted by fluorescence-activated cell sorting (FACS) from paired peripheral blood. Venn diagrams depict the total number of clones from nasal samples (gray), SARS-CoV-2-specific CD8+ T cells in peripheral blood mononuclear cells (PBMC) based on activation after peptide stimulation (red) or tetramer binding (yellow), including clones accounting for <1% of total reads. Star indicates that this clone was found in both the tetramer-sorted and stimulation-sorted fraction. f, Bar plot showing the average nasal CD8+ T cell composition from the convalescent samples of four individuals with TCR sequencing. Naive and TEM cells were too infrequent to be visible and are not labeled.

To demonstrate antigen specificity, we performed TCR sequencing of nasal samples during convalescence. These sequences were compared with SARS-CoV-2-specific CD8+ and CD4+ T cells sorted by FACS from paired peripheral blood based on either induction of activation markers CD137 or CD154 after stimulation with SARS-CoV-2 peptide pools or following binding to SARS-CoV-2-specific peptide–major histocompatibility complex (pMHC) tetramer complexes (Extended Data Fig. 9d,e). For four convalescent individuals (two with and two without a prior ICU stay), we obtained >10 unique TCRs from both nasal cells and sorted SARS-CoV-2-specific peripheral blood cells. In all individuals, nasal TCRs overlapped with TCRs from FACS-sorted SARS-CoV-2-specific CD8+ T cells from peripheral blood (Fig. 5e and Supplementary Table 7). The number of paired T cell clones ranged from 1 to 12 per individual and included both α- and β-chains that overlapped with CD8+ T cells sorted based on activation following SARS-CoV-2 stimulation and/or tetramer binding. Two individuals had strongly immunodominant nasal SARS-CoV-2-specific TCRs, with each clone accounting for >10% of all TCR reads. Samples were collected 36–70 d after viral clearance, and, on average, 60.7% of the nasal CD8+ T cells in these convalescent samples were of a TRM phenotype (Fig. 5f), indicating that antigen-specific tissue-resident memory was induced. Of note, the number of unique SARS-CoV-2-specific T cell clones detected and their overlap between nasal mucosa and peripheral blood might be underestimated as (1) a limited number of pMHC tetramers could be included, and (2) in the activation-induced assay, blood CD8+ T cells were isolated based on reactivity toward structural proteins, while CD8+ T cell reactivity is also directed against non-structural proteins in individuals with severe COVID-19 (refs. 28,29,30). Indeed, matching nasal TCR sequences to publicly known SARS-CoV-2-specific T cell clones identified another specific T cell clone (Supplementary Table 7). In contrast to CD8+ T cells, there were no overlapping nasal TCRs and SARS-CoV-2-specific CD4+ T cells from peripheral blood (Supplementary Table 7), which may reflect a better induction or maintenance of nasal CD8+ T cells. Alternatively, it may be that too few nasal cells were obtained during sampling, missing nasal SARS-CoV-2-specific CD4+ T cells during convalescence as CD8+ T cells are more abundant than CD4+ T cells in the nasal mucosa. Taken together, we demonstrated that SARS-CoV-2-specific CD8+ T cells in the nasal mucosa can persist for months after viral clearance. This suggests the establishment of local protective immune memory responses that could rapidly control and attenuate reinfections by SARS-CoV-2.

In conclusion, we provide an in-depth analysis of how COVID-19 affects nasal mucosal immunity during acute infection, early recovery and convalescence. We confirmed that acute COVID-19 leads to transient increases of respiratory granulocyte/monocyte subsets and extended this to CD11c+ NK cells and CD4+ Teff cells. This proinflammatory response positively correlated with viral load and inversely correlated with time since admission, while specifically nasal neutrophils showed a strong negative association with oxygen saturation in blood. Furthermore, increased blood HLA-DRlow monocyte numbers correlated with their emergence at the nasal mucosa, and blood CD38+PD1+CD4+ Teff cells positively associated with nasal CD4+ Teff abundance. Strikingly, unlike in blood, no general lymphopenia was observed in the nasal mucosa. Simultaneously, lymphocyte populations depleted from blood did not show a concurrent increase in nasal mucosa, suggesting that lymphopenia is not caused by nasal mucosal homing, but other mechanisms might be involved. Although most cell numbers returned to levels comparable to age-matched healthy donors during early and later recovery stages, several phenotypic changes in nasal immune populations persisted. For example, a subset of CD127+ granulocytes remained elevated, and continued increased activation of CD8+ TRM cells was observed during convalescence, while SARS-CoV-2-specific CD8+ T cells persisted for at least 2 months after viral clearance in the nasal mucosa. Altogether, this study provides unique insights into mucosal and systemic immune cell dynamics both during acute infection and recovery of COVID-19.

Methods

Study design and ethics

In this prospective observational cohort study, adult individuals with PCR-confirmed COVID-19 who were admitted to our academic hospital were recruited. All hospitalized individuals had hypoxia. The study was performed at the Leiden University Medical Center (LUMC), Leiden, the Netherlands, in individuals included from April 2020 to December 2020. All participants provided written informed consent. After informed consent was obtained, longitudinal sampling was performed for the duration of the hospital admission, and one convalescent sample was obtained at the outpatient follow-up appointment, which was scheduled 6 weeks after hospital discharge. Ethical approval was obtained from the Medical Ethical Committee Leiden-Den Haag-Delft (NL73740.058.20). The trial was registered in the Dutch Trial Registry (NL8589). As individuals in the ICU had substantial breathing support, we were unable to collect nasal mucosal cells from individuals in the ICU. Twelve healthy donors were included in the study. These individuals were all 60 years or older and with a male:female ratio of 2:1 to match the participant population. The healthy donors had no recent history of symptoms of airway infection (fever, cough, hypoxia, rhinorrhea, myalgia, anosmia and/or ageusia or fatigue) and were included after confirmed negative SARS-CoV-2 IgG.

Nasal cell collection and storage

Nasal cells were collected by gently scraping the nasal inferior turbinate using curettes (Rhino-Pro, Arlington Scientific) and placing them in a 15-ml Falcon tube (Corning) containing 8 ml of precooled sterile PBS containing 5 mM EDTA (Life Technologies). Per individual and time point, two curettes from one nostril were collected. Cells were dislodged by pipetting liquid up and down the tip of the curette, and cells were centrifuged at 300g for 10 min at 4 °C. The supernatant was completely removed, and cells were resuspended in 500 µl of PBS. For fixation, an equal amount of freshly prepared 8% formaldehyde (Fisher Scientific) was added, followed by a 30-min incubation at room temperature. Cells were then centrifuged at 800g for 10 min. The supernatant was completely removed, and the pellet was resuspended in 1 ml of heat-inactivated fetal bovine serum (FBS) containing 10% DMSO and moved to a cryovial. Cryovials were frozen in a Mr. Frosty freezing container (Thermo Fisher Scientific) at –80 °C and moved to liquid nitrogen within 3 d.

Whole blood cryopreservation

Per sample, 200 µl of whole blood that was collected in an EDTA tube was added to 1 ml of cold Cryostor CS10 (Stemcell Technologies) in a 1.8-ml cryovial and mixed by gently inverting. Cryovials were then placed in a cold Mr. Frosty and moved after 10 min to –80 °C. The next day, cryovials were moved to liquid nitrogen.

CyTOF staining

Samples were barcoded and measured in batches. In every batch, one aliquot of PBMCs from a reference sample was included to normalize staining between batches. Nasal cells were thawed in a water bath at 37 °C, followed by dropwise addition of 2 ml of RPMI + 50% FBS and centrifugation for 10 min at 1,600 r.p.m. at room temperature. Supernatant was discarded by pipetting. Reference PBMCs were washed with 2 ml of PBS and then fixed with 4% formaldehyde for 15 min at room temperature. Reference PBMCs were washed two times with 2 ml of BD Perm/Wash (BD). Nasal cells were washed one time with 1 ml of BD Perm/Wash, and if clumps were visible, cells were filtered through a 100-µm cell strainer (Thermo Fisher Scientific). Nasal cells and reference PBMCs were resuspended in 50 µl of Perm/Wash, and 50 µl of barcode mix targeting β2-microglobulin (B2M) was added to each individual sample in a 6-choose-3 scheme using cadmiums 106, 110, 111, 112, 114 and 116 (refs. 31,32). Samples were incubated for 30 min at room temperature and then washed with 4 ml of Cell Staining Buffer (Fluidigm). Cells were centrifuged for 5 min at 800g, and the supernatant was removed, resuspended and combined into 3 ml of Perm/Wash. Cells were centrifuged again for 5 min at 800g and were resuspended in 45 µl of Perm/Wash. FcR block (Biolegend; 5 µl) and sodium heparin (0.5 µl, 100 U ml–1) were added to prevent aspecific binding of antibodies, and cells were incubated for 20 min at room temperature33. Then, 50 µl of antibody cocktail (Supplementary Table 2) was added, followed by a 45-min incubation at room temperature. Cells were then washed twice with 2 ml of Cell Staining Buffer and centrifuged for 5 min at 800g. DNA was then stained overnight at 4 °C using 1 ml of Fix and Perm buffer (Fluidigm) containing 1,000× diluted Intercalator-Ir (Fluidigm). Cells were then washed with Cell Staining Buffer, counted and divided into tubes of 1 × 106 cells and pelleted. Pellets were washed and resuspended in cell acquisition solution (CAS, Fluidigm) with EQ Four Element Calibration Beads (Fluidigm) and acquired on a Helios mass cytometer (Fluidigm) with CyTOF Software (v7.0.8493) at the Flow Cytometry Core Facility (FCF) of LUMC in Leiden, the Netherlands (https://www.lumc.nl/research/facilities/fcf). Whole-blood samples were thawed in a water bath at 37 °C, followed by dropwise addition of 4 ml of thawing medium (RPMI 1640, penicilin/streptomycin, pyruvate, l-glutamine with 20% heat-inactivated FBS, 2 mM Mg2+ and 25 U ml–1 benzonase). Cells were centrifuged at 400g for 10 min and resuspended in 2 ml of 1× RBC lysis buffer (Biolegend). After 10–15 min, the cells were centrifuged at 400g for 10 min, resuspended in medium and filtered through a 100-µm cell strainer if clumps were visible. Staining with barcodes and surface antibodies and acquisition was performed as for nasal cells, except the sodium heparin blocking step was omitted and Cell Staining Buffer was used for whole blood instead of Perm/Wash.

Granulocyte flow cytometry analysis

Cryopreserved whole blood (one donor) or fixed and stored nasal cells (two donors) were thawed at 37 °C. Whole-blood cells were then fixed as described above for nasal cells. Then, 4 ml of BD Perm/Wash buffer was added to each sample, followed by a 5-min centrifugation at 800g. Cells were washed again in 4 ml of BD Perm/Wash and centrifuged again. Cells were resuspended in 100 µl of antibody cocktail containing CD45-BV785 (1:200; clone HI30, Biolegend, 304047), CD66b-APC-Vio770 (1:800; clone REA306, Miltenyi, 130-120-146), CD127-PercpCy5.5 (1:50; clone A019D5, Biolegend, 351321), CD16-BV711 (1:100; clone 3G8, Biolegend, 302043), Siglec8-PEDazzle594 (1:200; clone 7C9, Biolegend, 347109), CD10-PECy7 (1:400; clone HI10a, Biolegend, 312213) and CXCR4-APC (1:3,200; clone 12G5, Biolegend, 306509). Cells were incubated for 15 min and washed with 3 ml of FACS buffer (PBS, 0.5% bovine serum albumin and 2 mM EDTA). Cells were centrifuged for 5 min at 800g, resuspended in 200 µl of FACS buffer and acquired on an Aurora three-laser spectral cytometer (Cytek) using Spectroflo (v1.1) software. Single-stain controls on whole blood were used for unmixing in addition to unstained fixed nasal cells. Exported FCS files were analyzed using Flowjo X (BD).

Peripheral blood flow cytometry

Peripheral blood samples were handled according to a standard sample processing protocol for flow cytometry (for a detailed protocol, see www.EuroFlow.org and ref. 34). The combination of markers used for cell surface staining (Primary Immunodeficiency Orientation Tube (PIDOT)) has been designed along with a set of fully standardized methods and tools by EuroFlow consortium for the study of primary immunodeficiencies and was modified here by the addition of CD38 for a better separation of the circulating plasmablasts. Briefly, the procedure consisted of the bulk lysis of erythrocytes in up to 2 ml of fresh blood samples and surface staining of 2.5 × 106 white blood cells (100 µl final staining volume) with a reconstituted PIDOT lyophilized antibody cocktail (CYT-PIDOT, Cytognos) containing CD8-FITC (UCHT-4), IgD-FITC (IADB6), IgM-PerCPCy5.5 (clone MHM-88), CD16-PE (3G8), CD56-PE (C5.9), CD4-PerCPCy5.5 (RPA-T4), CD19-PeCy7 (IADB6), TCRγδ-PE-Cy7 (11F2), CD3-APC (IADB6), CD45-APC-C750 (HI30) and a drop-in antibody cocktail containing 2 µl of CD27-BV421 (clone M-T271, 562513, BDBiosciences), 2.5 µl of CD45RA-BV510 (clone HI100, 563031, BDBiosciences), 2 µl of CD38-BV605 (clone HIT-2, 303532, BioLegend) and 0.6 µl of pure CD38 (clone HIT-2, 303502, BioLegend) per test. At least 1 million events were acquired using a three-laser Cytek Aurora instrument (Cytek Biosciences; acquisition software, SpectroFlo, v1.1) from the FCF of LUMC in Leiden, the Netherlands (https://www.lumc.nl/research/facilities/fcf). The main circulating leukocyte subsets identified by flow cytometry were assigned manually using the data analysis software Infinicyt v2.04 (Cytognos SL) based on a standardized gating strategy35 with an adaptation for the CD38 inclusion (Supplementary Fig. 1). The absolute counts per microliter of fresh blood were determined by a double platform approach using the absolute fresh leukocyte counts determined before sample processing with a hematological analyzer (Sysmex) to the Statistics Configure tool of the Infinicyt software.

CyTOF data preprocessing and clustering

An outline of data preprocessing steps is shown in Extended Data Fig. 10. Debris and normalization beads were filtered from FCS files using the CyTOFclean package (v1.0.1). Single cells were then manually gated based on DNA stain and the CATALYST package (v1.12.2), and single-stain controls were used to compensate data using the non-negative linear least squares method. One-by-one plots were used to confirm correct compensation of data. Epithelial and immune cells were manually gated based on CD45 and EpCAM expression, with exclusion of cPARP+ apoptotic cells and immune doublets (CD14+CD3+, CD66b+CD3+, CD14+CD66b+). Subsequently, the CATALYST package (v1.12.2) was used to debarcode immune and epithelial cells individually per batch. FCS files were then normalized using the reference PBMCs and the CyTOFBatchAdjust function with 99th percentile scaling for each marker individually36. The marker CD69 was not present in the reference PBMCs at sufficient levels to scale and was thus not normalized. Signal intensity and clustering of reference samples before and after normalization was used to verify appropriate normalization. Clustering of cells into populations was done using hierarchical SNE or t-SNE with Cytosplore software (v2.3.0; https://www.cytosplore.org/) using all markers except EpCAM and cPARP. All t-SNE analyses were performed with complexity = 30. A diffusion map of monocytes was created using the destiny package (v.3.2.0) with k = 1,000 using the following markers: HLA-DR, CD11c, CD163, ACE-2, CD45RO, CD14, CD38, CD127, CD206, CD86, CD4, CD123 and CD45RA.

SARS-CoV-2-specific T cell isolation

PBMCs from convalescent individuals with prior COVID-19 were isolated from fresh whole blood using Ficoll-Isopaque and were cryopreserved until further use. PBMCs were thawed, and 80% was used for overnight stimulation assays. The remaining 20% was immediately used for FACS based on pMHC tetramers. For the stimulation assay, ~8 × 106 PBMCs were seeded in 1 ml of IMDM (Lonza) containing 10% heat-inactivated fetal calf serum (FCS; Sigma), 2.7 mM l-glutamine (Lonza), 100 U ml–1 penicillin (Lonza), 100 µg ml–1 streptavidin (Lonza) and in the presence of 1 µg ml–1 SARS-CoV-2 peptide pool or 1% DMSO (negative control). The SARS-CoV-2 peptide pool consisted of 15-mer peptides for all individuals and also consisted of HLA-matched 9-mer peptides when the HLA type of the individual was known (see Supplementary Table 8). The 15-mer peptides were derived from nucleocapsid (Miltenyi, 130-126-699), membrane (Miltenyi, 130-126-703) and most immunogenic sequences from the spike protein (Miltenyi, 130-126-701). Peptides were dissolved and used according to manufacturer’s protocol. Nine-mer SARS-CoV-2 peptides were predicted high binders for the 10 most common European ancestry HLA class 1 alleles and were kindly provided by P. Kvistborg (Netherlands Cancer Institute, Amsterdam, the Netherlands)37. See Supplementary Table 9 for a list of all peptides used. After 24 h, the PBMCs were washed and stained for CD4-FITC (1:30; BD, 555346, clone RPA-T4), CD8-PeCy7 (1:320; BD, 557746, clone RPA-T8), CD154-Pacific Blue (1:50; Biolegend, 310820, clone 24–31) and CD137-APC (1:100; BD, 550890, clone 4B4-1) in phenol red-free DMEM (Gibco), 2% FCS and penicillin/streptomycin for 30 min at 4 °C. Cells were sorted with a BD Aria III with BD FACSDiva software (v9.0) directly into lysis buffer and further processed as described below.

For pMHC tetramer FACS, PBMCs (2 × 106) were first incubated with a mix of in-house-produced38, PE-labeled pMHC complexes for 30 min at 4 °C before adding APC-labeled CD8 (1:64; BD, 555369, clone RPA-T8) and fluorescein isothiocyanate (FITC)-labeled CD4 (1:30; BD, 555346) antibodies at 4 °C for 30 min. The pMHC tetramers used are shown in Supplementary Table 8. Tetramer-positive CD8+CD4 T cells were sorted with a BD Aria III with BD FACSDiva software (v9.0) into 15-ml Falcon tubes and collected in lysis buffer.

TCR identification

TCRαβ sequences of T cell populations were identified as previously described with minor modifications39. Total RNA was isolated from 1 × 102 to 1 × 106 cells using the ReliaPrep RNA cell Miniprep system (Promega). The entire total RNA yield of each sample (10 μl) was mixed with 2 μl of anchored oligo(dT) primer (10 μM; Eurogentec) and 1.7 μl of SS2m_TSO primer (10 μM; Eurogentec), denatured at 72 °C for 3 min and immediately placed on ice afterward (see Supplementary Table 10 for a full list of primers used). To each sample, 7.8 μl of the first-strand reaction mix containing 0.9 μl of SMARTScribe RT enzyme (Takara, Clontech; 100 U μl–1), 0.4 μl of RNAsin (Promega; 40 U μl–1), 5 μl of 5× first-strand SMSRTScribe buffer (Takara, Clontech), 0.9 μl of DTT (100 mM; Invitrogen) and 1.7 μl of dNTPs (Promega; 10 mM) was added. The reverse transcription reaction was performed by incubating at 42 °C for 90 min followed by 10 cycles of 50 °C for 2 min and 42 °C for 2 min. Finally, the reverse transcriptase was inactivated by incubation at 70 °C for 15 min. Preamplification of the cDNA was performed on samples containing RNA from 500 or fewer cells; 12.5 μl of Phusion Flash (Thermo Fisher Scientific), 0.63 μl of IS PCR primer (10 μM; Eurogentec) and 1.87 μl of water (UltraPure water, Invitrogen) were added to 10 μl of the previously synthesized cDNA. The IS primer anneals both to the SA.rt and oligo(dT) IS region to amplify the complete cDNA40. The reaction was incubated at 95 °C for 2 min and cycled 18 times between 95 °C for 1 s, 69 °C for 15 s and 72 °C for 2 min, with a final extension at 72 °C for 5 min. Barcoded TCR PCR product was generated in two rounds of PCR. In the first PCR, TRA and TRB product was generated in separate PCR reactions for α- and β-chains as follows: 25 μl of Phusion Flash (Thermo Fisher Scientific), 1.25 μl of SS2m_For PCR primer (10 μM; Eurogentec), 1.25 μl of TRAC_rev or 1.25 μl of TRBC1 + 2 mix PCR primer (10 μM; Eurogentec) and 17.5 μl of water (UltraPure water, Invitrogen) added to either 5 μl of cDNA or 5 μl of preamplified cDNA. The reaction was incubated at 98 °C for 2 min and cycled 30 times between 98 °C for 1 s, 67 °C for 15 s and 72 °C for 15 s, with a final extension at 72 °C for 1 min. The PCR product was 96-well plate purified with the Wizard SV 96 PCR Clean-Up System (Promega) and eluted in 70 μl of water. In a second PCR, the first purified PCR product was used to include a two-sided six-nucleotide barcode sequence that allows for discrimination between TCRs of different T cell populations. Then, 20 μl of Phusion Flash (Thermo Fisher Scientific), 1.6 μl of BC_R7xx_For bc PCR primer (10 μM; Eurogentec), 1.6 μl of BC_R7xx_Rev bc PCR primer (10 μM; Eurogentec) and 6.8 μl of water (UltraPure water, Invitrogen) were added to 10 μl of purified PCR product. The reaction was incubated at 98 °C for 2 min and cycled 10 times between 98 °C for 1 s, 65 °C for 15 s and 72 °C for 30 s, with a final extension at 72 °C for 2 min. The barcoded PCR product was also purified with the Wizard SV 96 PCR Clean-Up System (Promega) and eluted in 70 μl of water. PCR products of different T cell populations were pooled, after which TCR sequences were identified by NovaSeq (GenomeScan). NovaSeq data were analyzed using MiXCR software (v3.0.13) to determine the Vα and Vβ family and CDR3 regions using annotation to the IMGT library (http://www.imgt.org; v6). CDR3 regions were analyzed in RStudio, and CDR3 sequences with ≤50 reads that were non-functional or occurred on all samples were excluded from the analysis.

Statistics

Statistical differences in cellular abundance between groups were compared with a linear mixed model in which individuals were included as random effect and groups (acute, ERS or all hospitalized individuals, convalescent and healthy donors) as fixed effect using the lme4 (v.1.1–23) and lmerTest (v3.1–2) packages. Post hoc comparison of all groups included in a linear mixed model was conducted with the emmeans package (v1.4.8) using the Tukey correction for multiple comparisons. For the comparison of multiple subsets or lineages, Benjamini–Hochberg multiple testing correction was subsequently used. Correlation values were calculated with Pearson tests. MDS was performed in R using the dist function in the stats package (4.0.1). Features were normalized before MDS with labels shown for selected clusters. If viral load data were not available for the time of nasal cell sampling, the closest time point with available data was used. Viral load is calculated as 40 – Ct of SARS-CoV-2 PCR, meaning that higher values indicate more virus. MOFA was conducted using the MOFA package (v1.4.0) using standard parameters, with the exception of DropFactorThreshold of 0.02 and Tolerance of 0.01 (ref. 22). Features were scaled before use in MOFA, and the model was repeated three times to qualitatively assess stability of LFs and defined features. All analyses were performed with R version 4.0.1, except for the cytofclean package performed in R3.6.3 using RStudio (v1.2.5033).

Reporting Summary

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