The ethical approval for the study was obtained from CMO Arnhem-Nijmegen (NL32 357.091.10). All participants gave written consent for their participation in the study. In total, 14 convalescent COVID-19 patients and 15 sex-matched healthy donors (controls) from Radboudumc, the Netherlands, participated in this study. Convalescent patients were included in the study based on self-reported symptoms and subsequent confirmation of the presence of IgA, IgG, and IgM antibodies against the SARS-CoV-2 Spike protein in the circulation. All controls were at healthy status during sampling and did not experience COVID-19 based on antibody level ( Figure S1 ). At the time of inclusion, patients had no reported COVID-19 symptoms or any other infections. Inclusions took place between April and June of 2020, during the first wave of the pandemic in the Netherlands. The demographic and clinical characteristics of the participants are provided in Table 1 .

Whole blood was diluted with phosphate-buffered saline (PBS), and PBMC isolation was performed by density gradient centrifugation using Ficoll-Paque (GE Healthcare, IL, USA). The PBMC fraction was collected and washed three times with PBS. The cells were resuspended in RPMI 1640 Medium (Dutch modification) (Thermo Fisher Scientific, MA, USA) supplemented with 1 mM sodium pyruvate (Thermo Fisher Scientific), 2 mM GlutaMAX supplement (Thermo Fisher Scientific), and 5 µg/mL gentamicin (Centraform, Netherlands). For DNA methylation, single cell RNA sequencing (scRNA-seq), and single cell ATAC sequencing (scATAC-seq) analysis, and future stimulation assays, cells were cryopreserved in RPMI containing 40% bovine calf serum (Cytiva, MA, USA) and 15% dimethyl sulfoxide (VWR, PA, USA).

After erythrocyte lysis in isotonic NH4CL buffer and washing twice with PBS, total leukocytes were obtained. White blood cell counts were determined by a cell counter (Coulter Ac-T Diff^® cell counter; Beckman Coulter, CA, USA) and used to calculate the absolute numbers of CD45+ leukocytes identified by flow cytometry as previously described in detail (18). Briefly, half a million of total leukocytes per staining panel were used to analyze surface markers with the Navios™ flow cytometer (Beckman Coulter). Cells were transferred to a V-bottom 96-well plate and washed twice with PBS containing 0.2% bovine serum albumin (BSA; Sigma-Aldrich, St. Louis, USA), stained for 20 minutes at room temperature in the dark with the staining panels of interest and washed twice with PBS containing 0.2% BSA. The conjugated antibodies specific for human cells are given in Table S1 . The gating strategy used in flow cytometry analyses were provided in Figure S2 . Data analysis was performed using the Kaluza 2.1^® software (Beckman Coulter).

Frozen cells were rapidly thawed at 37 °C and transferred into 50 mL centrifuge tubes. Subsequently, cells were counted using an automated cell counter (Thermo Fisher Scientific). An equal number of cells (3,300 per individual) from 4 different individuals were pooled together and loaded into a 10X Chromium Controller to generate Gel Beads-in-emulsion (GEMs). scRNA-seq libraries were generated according to the manufacturer’s instructions (Chromium Next GEM Single Cell 3’ Reagent Kits v3.1 (Dual Index) User Guide, Rev A, CG000315 Rev A). In brief, all the loading cells were separated into nanoliter-scale droplets. Within each droplet, cDNA was generated via reverse transcription, adding a cell barcoding sequence and unique molecular identifier (UMI) to each cDNA molecule. Library quality per pool was examined using the Agilent Bioanalyzer High Sensitivity DNA kit. Sequencing was carried out on NovaSeq 6000, having a 50,000 reads depth per cell. DNA was isolated from PBMCs and then subjected to genotyping by Illumina Global Screening Array Beadchip to demultiplex the pooled samples.

