MATERIALS AND METHODS

Human subjects
The study was approved by the Ethics Committee of Zhongshan Ophthalmic Center, China and the Ethics Committee of Wuhan Hankou Hospital, China. A written informed consent was routinely obtained from all individuals participating in the study and all relevant ethical regulations regarding human research participants were followed. Healthy non-frail individuals were recruited in the Zhongshan Ophthalmic Center, and divided by age into two groups in cohort-1: young adults (YA) and aged adults (AA). The YA group ranged from ages 20 to 45 years old and the AA group ranged from ages 60 to 80 years old. COVID-19 patients diagnosed by real-time fluorescent quantitative reverse transcription polymerase chain reaction (RT-qPCR) and CT images were enrolled in the Wuhan Hankou Hospital, China. Based on their clinical history, patients were divided into incipient and recovered groups in cohort-2 and cohort-3 respectively, and the incipient hospitalized patients were further divided by age into young COVID-19 patient onset (YCO) and aged COVID-19 patient onset (ACO). Enrolled patients that tested negative with nucleic acid transfer in 7–14 days were further divided into young COVID-19 patient recovered (YCR) and aged COVID-19 patient recovered (ACR). Individuals with comorbid conditions including cancer, immunocompromising disorders, hypertension, diabetes and steroid usage were excluded. No significant gender differences were detected between YA group and AA group in cohort-1 (Table S1C–E), between YH, AH, YCO and ACO group in cohort-2 (Table S1F), between YH, AH, YCR and ACR group in cohort-3 (Table S1G).

Antibodies and reagents
Antibodies against the following markers in flow cytometric analysis were purchased from Biolegend, BD biosciences and Abcam: CD3 (clone SK7) BV785 (Cat. 344842), CD19 (clone HB19) APC (Cat. 302212), CD88 (clone S5/1) PE/Cy7 (Cat. 344307), CD89 (clone A59) PE/Cy7 (Cat. 354107), HLA-DR (clone L243) FITC (Cat. 307604), CD11c (clone 3.9) BV421 (Cat. 301627), FcεRIa (clone AER-37) PercP/Cy5.5 (Cat. 334622), CD1c (clone L161) BV650 (Cat. 331541), CD371 (CLEC12A) (clone 50C1) PE (Cat. 353603) were purchased from Biolegend, CD147 (clone HIM6) PE (Cat. 562552) was purchased from BD biosciences. Fetal bovine serum (FBS) (Cat. 10270-106), penicillin/streptomycin (Cat. 15140-122), and Trypsin-EDTA (0.25%) (Cat. 25200-072) were purchased from GIBCO. RT-qPCR kit (Cat. 25200-072) was purchased from TaKaRa.

Detection of SARS-Cov-2 with RT-qPCR
Samples used for RT-qPCR were blood, upper respiratory tract sputum and throat swab obtained from patients at specified time-points during hospitalization. The patient samples were collected, processed and analyzed following the guideline stipulated by the WHO. To extract viral RNA, the specimens were treated with the QIAamp RNA Viral Kit (Qiagen, Heiden, Germany) following the manufacturer’s guidelines. The presence of SARS-CoV-2 infection was confirmed with a China CDC recommended RT-qPCR kit (TaKaRa, Dalian, China). qPCR was performed as previously described (Zhang et al., 2019; Bi et al., 2020; Li et al., 2020).

Isolation of PBMCs for mass cytometry, scRNA-seq and scATAC-seq
For pipeline analysis, venous blood samples were derived from each healthy donor or patient using Ficoll-Hypaque density solution, heparinized and then processed by standard density gradient centrifugation methods to isolate PBMCs. The viability and quantity of PBMCs in single-cell suspensions were determined using Trypan Blue. For each sample, the cell viability exceeded 90%. For each sample with more than 1 × 107 viable cells, a fraction of PBMCs was extracted for scRNA-seq analysis, a fraction of PBMCs was allocated for scATAC-seq and mass cytometry.

