Single‐cell RNA sequencing reveals characteristics of myeloid cells in post-acute sequelae of SARS-CoV-2 patients with persistent respiratory symptoms

Human venous peripheral blood samples were taken at the state's main tertiary care hospital-Ambulatory Post-COVID Clinic (Queen's Medical Center, Honolulu, Hawaii) and collected in ethylenediaminetetraacetic acid (EDTA) tubes (BD, Vacutainer) by venipuncture. In brief, venous blood was diluted with an equal volume of phosphate-buffered saline (PBS) and layered on top of Ficoll-Paque Plus (GE Healthcare Biosciences, Piscataway, NJ) following the manufacturer's protocol. PBMC were separated by centrifugation at 400 × g for 30 minutes at room temperature (RT). PBMCs were collected from the buffy coat, red blood cells were lysed, and the remaining pellet was washed twice in PBS supplemented with 2% fetal bovine serum (FBS). Cells were then counted, viability determined, and cryopreserved at 5 million cells/vial. One participant naïve to SARS-CoV-2 infection and two participants with prolonged COVID-19 pulmonary symptoms, confirmed via pulmonary function tests (PFTs), were selected for single-cell sequencing analysis. This study was approved by the Queen's Medical Center Research & Institutional Review Committee RA-2020-053 and by the University of Hawaii Institutional Review Board 2020-00406. Written informed consent was obtained from all participants and all assays were performed according to institutional guidelines and regulations.

PFTs was performed on individuals with PPASC. All PPASC participants underwent PFTs (Vyaire) with the measurements of forced vital capacity (FVC), forced expiratory volume in 1 second (FEV1), total lung capacity (TLC), and diffusion capacity corrected for hemoglobin percent predicted (DLCOc%) interpreted in accordance with American Thoracic Society (ATS) guidelines (58).

For live PBMC purification for scRNA-seq, PBMC were stained with BV711-CD45 (Clone H130, 1:200 dilution, Biolegend, San Diego, CA) for 30 minutes at RT after adding Human TruStain FcX (1:200 dilution, Biolegend, San Diego, CA) for 15 minutes. Propidium iodide (PI) was added just before cell acquisition to assess cell viability. UltraComp eBeads Compensation Beads (Thermo Fisher Scientific, 01-2222-42) were used for compensation. Approximately 98-99% of total PBMC were identified as live and CD45⁺ and were sorted into 1% bovine serum albumin (BSA) in DPBS by a FACSAria IIu Cell Sorter (BD Biosciences, Franklin Lakes, NJ). Sorted cell suspensions were pelleted at 350 x g for 5 minutes at 4°C, and cells were resuspended in 0.1% BSA in Dulbecco's Phosphate Buffered Saline (DPBS) to approximately 1,000 cells/μl.

Cell concentration and viability were confirmed using an automated cell counter Countess II (Thermo Fisher Scientific) with 0.4% trypan blue solution and samples with viability > 70% were further processed. To obtain single-cell gel beads-in-emulsion (GEM), cell suspension was pelleted by 400 x g for 5 minutes at 4°C, and cells were resuspended at a concentration of 1,000 cells/μl in 0.1% BSA in DPBS. scRNA-seq libraries were prepared using Chromium Next GEM Single-Cell 5' Reagent Kit (10x Genomics). Briefly, GEMs were generated in Chromium Controller by combining barcoded Single-Cell VDJ 5' Gel Beads v1.1, a Master Mix with mixture of single cells, and Partitioning Oil on Chromium Next GEM Chip G. To achieve single-cell resolution, cells were delivered at a limiting dilution, such that the majority (~90-99%) of generated GEMs contains no cell, while the remainder largely contain a single cell.