The Cellranger (v.4.0.0) utility of 10X genomics was used to process the scRNA-seq data. The standard protocol of the Cellranger pipeline consists of mapping sequenced reads and measuring gene expression (19). The human reference genome (GRCh38) was used for mapping sequence reads. Later, demultiplexing of individuals mixed in the same pool was done by using the genotype-free clustering method souporcell (20). Further downstream analyses were performed using the Seurat package V4 (21) in R following its default pipelines, i.e., quality filtering, data normalization, and scaling. During quality filtering, mitochondrial and ribosomal genes along with cells having mitochondrial reads > 25% and expressed genes > 5000 or < 250 were removed. Data normalization was done using a global-scaling normalization method, namely “LogNormalize”. Subsequently, data was scaled prior to dimensionality reduction through principal component analysis (PCA) using the top 2,000 variable features detected using the vst method. Uniform Manifold Approximation and Projection (UMAP) technique was utilized for cell clustering. The “FindAllMarkers” function was used for identifying marker genes for each cluster. Using publicly available information, each cluster was annotated based on top marker genes (having the lowest adjusted p-value and positive logFC). Differentially expressed genes (DEGs) between convalescent COVID-19 patients and controls were detected through the “FindMarkers” function using the MAST test (22) adjusted by age and sex. Genes that were expressed in at least 10% cells and with adjusted P-value (after Bonferroni post-hoc correction) < 0.05 were considered significant. Gene Ontology (GO) enrichment analysis of significant DEGs was done using the “clusterProfiler” package in R (23). Later, for cross-validation, a matrix containing the total gene count per sample was generated from Seurat by summing up all the read counts per gene, and the pseudo-bulk RNA analysis was performed through DESeq2 (24).

Cryopreserved PBMCs were thawed at 37°C and transferred into 50 mL centrifuge tubes. PBMCs were washed in PBS through centrifugation at 400G for 5 minutes at 4 °C and lysed for 3 minutes on ice. After discarding the supernatant, lysed cells were diluted within 1× diluted nuclei buffer (10X Genomics) prior to counting with trypan blue using a Countess II FL Automated Cell Counter to validate lysis. An equal number of nuclei (approximately 3000 nuclei per sample) from four individuals were pooled together and then loaded into the Chromium Next GEM Chip H based on the user guides from 10X genomics (Chromium Next GEM Single Cell ATAC Reagent Kits v1.1 User Guide, CG000209 Rev D). After breaking the emulsion, the barcoded tagmented DNA was purified and amplified for sample indexing and generation of scATAC-seq libraries. The final libraries were quantified using the Agilent Bioanalyzer High Sensitivity DNA kit. Sequencing was performed on NovaSeq 6000 with a depth of 25,000 reads per nuclei.

Using the cellranger-atac mkfastq (10X Genomics, v.1.2.0) utility, raw sequencing data were converted to FASTQ format. Reads were aligned based on the human reference genome (GRCh38). Demultiplexing and doublets removal was done using souporcell (20). Each scATAC-seq library fragment file was utilized further for downstream analysis. Employing the ArchR package (v.1.0.1) within R, quality control (cells having a transcription start site (TSS) enrichment score < 4 and unique nuclear fragments < 1000 were filtered), dimensionality reduction (iterativeLSI function, default parameters), harmony batch correction (25) (addHarmony function, default parameters), clustering (addClustersfunction, method = “Seurat”, resolution = 1, UMAP) and cluster visualization (addUMAP function) were done (26). The default parameter of ArchR was used to calculate gene expression by computing a “gene score’’ based on chromatin accessibility within a gene body and distally & proximally from the transcription start site (TSS). A gene score matrix was created using the gene score profiles for all cells. Later, the wilcoxon rank-sum test was used to identify cell type specific marker genes. Gene having adjusted P-value (Benjamini-Hochberg) < 0.05 and logFC >= 0.05 were considered as marker genes. Each cluster was annotated based on cell type specific marker genes and cluster specific marker genes retrieved from scRNA-seq data. Further, using the method (“addGeneIntegrationMatrix” function) developed by Stuart (27), the gene score matrix was integrated with scRNA expression data. MACS2 (28) strategy (getMarkerFeatures function, useMatrix = “PeakMatrix”) was adapted on an integrated dataset to call peaks in each cluster. Peaks having logFC greater than 0.05 and P-value < 0.05 were considered as marker peaks. For cross-validation, a peak matrix per sample was generated, and the pseudo-bulk ATACseq analysis was performed through ArchR.