Flow cytometric analysis
PBMCs suspended in phosphate buffered saline (PBS) were cultured with Live/Dead yellow dye (Invitrogen) at 4 °C for 30 min and then washed once with 1 mL of PBS containing 1% FBS (GIBCO, Grand Island, NY, USA). Subsequently, cells were treated with antibodies for 30 min at 4 °C. These antibodies included: CD3-BV785 (clone SK7, Biolegend), CD19-APC (clone HB19, Biolegend), CD88-PE/Cy7 (clone S5/1, Biolegend), CD89-PE/Cy7 (clone A59, Biolegend), HLA-DR-FITC (clone L243, Biolegend), CD11c-BV421 (clone 3.9, Biolegend), FcεRIa- PercP/Cy5.5 (clone AER-37, Biolegend), CD1c-BV650 (clone L161, Biolegend), CD147-PE (clone HIM6, BD biosciences), CD371 (CLEC12A)-PE (clone 50C1, Biolegend). Analysis of PBMCs with flow cytometry was conducted with BD Fortessa (BD Biosciences) and the results were evaluated with FlowJo (version 10.0.7, Tree Star, Ashland, OR, USA).

Mass cytometry live cell barcoding and surface staining
We made use of a live cell barcoding approach to minimize inter-sample staining variability, sample handling time and antibody consumption. After incubating with anti-human CD45 loaded with different isotopes (89Y, 162Dy, 165Ho, 169Tm, 175Lu), all the samples were then pooled for surface staining. The Maxpar Direct Immune Profiling Assay (Fluidigm) was used for cell surface staining and the monoclonal anti-human antibodies in the assay kit are listed as Table S2.
Barcoded and combined samples were washed and stained with viability dyes cisplatin-195pt (0.5 μmolL) (Fluidigm, 201064) and vortexed to mix thoroughly for 2 min at room temperature for cell viability, terminated with Maxpar Cell Staining buffer at room temperature (400 rcf.), washed, fixed with 1.6% paraformaldehyde (PFA; Electron Microscopy Sciences) in PBS for 10 min at room temperature on a rotary shaker (500 rpm). The fixed cells were resuspended in pre-cooling Maxpar Cell Staining to slow fix reaction. Fixed samples were washed twice with PBS/bovine serum albumin and once with double-distilled water before resuspended in 400 μL of surface-antibody mixture. Surface staining was performed for 30 min at 37 °C on a rotating shaker (500 rpm). The samples then stored in freshly diluted 2% formaldehyde (Electron Microscopy Sciences) in PBS containing 0.125 nmol/L iridium 191/193 intercalator (Fluidigm, 201192) at 4 °C overnight.