Immediately following GEM generation, the Gel Beads were dissolved and any co-partitioned cell was lysed. Oligonucleotides containing (i) an Illumina R1 sequence (read 1 primer sequence), (ii) a 16 nucleotide (nt) 10x Barcode, (iii) a 10 nt unique molecular identifier (UMI), and (iv) 13 nt template switch oligo (TSO) were released and mixed with the cell lysate and a Master Mix containing reverse transcription (RT) reagents and poly(dT) RT primers, resulting in full-length cDNA from poly-adenylated mRNA. GEMs were broken and pooled after GEM-RT reaction mixtures were recovered. Silane magnetic beads were used to purify the 10x Barcoded first-strand cDNA from the post-GEM-RT reaction mixture, which includes leftover biochemical reagents and primers. After cleanup, 10x Barcoded, full-length cDNA was amplified via PCR with primers (Forward primer [P5]: 5'-AATGATACGGCGACCACCGAGA-3' and Reverse primer [P7]: 5'-CAAGCAGAAGACGGCATACGAGAT-3'against common 5' and 3' ends added during GEM-RT). Amplification generated sufficient material to construct 5' Gene Expression libraries. Enzymatic fragmentation and size selection were used to optimize the cDNA amplicon size prior to the 5' Gene Expression library construction. P5, P7, a sample index, and Illumina R2 sequence (read 2 primer sequence) were added via End Repair, A-tailing, Adaptor Ligation, and Sample Index PCR. The final libraries contain the P5 and P7 priming sites used in Illumina sequencers. Library quality and size distribution were confirmed on Agilent Bioanalyzer High Sensitivity DNA chips (Agilent). Library quantification was performed with qPCR using the KAPA Library Quantification Kit for Illumina Platforms (KAPA/Roche). Libraries were normalized, denatured, and diluted to obtain 1.5 pM molarity. Libraries were sequenced on the NextSeq500 instrument at a depth of a minimum of 20,000 reads/cell using 26 x 91 bp read settings.