Six healthy controls and six convalescent COVID-19 patients were used to assess interferon response. Cryopreserved PBMCs were stimulated with 10 ng/ml recombinant human IFN-α2 (R&D Systems), 1000 U/ml IFN-β (R&D Systems), and 50 ng/ml IFN-γ1b (Miltenyi Biotec, Germany) for 24 hours. Stimulations were performed in U-bottom 96-well tissue culture plates with 5x10⁵ PBMCs per well. RNA isolation was carried out using the RNeasy Mini kit (Qiagen, Germany), and cDNAs were generated with the iScript cDNA synthesis kit (Bio-Rad Laboratories, CA, USA) according to the instructions of the manufacturers. RT-qPCR for the interferon-response genes IFI44L, IFI6, IRF3, IRF7, IRF9, ISG15, and OAS2 was performed with StepOnePlus PCR System (Applied Biosystems, MA, USA). HPRT1 was used as the housekeeping control gene. Primers and reaction conditions are provided in Tables S2, S3 .

DNA was isolated from whole blood using QIAamp DNA Micro Kit (Cat: 56304, Qiagen, Hilden, Germany), and the concentration was determined using a NanoDrop spectrophotometer at 260 nm. The high-quality DNA from 20 samples (11 convalescent COVID-19 and 9 controls) were obtained successfully for the genome-wide DNA methylation profiling through the Illumina Infinium^© MethylationEPIC array (~850,000 CpG sites). The DNA methylation values were gained from the raw IDAT files using the minfi package in R (v.4.0.3) (29). Initially pre-processing was performed to filter bad quality probes with a detection P-value > 0.01, cross-reactive probes, polymorphic probes (30), and probes in the sex chromosome. Subsequently, after stratified quantile normalization (31), the comparison was made between convalescent COVID-19 individuals and controls to detect the relative proportion of cell types and differentially methylated CpG sites. Based on methylation value, cell proportion was estimated using modified Housman’s method implemented in the estimateCellCounts2 function of the FlowSorted.Blood.EPIC R package (32). Later, differential analysis of cell-type proportions was done using the regression model of the DirichletReg package of R. Differentially methylated CpG sites were detected by a linear regression model using the limma R package (33), with age and sex as covariates. The Spearman correlation was calculated between the cell proportion and the top 150 CpG differentially methylated sites. DMRs analysis was done by using combp R package (34).

5x10⁵ PBMCs per well were seeded in U-bottom 96-well tissue culture plates (Greiner Bio-One, Austria) and stimulated with 10 ng/ml lipopolysaccharides (LPS) derived from Escherichia coli O55:B5 (Sigma-Aldrich, MO, USA), 10⁶/ml heat-killed Staphylococcus aureus, 10⁶/ml heat-killed Candida albicans, R848 (InvivoGen), or heat-inactivated(30 minutes, 60°C) SARS-CoV-2 Wuhan-Hu-1 strain (NRW-42 isolate, GISAID accession number: EPI_ISL_425126) (35) with 1 µg/ml soluble purified mouse anti-human CD28 (BD Biosciences, New Jersey, USA), for 24 hours and 7 days in the presence of 10% pooled human serum. A control condition without any stimulant was included. Cells were incubated at 37°C with 5% CO2. Stimulations with R848 and SARS-CoV-2 were performed later with cryopreserved cells, while the rest of the stimulations were done with freshly isolated PBMCs. Concentrations of TNF-α, IL-6, IL-1β, IL-1Ra, and IFN-γ in culture supernatants were determined using DuoSet ELISA kits (R&D Systems, MN, USA). IFNα concentrations were determined using VeriKine Human Interferon Alpha ELISA Kit (PBL Assay Science, NJ, USA) according to the manufacturer’s instructions. Supernatants from 7 days stimulations were used for IFN-γ, and the remaining cytokines were measured in 24-hour supernatants.

To characterize the immunological features of convalescent COVID-19 patients, we collected blood from 14 convalescent COVID-19 patients and 15 healthy donors. At the time of sampling, none of the control individuals experienced any symptoms or signs of infection, including COVID-19. All convalescent COVID-19 patients had recovered from SARS-CoV-2 infection. Out of the 14 patients, only one convalescent COVID-19 patient had experienced a severe COVID-19 and was admitted to an intensive care unit (ICU), while the remaining 13 experienced mild infections. The recovery time for most of the convalescent COVID-19 patients (9 out of 14) was ca. one month, only one patient with one-week recovery time and two patients with around two-month recovery time. Sex was distributed almost equally in both convalescent and healthy individuals. The mean age of convalescent COVID-19 patients and control individuals were 60 ± 11 years and 50 ± 14 years, respectively ( Table 1 ).