scRNA-seq data alignment, processing and sample aggregation
The Chromium Single Cell 5′ Library (the 10x Genomics chromium platform Illumina NovaSeq6000), Gel Bead and Multiplex Kit, and Chip Kit (10x Genomics) were used to convert single-cell suspension samples to barcoded scRNA-seq libraries. Single-cell RNA libraries were prepared using the Chromium Single Cell 5′ v2 Reagent (10x Genomics, 120237) kit as per the manufacturer’s protocols. The quality of the libraries was checked using the FastQC software. Initial processing of the sequenced data was performed using CellRanger software (https://support.10xgenomics.com, version 3.1.0).
The command Cell Ranger count in CellRanger Software Suite (10x Genomics) was used to demultiplex and barcode the sequences derived from the 10x Genomics single-cell RNA-seq platform. The data was filtered, normalized, dimensionality was reduced, clustered, and differential gene expression analysis were performed after calculation of the single-cell expression matrix by CellRanger using Python (version 3.7.7) Scanpy (https://scanpy.readthedocs.io/en/stable/index.html, version 1.4.6). Data collection and the subsequent analyses were performed in an unsupervised manner, but not blinded to the conditions of the experiments. For quality control, the filtered cell population was mainly those cells that express HBB, HBA1, and several light and heavy chain transcripts, which identified as the RBC-contaminated cell population. Likewise, several clusters expressing genes has no significance (P ≥ 0.1, calculate by 10x genomics Loupe Cell Browser with it default algorithm. P values are adjusted using the Benjamini-Hochberg correction for multiple tests) were removed. A total of 16 libraries were sequenced, and 166,609 cells (YA 77,652 cells, AA 88,957 cells) were analyzed after quality control in cohort-1. For cohort-3, 22 libraries and 205,434 cells (YH 79,039 cells, AH 88,750 cells, YCR 19,533 cells, ACR 18,112 cells) were remained for the subsequent analysis. The genes used in principal component analysis (PCA) analysis have eliminated mitochondria (MT), and ribosomes (RPL and RPS) genes with 50 principal components, and then aligned together, followed by t-distributed stochastic neighbor embedding (t-SNE) are both used after the results of the aligned. And using the run_harmony function (in pyharmony package, version 1.0.7) and combat function (in Scanpy) methods to deal with batch effect issues if batch effect existing in dataset. Genes not detected in any cell were removed from subsequent analysis.

Dimensionality reduction and clustering analysis of scRNA-seq datasets
To analyze the scRNA-seq data, we log normalized data (1 + counts per 10,000) with the ‘‘sc.pp.normalize_total’’ function before clustering, reduction and performing 2-dimensional t-SNE algorithm clustering with the first 50 principal components. This was done following PCA on top 5,000 most variable genes by using “sc.pp.highly_variable_genes” function in Scanpy with the default parameters. Dimensionality method and identification of significant clusters and was performed using Leiden clustering algorithm which uses a shared nearest neighbour modularity optimization-based clustering algorithm. Marker genes for each significant cluster were found using the function sc.tl.rank_genes_groups with default parameters.

Differential expression analysis
Differential expression analysis for each cell type between different groups (YA and AA in cohort-1 and YH, AH, YCR and ACR in cohort-3) was performed using the t-test as implemented in the ‘‘sc.tl.rank_genes_groups’’ function of the Scanpy package. For each cluster, differentially-expressed genes (DEGs) were performed using the t-test and generated relative to all of the other cells. Before executing the differential expression analysis, we filtered out the cell types that were missing or had fewer than three cells in the comparison groups. An aging-associated and disease-related DEG dataset was established (adjusted P value < 0.05, |Log2FC| > 0.25) after identification of DEGs between AA and YA groups in cohort-1, AH and YH groups in cohort-3, ACR and AH groups in cohort-3, YCR and YH groups in cohort-3. The ‘‘upregulated DEGs during aging’’ were defined as the DEGs that increased in AA group and decreased in YA group. The ‘‘downregulated DEGs in aging’’ were defined as the DEGs that decreased in AA group and increased in YA group.

Gene functional annotation
The Metascape webtool (www.metascape.org) (Zhou et al., 2019) that allow visualization of functional patterns of gene clusters and statistical analysis was used to conduct DEGs gene ontology, pathway enrichment analyses. Among the top 30 enriched GO terms or pathways across various types of cells and tissues, 10 GO terms or pathways which were associated with aging were visualized. Gene expression profile cluster plots and heatmaps were established using the pheatmap R package (https://cran.r-project.org/web/packages/pheatmap/index.html, version 1.0.12).

Aging score analysis
To assess the impact of aging in circulating immune cells, we selected the top 20 genes out of 60 common upregulated genes in all immune cells. Aging scores were estimated for all cells as the average of the scaled (Z-normalized) expression of the genes in the list. The score of aging for all immune cell types can be used to predict the effect of aging on single cells and cell subtype levels. Calculation of the scores was done as follows: the score of the aging gene set in the given cell-subset (named as X) was computed as the sum of all UMI for all the aging genes expressed in X cell, divided by the sum of all UMI expressed by X cell (Pont et al., 2019).

Sequencing and analysis of TCR and BCR V(D)J
PCR amplification was done to enrich the full-length TCR/BCR V(D)J segments for the amplified cDNA from 5′ libraries with a Chromium Single-Cell V(D)J Enrichment kit (10 Genomics). The TCR/BCR sequences of each T/B cell were clustered using the CellRanger vdj pipeline (version 3.1.0, allowing identification of CDR3 sequence and the rearranged TCR/BCR gene. Analysis was performed using Loupe V(D)J Browser version 2.0.1 (https://support.10xgenomics.com, 10x Genomics). In summary, barcode information a containing clonotype frequency and TCR/BCR diversity metric were obtained. We projected T /B cells with dominant TCR/BCR clonotypes on a t-SNE plot using barcode information (Wen et al., 2020).

Determination of cell-cell interaction
We employed the expression of immune-related ligands and receptors to assess the cell-cell interactions (Ma et al., 2020). The possible ligand-receptor interactions between one set of receptor-expressing cells and then next ligand-expressing cells were determined as the average of the product of receptor and ligand expression (respectively from set one and two) across all single-cell pairs:
\documentclass[12pt]{minimal} 				\usepackage{amsmath} 				\usepackage{wasysym} 				\usepackage{amsfonts} 				\usepackage{amssymb} 				\usepackage{amsbsy} 				\usepackage{mathrsfs} 				\usepackage{upgreek} 				\setlength{\oddsidemargin}{-69pt} 				\begin{document}$$I = \mathop \sum \limits_{i}^{n} I_{i} \times \mathop \sum \limits_{j}^{m} r_{j} \left(\frac{1}{m*n}\right)$$\end{document}I=∑inIi×∑jmrj1m∗nwhere I refers to the interaction score between receptor expressing cells in set one and ligand-expressing cells in set two, Ii stands for the ligand expression of cell i in cell set one, rj represents the receptor expression of cell j in cell set two, n stands for the number of cells in set one and m denotes the number of cells in set two.
In the gene list, there were 168 pairs of well-annotated ligands and receptors, among which were co-stimulators, chemokines and cytokines. The possible interactions between two cell types were orchestrated by receptor-ligand pairs by the product of the average expression levels of the ligand in one cell type and the respective receptor in the other cell type.

Mass cytometry processing and quality control
CyTOF data were acquired with a CyTOF2 system using a SuperSampler fluidics system (Victorian Airships) at an event rate of < 400 events per second and normalized with Helios normalizer software (Fluidigm version 6.7.1016). Acquisitions from different days (three independent acquisitions were performed) were normalized using five-element beads (Fluidigm). Barcoded samples were deconvoluted and cross-sample doublets were filtered using cytobank application. CyTOF data was pre-processed with Cytobank (https://mtsinai.cytobank.org; Cytobank, 7.0) to sequentially remove calibration beads, dead cells, debris and barcodes for CD45+ PBMCs based on event length, DNA (191Ir and 193Ir) and live cell (195Pt) channels and then export the FCS files. We analyzed 200,000 PBMCs in cohort-1 and 160,000 PBMC in cohort-2, with an average of 20,000 cells per sample.

Mass cytometry visualizing and clustering
We created mass cytometry datasets for analysis by concatenating cells from all individuals for each cell type. In this way, we created downsampled datasets of 95,316 TCs, 35,254 NKs, 22,042 BCs, 39,144 MCs and 8,244 DCs in cohort-1 and 57,910 TCs, 34,857 NKs, 13,812 BCs, 45,431 MCs and 7,990 DCs in cohort-2 for analysis. We used FlowCore (65 flowCore: Basic structures for flow cytometry data.) to read and process FCS files for further analysis. For sample with more than 20,000 cells, we randomly selected 20,000 cells to ensure that samples were equally represented. At last, we run the t-SNE dimensionality reduction algorithm on a combined data sample using the Seurat package based on harmony embedding (https://github.com/immunogenomics/harmony, version 1.0.0).

Batch correction of mass cytometry data
PBMC mass cytometry data from 10 subjects of cohort-1 or 8 subjects of cohort-2 was combined and batch normalized using harmony respectively. First, mass cytometry data from each cohort all subjects was combined into a single dataset. Second, harmony batch correction was performed using one of the samples. Third, mass cytometry data were lognormalized in the Seurat’s NormalizeData function across the aggregated dataset.

Single-cell assay for transposase-accessible chromatin sequencing (scATAC-seq)
scATAC-seq targeting 4,000 cells per sample was performed using Chromium Single Cell ATAC Library and Gel Bead kit (10x Genomics, 1000110). Each sample library was uniquely barcoded and quantified by RT-qPCR. Libraries were then pooled and loaded on an Illumina Novaseq 6000 sequencer (3.5 pmol/L loading concentration, 50 + 8 + 16 + 49 bp read configuration) and sequenced to either 90% saturation or 30,000 unique reads per cell on average. All protocols to generate scATAC-seq data on the 10x Chromium platform, including sample preparation, library preparation and instrument and sequencing settings, are available here: https://support.10xgenomics.com/single-cell-atac.

scATAC-seq data processing

scATAC-seq processing
scATAC-seq reads were aligned to the GRCh38 (hg38) reference genome and quantified using CellRanger-ATAC count (10x Genomics, v.1.0.0).

scATAC-seq quality control
To ensure that each cell was both adequately sequenced and had a high signal-to-background ratio, we filtered cells with less than 1,000 unique fragments and enrichment at TSSs < 8. To calculate TSS enrichment > 2, genome-wide Tn5-corrected insertions were aggregated ± 2,000 bp relative (TSS-strand-corrected) to each unique TSS. This profile was normalized to the mean accessibility ± 1,900–2,000 bp from the TSS, smoothed every 51 bp and the maximum smoothed value was reported as TSS enrichment in R. To construct a counts matrix for each cell by each feature (peaks), we read each fragment.tsv.gz fill into a GenomicRanges object. For each Tn5 insertion, which can be thought of as the “start” and “end” of the ATAC fragments, we used findOverlaps to find all overlaps with the feature by insertions. Then we added a column with the unique id (integer) cell barcode to the overlaps object and fed this into a sparseMatrix in R. To calculate the fraction of reads/insertions in peaks, we used the colSums of the sparseMatrix and divided it by the number of insertions for each cell id barcode using table in R.

scATAC-seq visualization in genomic regions
To visualize scATAC-seq data, we read the fragments into a GenomicRanges object in R. We then computed sliding windows across each region we wanted to visualize for every 100 bp “slidingWindows (region, 100, 100)”. We computed a counts matrix for Tn5-corrected insertions as described above and then binarized this matrix. We then returned all non-zero indices (binarization) from the matrix (cell × 100-bp intervals) and plotted them in ggplot2 in R with “geom_tile”. For visualizing aggregate scATAC-seq data, the binarized matrix above was summed and normalized. Scale factors were computed by taking the binarized sum in the global peak set and normalizing to 10,000,000. Tracks were then plotted in Loupe Cell Browser, an interactive visualization software that shows scATAC-seq peak profiles for scATAC-seq cell clusters, similar to the analysis done in this manuscript and described at https://support.10xgenomics.com/single-cellatac/software/visualization/latest/what-is-loupe-cell-browser.

chromVAR
We measured global TF activity using chromVAR15. We used the cell-by-peaks and the Catalog of Inferred Sequence Binding Preferences (CIS-BP) motif (from chromVAR motifs “human_pwms_v1”) matches within these peaks from motifmatchr. We then computed the GC-bias-corrected deviation scores using the chromVAR “deviationScores” function.

Statistical analysis
The GraphPad Prism Software (version 8.0.2) was employed for data analysis and presentation. All results are presented as means ± SEM. Groups were compared with two-tailed Mann-Whitney-Wilcoxon tests and FDR was corrected using the Benjamini-Hochberg procedure. P value was derived by a hypergeometric test in representative GO terms and pathways.