The scRNA-seq raw fastq files were aligned to the human reference genome (ver. GRch38) using the Cell Ranger pipeline (59) to generate a raw gene-by-cell count matrix. To integrate scRNA data generated from different batches, Harmony (60) was used to mitigate batch effects by merging and integrating the data. We chose a publicly available scRNA-seq dataset (GSM4509024) (61) that contained SARS-CoV-2 naïve and Severe COVID-19 PBMC samples with the study participants of similar age to our own samples. Our scRNA-seq data of PBMC (PPASC; n=2, Control; n=1) or PPASC data were combined with a naïve participant and then severe COVID-19 patients (n=4). DropletUtils was utilized to assess deviations from the ambient cell profile and empty droplets were identified (62), and cells with an FDR < 0.01 were considered statistically significant and retained as "real cells," while those with an false discovery rates (FDR) > 0.01 were deemed "empty droplets" and removed. The Scater package was used for the calculation of quality control (QC) metrics, visualized through principal component analysis (PCA) (63). A total of 34,139 cells passed quality control (QC) for the integrated data of PPASC and Controls. A total of 34,070 cells passed QC for the integrated data of PPASC and Severe (QC; Removing empty droplets, cells with low UMI/feature counts, high mitochondrial gene expression indicative of dying cells). Via these methods, 5,861 cells were removed (Supplementary Figures 1-4). The Seurat package's merge function was used to combine individual samples and remove cell-specific biases. The scran package's quickCluster function was utilized for grouping cells. Cell-specific size factors were calculated using the computeSumFactors function within scran. Raw counts of each cell were normalized by cell-specific size factor and log2-transformed with a pseudo-count of 1; highly variable genes (HVG) were defined based on the variance of the log expression profiles of each gene for decomposed technical and biological components by fitted mean-variance trends through the modelGeneVar function within scran. HVG were defined as those with FDR < 0.05. Downstream analysis involved computing 20 principal components (PCs), constructing a shared nearest neighbor (SNN) graph, and performing cell clustering using the Seurat package. Batch effects in the combined individual samples were addressed by applying 50 PCs and the RunHarmony function from the Harmony package. The detailed pipeline code is available on GitHub (https://github.com/lagom2728/PPASC-SARS-CoV-2-Paper-2023↗).

Cells were identified based on differential expressions of canonical marker genes and clustered accordingly (64). Additionally, cell types were confirmed by visual confirmation of the average expression values of canonical marker genes (64). Cell proportions were confirmed by generating cell type percentages, followed by a permutation test and bootstrapping to validate the statistically significant differences in proportions. As a result, cells were categorized into statistically significant groups based on an increase or decrease in the proportion of cells meeting the following criteria: relative differences in cell proportions for Log2 Fold Distribution (FD) in each group > 0.3, Log2FD < -0.3, and FDR < 0.05.

To assess the differential gene expression in the scRNA data of two groups (each PPASC and Control), we conducted differential expressed gene analysis with the built-in framework of MAST (65). In the analysis results, significant up- and down-regulated genes were classified based on statistical criteria of avg_log2FC > 0.25, < -0.25, and p-value < 0.05.

To conduct the analysis, Hallmark gene sets (H), Curated gene sets (C2), and Ontology gene sets (C5) databases of MsigDB (66) were selectively used and analyzed (67). Significant regulated pathways were categorized by their normalized enrichment score (NES) as enriched (>0) or depleted (<0), and p-values < 0.05.

The DoRothEA database (68) was used to identify transcription factors for differentially expressed target genes. The transcription factors from DoRothEA are complemented by an empirical measure of confidence, reflecting the certainty in their regulons. This confidence level is categorized into grades ranging from A (highest confidence) to E. In this particular benchmark, we specifically considered only the transcription factors (TFs) with confidence levels A and B, collectively referred to as DoRothEA (AB), used for transcription factor enrichment analysis.

TF enrichment analysis was performed using the msviper function within the viper package (69). Using the control patient as a reference, genes that were differentially expressed in each cell type via the Log2FD were identified and NES was calculated. TFs were classified as statistically significantly upregulated or downregulated according to the NES criteria and p-adjusted value thresholds (NES > 0, < 0, p-value < 0.05).

The human ligand-receptor (only Secreted signaling) database was used to analyze the ligand-receptor interaction between cells (70). Cells expressing less than 5% of genes were excluded, and only interactions with statistical significance (p-value < 0.05) among autocrine or paracrine interactions were classified. To compare the interaction signals between the two groups, the communication probability value was calculated based on the geometric mean expression level.

To make scRNA into pseudo-bulk RNA, we aggregated transcript count tables of all single-cells per sample to verify possible false positives in single-cell differentially expressed genes (DEGs). Differential expression analysis was performed (71), and genes satisfying LogFC > 0.25 or < -0.25, along with a p-value < 0.05 were classified as up- or down-regulated genes, respectively.

Integrated scRNA data was utilized to calculate Pearson correlation values (72). Subsequently, correlation networks were created for specific cell types, and the identification and validation of differentially expressed network modules were accomplished through enrichment analysis (73).

For metabolic analysis, we utilized the pipeline provided by the Compass package (74), specifically employing Single-Cell Flux Estimation Analysis and Single-Cell Flux Balance Analysis. These analytical tools were used to estimate cell-specific metabolic flux to infer the metabolic states of immune cells. For analysis of metabolite communication, we employed MEBOCOST (67). This model estimates the relative abundance of metabolic products based on the expression of genes encoding metabolic reaction enzymes. It additionally identifies cell-to-cell communication of metabolites and sensors by gathering enzyme-related genes from the Human Metabolome Database (75) across different cells.

Circulating immune cell levels and their functional status hold promise as biomarkers for assessing the severity of COVID-19. Analysis of blood immune cells in COVID-19 convalescents has been used to provide insight into the long-term consequences of host immune responses after SARS-CoV-2 infection. To explore immune system alteration associated with persistent lung sequelae after recovery from acute illness of SARS-CoV-2, we selected two participants who had reduced DLCOc% (<80%) by PFT among individuals with PPASC. They had persistent pulmonary symptoms greater than 5 months after onset of SARS-CoV-2 infection, were hospitalized for COVID-19 disease, and were fully vaccinated against SARS-CoV-2 (Supplementary Table 1). PASC symptoms included shortness of breath, chest pain, joint pain, brain fog, and depression. As a comparator, an age-matched participant naïve to SARS-CoV-2 infection (Control) was selected and confirmed negative via a SARS-CoV-2 nucleocapsid antibody test (Supplementary Table 1). Live CD45⁺ cells from two PPASC patients and this SARS-CoV-2 naïve individual were subjected to scRNAseq analysis utilizing the 10X Genomics Chromium system which yielded an integrated object of 34,139 cells after QC.

We integrated our scRNA-seq data with a naïve participant from a public scRNA-seq dataset (GSM4509024) (61) to compare the transcriptome of PPASC (n=2) and controls (n=2). Each of the four independent scRNA-seq samples was subjected to uniform manifold approximation and projection (UMAP) based on RNA expression (Figure 1A). A total of 11 cell types were identified, including CD4⁺/CD8⁺ T cells, natural killer (NK) and natural killer T (NKT) cells, CD14⁺ and CD16⁺ monocytes, dendritic cells, platelets, hematopoietic stem cells (HSC), B cells, and erythrocytes (Figure 1A). Canonical marker genes of known cell types were enriched in each cluster and the marker expressions were used to annotate the clusters (Figure 1B). All identified immune cell populations were present in both SARS-CoV-2-infected and naïve participants (Supplementary Figure 5), albeit in differing proportions between PPASC and controls (Figure 1C). We further attempted to determine a vital feature for reflecting immune cell alterations associated with PPASC via the relative proportions of peripheral immune cells. We observed significant increases in myeloid cells (CD14⁺/CD16⁺ monocytes and dendritic cells), and NKT cells and B cells, whereas the proportions of NK cells, HSC, and erythrocytes were significantly decreased in PPASC compared to controls (Figure 1D).

To explore the characteristics of immune cells and identify the molecular changes associated with PPASC, we performed a detailed analysis of the DEGs of immune cells from PPASC compared with those from the controls. Our scRNA-seq analysis identified pronounced alterations in myeloid cells. Furthermore, we found that elevated circulating monocyte levels and their activation were observed in COVID-19 convalescents (76). Thus, we focused our attention on myeloid cells, particularly CD14⁺ and CD16⁺ monocytes, as well as dendritic cells, for downstream analysis hereafter referred to as myeloid-lineage cells (MLCs). We found a total of 872 DEGs in monocytes and dendritic cells: CD14⁺ monocytes (228 upregulated and 132 downregulated), CD16⁺ monocytes (303 upregulated and 226 downregulated), and dendritic cells (332 upregulated and 217 downregulated) (Supplementary Figures 6A, B). In addition, we demonstrated shared DEGs across three cell types (118 upregulated and 54 downregulated) (Supplementary Figures 6A, B). Interestingly, our data showed that genes known to be associated with pulmonary fibrosis (ANKRD11, CTNNB1, CXCR4, HIF1A, HMGB1, ITSN2, LITAF, NEAT1, VEGFA, and DSE) (77 –84) were upregulated in MLCs. However, genes associated with glycolytic metabolism (ALDOA, PGK1, TPI1, and MYL6) were downregulated in these cell populations (Figures 2A, B).

Next, we performed gene set enrichment analysis (GSEA) across each cell type and conducted direct comparisons of the PPASC versus the control groups, respectively. Overall, compared to controls, PPASC displayed upregulated pathways related to pulmonary fibrosis, such as VEGF, WNT, and TGF-β and cell apoptosis, while pathways related to cytokine/inflammation (including interferon) were downregulated in MLCs (Figure 2C). The VEGF pathway was upregulated in CD14⁺ monocytes and dendritic cells, while the WNT pathway was upregulated only in CD16⁺ monocytes (Figure 2C).

To identify potential direct targets involved in the induction of profibrotic features in myeloid cells, we performed transcription factor enrichment analysis (TFEA) to identify transcription factors inferred from the perturbation of gene signatures between PPASC and controls. We found that a total of 27 transcription factors were significantly dysregulated in individual cell types between groups (Figure 2D). Consistent with earlier analyses, we found that the downstream transcription factors of VEGFR (HIF1A; CD14⁺ monocytes) and TGFBR (SMAD and STAT3; CD14⁺ monocytes and dendritic cells, ATF2; CD14⁺/CD16⁺ monocytes and dendritic cells) signaling were upregulated in MLCs (Figure 2D). These observations demonstrate that PPASC display enriched genes driving pulmonary symptoms and profibrotic features in MLCs and that increased MLC proportions and the altered gene signatures likely contribute to development of pulmonary sequelae after SARS-CoV-2 infection.

We were intrigued as to how our findings in PPASC patients, who had recovered from acute COVID-19 disease and suffered ongoing symptoms, compared with patients who suffered severe COVID-19. We integrated a publicly available scRNA-seq dataset from severe COVID-19 patients, hereafter referred to as severe COVID-19 (n=4) (Supplementary Figures 7A-F) and observed 10 of the previously identified cell types via UMAP visualization of the combined dataset; however, HSCs were unidentifiable (Supplementary Figures 8A, B). The cell proportion makeup differed drastically between PPASC and severe COVID-19, with CD16⁺ monocytes and dendritic cells found to be significantly increased and CD14⁺ monocytes significantly decreased in PPASC compared to severe COVID-19 (Supplementary Figures 8C, D). MLCs demonstrated an up-regulation of fibrosis-associated genes, such as VEGFA and HIF1A in PPASC compared to severe COVID-19 (Supplementary Figure 9A). However, genes involved in glycolytic metabolism were not found to be downregulated in PPASC, compared to severe COVID-19 (Supplementary Figure 9B). Pulmonary fibrosis-related pathways and cell apoptosis and stress-related pathways were found to be upregulated in PPASC compared to severe MLCs (Supplementary Figure 9C). Immune response pathways were found to be downregulated in CD14⁺ monocytes but not in CD16⁺ monocytes or dendritic cells in PPASC compared to severe COVID-19 (Supplementary Figure 9C). TF enrichment analysis revealed upregulation of fibrotic-associated ATF2, STAT3, and SMAD family TFs across MLCs in PPASC compared to severe COVID-19. STAT3 was specifically upregulated in CD16⁺ monocytes, while ATF2 was upregulated in CD14⁺ monocytes (Supplementary Figure 9D). DEGs in MLCs of PPASC compared to the severe COVID-19 demonstrated a larger overlapping proportion of up-regulated genes (Supplementary Figure 10A). Few overlapping downregulated DEGs were noted in MLC populations (Supplementary Figure 10B) and were not inclusive of glycolytic genes as found in PPASC vs control MLCs (Data not shown). Specifically, via analysis of cell-to-cell interactions in MLCs, we found that VEGF-VEGFR interaction was only found in PPASC but not in severe COVID-19 (Supplementary Figure 11A). VEGFA and its ligand expression was minimal in CD14⁺ and CD16⁺ monocytes in the severe COVID-19 but was over 40% in dendritic cells (Supplementary Figure 11B). Pseudo-bulk DEG revealed no differential expression of VEGF in CD16⁺ monocytes but confirmed up-regulation of VEGF in CD14⁺ monocytes and dendritic cells in PPASC (Supplementary Figure 11C). These data suggest that VEGF expression patterns in MLCs clearly distinguished PPASC from severe COVID-19.

To gain a comprehensive view of how the inferred targets may interact regarding PPASC development, we then identified potential cell–cell interactions and differential ligand-receptor interactions that are conserved in PPASC. The pathway of VEGF interaction (VEGFA-FLT1/KDR) was increased among MLCs, compared to controls (Figure 3A). We observed high-fidelity communications via VEGF pathways in MLCs via an autocrine or paracrine manner (Figure 3A). Furthermore, using pseudo-bulk single-cell data, we validated the expression of genes identified from ligand-receptor interaction analysis. VEGFA percent expression was increased in MLCs from PPASC compared to those cells from controls (Figure 3B). Notably, VEGFA was highly expressed in CD14⁺ and CD16⁺ monocytes, as well as dendritic cells (Figures 3B, C).

Within the monocyte clusters (CD14⁺ and CD16⁺ monocytes), we further sub-grouped the signatures into gene correlation based on network modules (Figure 4A). Within CD14⁺ monocyte populations, network module 5 (Net-M5) was specifically found to be upregulated, while CD16⁺ monocyte demonstrated an upregulation of network module 4 (Net-M4) (Figure 4B). To identify the organizing hub genes in Net-M4 and M5, the top 60 genes with high eigengene-based connectivity (kME) values were selected (Supplementary Figure 12). In the Net-M4, genes involved in VEGFA-VEGFR2 signaling (NR4A1, TPM3, FAM120A, and NUMB) and WNT signaling (HHEX and APC) were identified. The Net-M5 detected genes involved in lung fibrosis (CXCL8 and IL1B), VEGFA-VEGFR2 signaling (NR4A1 and CBL), and TGF-β signaling (JUNB and ATF3) (Supplementary Figure 7). Furthermore, the pathway enrichment analysis of Net-M4 and Net-M5 revealed that the CD16⁺ M4 module was associated with WNT, VEGF, and TGF-β-signaling pathways (Figure 4C). The CD14⁺ M5 module was associated with COVID-19 adverse outcome linked pathways as well as fibroblast growth factor stimulation pathways. Interestingly, this same module in CD14⁺ monocytes revealed increases in pathways for fibroblast-specific apoptosis and general apoptotic pathways (Figure 4C).

Reports of the metabolic impacts of SARS-CoV-2 infection and association with chronic fatigue and diabetic-like events during Long-COVID are well documented (85). These observations aligned with our findings of downregulated genes associated with glycolytic metabolism (Figure 2B) and prompted us to determine the impact of PASC on immunometabolism in individuals with persistent pulmonary sequelae. Using Gene Set Variation Analysis (GSVA) based on non-parametric unsupervised learning to infer the variability of gene set enrichment related to gluconeogenesis/glycolysis, we observed that the average score of glycolysis-related pathways in MLC populations was lower in PPASC than in controls (Figure 5A). The inferred metabolic state of individual cells using a database of integrated metabolic networks and metabolic flux balance analysis demonstrated that MLCs of PPASC exhibited downregulation of gluconeogenesis/glycolysis metabolic reactions compared to the control group (Figure 5B). Furthermore, inference of intercellular metabolite communication based on estimated metabolite scores and averaged sensor gene expression values by cell group demonstrated that the metabolite, D-Glucose, was directly associated with gluconeogenesis/glycolysis metabolism. Furthermore, these inferred interaction scores revealed decreased mean abundance values in MLCs of PPASC compared to the control group (Figure 5C). The calculated communication scores for metabolite-sensor partners revealed statistically significant interaction scores in MLCs of the control group but not in PPASC (Figure 5D). This suggests a lack of gluconeogenesis/glycolysis-related metabolite interaction or a diminished interaction in PPASC. In summary, these metabolic analysis observations suggest a decreased regulation of glucose-specific metabolic pathways in MLCs of PPASC.