To explore the effect of the SARS-CoV-2 infection on the immune system after recovery, we have performed a comprehensive set of assays and analyses, which were summarized in Figure 1A .

Since acute SARS-CoV-2 infection leads to striking changes in immune cell types in the blood, flow cytometry analysis was performed to assess whether the abundance of various immune cell subsets in convalescent COVID-19 individuals is still different from those of healthy controls. The absolute numbers and percentages of most immune cell subsets, such as neutrophils, natural killer (NK) cells, B cells, and T cells, were similar between convalescent patients and controls ( Figures S3A–C ). Although COVID-19 leads to a dramatic decline in CD4+ T and CD8+ T cells in the acute phase, these cell populations in convalescent patients were not significantly different from those of uninfected people, indicating the full recovery from lymphopenia ( Figure 1B ). We observed that regulatory T cells (Treg) and CD4+ naive/terminally differentiated effector memory cell (naïve/TEMRA) populations were somewhat more abundant in the convalescent COVID-19 individuals compared to controls, but these differences were not statistically significant ( Figure 1B ).

Next, we performed an integrative analysis of scRNA-seq and scATAC-seq in convalescent and healthy individuals. For scRNA-seq analysis, we examined 66,753 single cells from 21 individuals (12 convalescent and 9 controls) after QC ( Figure S4A ). Uniform manifold approximation and projection (UMAP) was used for cell clustering. In total, we detected 11 immune cell types, namely, CD4+ T cells (IL7R+), CD8+ T cells (CD8A+CD8B+GZMK+), B cells (CD79A+), NK cells (NKG7+GZMB+), classical monocytes (CD14+LYZ+), non-classical monocytes (FCGR3A+), platelets (PPBP+), and monocyte-derived dendritic cells (mDCs, HLA-DP^high, HLA-DR^high) ( Figures 2A, B ). Subsequently, we analysed 38,186 nuclei from 16 individuals (9 convalescent and 7 controls) in scATAC-seq analysis ( Figures 2C, D and Figure S4B–D ). It is pertinent to note that a strong correlation of the marker genes between scRNA-seq and scATAC-seq analysis was detected ( Figure S4E ). Subsequently, we compared the cell compositions of both scRNA-seq and scATAC-seq. No statistically significant difference was detected between the convalescent COVID-19 patients and controls ( Figures 2E, F ). This supports the findings of the flow cytometry measurements showing full recovery of the immune cell subsets.

To identify if any differences exist at the transcriptomic level between immune cells of convalescent COVID-19 patients and controls, we performed a differential expression analysis. After correcting for age and sex, 907 differentially expressed genes (DEGs) were identified with false discovery rate (FDR) < 0.05 and absolute log-fold change (logFC) > 0.05 ( Figure S5A ). DEGs were found predominantly in monocytes (both classical (DEGs=579) and non-classical (DEGs=62)) and CD4+ T (DEGs=172) cells. Moreover, we assessed the number of DEGs at different logFC thresholds: we identified 293 DEGs with FDR < 0.05 and absolute logFC > 0.1, and 30 DEGs with FDR < 0.05 and absolute logFC > 0.25. Although the number of DEGs decreased with higher logFC thresholds, most of them predominantly existed in monocytes and CD4+ T cells ( Figure S5A ). This indicates that the observed differences in the transcriptional profile of convalescent COVID-19 patients were mainly in monocytes and CD4+ T cells. Additionally, we compared the peak accessibility in each cell type based on scATAC-seq data. After multiple testing correction, none of the peaks were significant, which indicated that the open chromatin accessibility between convalescence patients and healthy controls did not show significant difference. The differential expressed genes were mainly found in classical monocytes. In classical monocytes, for peaks in the proximity (500kb) of up-regulated DEGs in convalescence, 7205 peaks were up-regulated (mean log2FC = 1.077) in convalescence while 7371 peaks were down-regulated (mean log2FC = -1.073). For peaks in the proximity (500kb) of down-regulated DEGs in convalescence, 10252 peaks were down-regulated (mean log2FC = -1.078) in convalescence while 9777 peaks were up-regulated (mean log2FC = 1.086). None of these peaks were significantly different between convalescence and controls (P-value _{FDR adjusted} > 0.5).

Differential expression analysis in the scRNA-seq revealed that MHC class II genes, such as HLA-DRB5, HLA-DPB1, HLA-DPA1, HLA-DRB1, and HLA-DRA, were downregulated in convalescent COVID-19 individuals ( Figure 3A ). Moreover, antigen processing and presentation via the MHC II pathways were downregulated in CD4+ T cells ( Figure 3B ). On the other hand, no significant difference in the surface expression of HLA-DR in CD4+ T cells and monocytes between healthy and convalescent individuals was observed ( Figure S5B ).

Since we identified DEGs mainly in classical monocytes ( Figure S5A ), we performed a sub-clustering analysis of these cells. Initially, using UMAP, we captured 6 sub-clusters in classical monocytes. Further inspection revealed that one healthy control, HC04, was markedly different from others ( Figures S5C, D ). Classical monocytes-specific pseudo-bulk RNA analysis also showed that HC04 was separated from the main group in a PCA plot ( Figure S5E ). Therefore, we temporarily removed HC04 while performing differential expression analysis within the classical monocytes. After discarding HC04, we re-performed sub-clustering and detected five sub-clusters ( Figure 3C ). The top markers for each sub-cluster are shown in Figure 3D . DEGs between convalescent individuals and controls mainly occurred in cluster 0 (S100A^high) and cluster 1 (IL7R^high) with FDR < 0.05 and logFC > 0.05 ( Figure S5F ). Subsequently, GO enrichment analysis was performed using the DEGs present in cluster 0 (812 DEGs) and cluster 1 (33 DEGs), respectively. The upregulated DEGs in convalescent individuals were enriched in antigen processing and presentation via MHC I. while the down-regulated DEGs were enriched in response to the interferon pathway in the S100A^high classical monocytes. No enriched pathway was found in DEGs from cluster 1 ( Figure 3E ). The genes engaged in the down-regulated interferon response pathway included interferon regulatory factor 3 (IRF3), IRF7, IRF9, 2’-5’-Oligoadenylate Synthetase 2 (OAS2), interferon-alpha inducible 6 (IFI6), interferon-alpha inducible ligand 44 (IFI44L), interferon-stimulated gene 20 (ISG20), and ISG15 ( Figure 3F ). These genes were thus selected for the in vitro validation of the sequencing results.

IFN production constitutes the major first line of defense against viruses. Type I and type II IFNs induce hundreds of antiviral effectors, or ISGs, to achieve a cell-intrinsic state of viral resistance (36). Since IFN response pathway genes were downregulated in monocytes of convalescent patients, we hypothesized that expressions of those genes would be different between healthy and convalescent COVID-19 individuals after stimulation with IFN-α, IFN-β, and IFN-γ. Therefore, PBMCs of 6 healthy and 6 convalescent individuals were stimulated with type I and II interferons, and genes involved in IFN response were analyzed via RT-qPCR. The baseline gene expressions were not different between controls and convalescent individuals ( Figure S5G ). On the other hand, we found that expression of three ISGs, namely, IRF3, ISG15, and IFI6, were significantly higher in convalescent COVID-19 individuals after stimulation with IFN-β, IFN-α/β, and IFN-β, respectively ( Figure 3G ), indicating that the degree of induction by type I IFNs is higher for some of the IFN regulated genes in convalescent individuals compared to healthy volunteers. We performed pseudo-bulk RNA analysis between convalescent COVID-19 individuals and controls and found no significant difference, in line with the RT-qPCR results from non-stimulated PBMCs ( Table 2 ). Consistently, pseudo-bulk ATAC analysis also revealed no significant difference in the chromatin accessibility of ISGs associated with the interferon response pathway in PBMC between conditions ( Table 2 ). These results suggest that the down-regulation of the interferon pathway is specific to monocytes, while at the PBMC level, this pathway is not changed.

Earlier studies have reported that the host cell epigenetic landscape of DNA methylation changes during SARS-COV2 infection (14). Therefore, we performed DNA methylome profiling from high-quality DNA of 20 whole blood samples (11 convalescent COVID-19 and 9 controls). Initially, pre-processing of the dataset scanned 794,745 good-quality probes. Cellular deconvolution analysis using modified Housman’s method (32, 37, 38) detected six cell-types, namely NK, CD8+ T, CD4+ T, B cells, neutrophils, and monocytes, which were present in both convalescent COVID-19 individuals and controls. The relative proportion of these six cell types did not significantly differ ( Figure 4A ). This result, consistent with our flow cytometry, scRNA-seq and scATAC-seq analyses, indicated that the cell proportions returned to normal levels after recovery from SARS-CoV-2 infection.

Furthermore, epigenome-wide association analysis (EWAS) was performed to assess the difference in genome-wide DNA methylation profiles between convalescent COVID-19 individuals and controls. No significant differentially methylated CpG sites (DMSs) (FDR < 0.05) were identified ( Figures S6A-C ). Subsequently, to detect any minor changes at the molecular level, we scanned the top 50, 75, 100, and 150 DMSs. A prevalence of hypermethylated and demethylated CpG sites was consistently observed at these top DMSs in convalescent COVID-19 patients compared with those in healthy controls ( Figure 4B ). Thus, genome-wide DNA methylation profiles of convalescent COVID-19 individuals show only minor differences compared to controls.

Next, we aimed to detect individual cell types where the methylation signal was the strongest. We performed a correlation analysis between the top 150 DMSs (94 downregulated and 56 upregulated) having the lowest P-value and the estimated cell proportions from the above-mentioned deconvolution analysis. We observed that the upregulated CpG sites were positively correlated with CD4+ T cells, monocytes, and neutrophils ( Figure 4C ). Conversely, the downregulated CpG sites showed a negative correlation with CD4+ T cells, monocytes, and neutrophils ( Figure 4D ), suggesting that majority of the top differentially DNA methylation signatures in convalescent COVID-19 individuals mainly originated from monocytes and CD4+ T cells, and somewhat less neutrophils. These results could be replicated using the top 100 or top 75 DMSs ( Figures S6D, E ). However, we could not identify enriched epigenetic changes around the DEGs (+/- 250kb window) in monocytes observed from the scRNA-seq, compared to randomly selected CpG sites ( Figure S6F ).

Differential methylation regions (DMRs) have identified 30 significant DMRs between convalescent COVID-19 individuals and controls ( Table S4 ). The genes annotated to these DMRs were enriched in the following biological pathways: glycosphingolipid biosynthesis, telomere maintenance, recognition of DNA damage, sialic acid metabolism, chromosome maintenance and TGF-Ncore.

Lastly, we measured the cytokine production capacity of PBMCs upon incubation with different viral, bacterial, and fungal stimuli to assess whether the differences in transcriptome influence immune responses upon recovery from COVID-19. PBMCs were either left unstimulated or stimulated with LPS, S. aureus, C. albicans, viral RNA mimic R848, and heat-inactivated SARS-CoV-2 (original Wuhan-Hu1 strain) for 24 hours to measure TNF-α, IL-6, IL-1β, IFNα and IL-1Ra and seven days to measure IFN-γ. Anti-CD28 agonistic antibodies were combined with inactivated SARS-CoV-2 to better activate T cell responses, particularly in convalescent COVID-19 patients. TNF-α and IL-6 production upon different stimuli were similar between healthy individuals and convalescent COVID-19 patients ( Figures 5A, B ). Another proinflammatory cytokine, IL-1β, was also produced similarly in healthy and convalescent individuals upon R848 and S. aureus stimulation ( Figure 5C ). Furthermore, R848 and SARS-CoV-2 induced IFNα production from PBMCs did not differ between the groups ( Figure 5D ). Intriguingly, baseline IL-1Ra secretion of convalescent COVID-19 patients was significantly higher than healthy controls, although IL-1Ra production following R848 or heat-inactivated SARS-CoV-2 stimulation was similar ( Figure 5E ). Furthermore, PBMCs from convalescent individuals produced significantly more IFN-γ without any stimulation ( Figure 5F ). As expected, incubation of convalescent COVID-19 patient PBMCs with heat-inactivated SARS-CoV-2 resulted in higher IFN-γ production, indicating the presence of T cell memory after recovery. Lastly, IFN-γ secretion did not differ between healthy and convalescent upon incubation with viral RNA mimic R848 ( Figure 5F ).