Differential chromatin accessibility in peripheral blood mononuclear cells underlies COVID-19 disease severity prior to seroconversion

SARS-CoV-2 infection triggers profound and variable immune responses in human hosts. Chromatin remodeling has been observed in individuals severely ill or convalescing with COVID-19, but chromatin remodeling early in disease prior to anti-spike protein IgG seroconversion has not been defined. We performed the Assay for Transposase-Accessible Chromatin using sequencing (ATAC-seq) and RNA-seq on peripheral blood mononuclear cells (PBMCs) from outpatients with mild or moderate symptom severity at different stages of clinical illness. Early in the disease course prior to IgG seroconversion, modifications in chromatin accessibility associated with mild or moderate symptoms were already robust and included severity-associated changes in accessibility of genes in interleukin signaling, regulation of cell differentiation and cell morphology. Furthermore, single-cell analyses revealed evolution of the chromatin accessibility landscape and transcription factor motif accessibility for individual PBMC cell types over time. The most extensive remodeling occurred in CD14+ monocytes, where sub-populations with distinct chromatin accessibility profiles were observed prior to seroconversion. Mild symptom severity was marked by upregulation of classical antiviral pathways, including those regulating IRF1 and IRF7, whereas in moderate disease, these classical antiviral signals diminished, suggesting dysregulated and less effective responses. Together, these observations offer novel insight into the epigenome of early mild SARS-CoV-2 infection and suggest that detection of chromatin remodeling in early disease may offer promise for a new class of diagnostic tools for COVID-19.


Results
A clinical cohort to study early COVID-19 with mild or moderate symptoms. To establish a cohort for the longitudinal study of early SARS-CoV-2 infection, subjects were selected from a large ongoing community-based prospective cohort, the Molecular Epidemiological Study of Suspected Infection (MESSI) in four experimental groups: (1) pre-pandemic healthy controls, (2) subjects exposed to SARS-CoV-2 as close contacts but negative for SARS-CoV-2 by qPCR and serology, (3) outpatients positive for SARS-CoV-2 by qPCR whose maximal illness was mild and (4) outpatients positive for SARS-CoV-2 by qPCR whose maximal illness was moderate. Demographic characteristics of study participants are shown in Table 1. Symptoms associated with COVID-19 illness were self-reported to clinic personnel using a COVID-19 symptom survey developed by Duke University Medical Center consisting of 38 symptoms. Each subject reported a symptom and perceived severity of that symptom from 0 (none) to 4 (very severe). Subjects with mild symptoms exhibited a mean score of 12.8 ± 1.9 and 12.2 ± 2.8 for total PBMC assays (bulk) and assays of individual PBMCs (single-cell assays), respectively, which corresponded to the World Health Organization (WHO) Ordinal 8-point scale 1 (OS1, ambulatory and with no limitation of activities). Subjects with moderate symptoms exhibited a mean score of 33.6 ± 2.4 and 36.0 ± 2.3 respectively which corresponded to the WHO Ordinal 8-point scale 2, (OS2, ambulatory and with limitation of activities). Interestingly, despite lacking any evidence of infection (negative SARS-CoV-2 qPCR and serological testing for at least two months after exposure), close household contacts who were exposed to SARS-CoV-2 but remained negative for infection showed broadly similar quantitative symptom scores to subjects with mild disease, with a mean score of 9.9 ± 2.9 and 3.3 ± 2.8 respectively ( Table 1). The sum of symptom severity across all symptoms reported was able to distinguish moderate subjects from mild or close contacts but did not distinguish mild from close contacts (Fig. S1A,B). However, the specific symptoms experienced, and their severity, were different between all three groups (Fig. S1C,D). Loss of taste and smell, headache, malaise, and fatigue were most common in the moderate group. Coughing was most common in the mild group and runny nose was most common in the close contacts.
Transcriptional profiles of PBMCs distinguish mild or moderate COVID-19 from healthy controls. To identify underlying molecular mechanisms of disease severity in mild or moderate COVID-19 both www.nature.com/scientificreports/ prior and after seroconversion, PBMCs were isolated and profiled as total PBMCs (bulk assays) or as single nuclei isolated from individual PBMCs (single cell assays). Longitudinal sample collection from each subject spanned the transition from seronegative (IgG−) to seropositive (IgG+) for SARS-CoV-2 spike protein in the mild and moderate symptom groups (Supplementary Table S1 and S2). For all pooled analysis of IgG− or IgG+  subjects, the dataset corresponding to the last IgG− timepoint or the first IgG+ timepoint, respectively, was used.  Differential gene expression analysis of bulk RNA-seq data from mild IgG− COVID-19 subjects compared to  healthy controls identified a transcriptional signature of 45 genes with p ≤ 0.05; moderate IgG− COVID-19 subjects compared to healthy controls identified 71 genes with p ≤ 0.05 (Fig. 1A). Using principal component analysis (PCA), differential gene expression between healthy controls and subjects with mild symptoms was sufficient to distinguish the groups Fig. 1B). Similar comparisons using RNA-seq from mild IgG+ subjects compared to healthy controls or moderate IgG+ subjects compared to healthy controls identified 250 and 43 genes, respectively, with p ≤ 0.05 (Fig. 1C). Using PCA, differential gene expression between healthy controls and subjects with moderate symptoms (Fig. 1D) was sufficient to distinguish the groups. However, both mild and moderate subjects clustered together and could not be easily distinguished. The gene symbols, gene names, fold change and p values adjusted for multiple hypothesis testing for the four comparisons are found in Supplemental Table S3.
Considering the level of gene expression (log fold change, LFC) in mild and moderate subjects compared to healthy controls, differentially expressed genes were identified in only mild, only moderate or both symptom cohorts (Fig. 1E). Comparing IgG− subjects, we identified 30 genes observed in both mild and moderate cohorts.
In contrast, we did not identify any genes solely associated with moderate IgG+ subjects and identified a further 43 genes observed in both IgG+ mild and moderate cohorts. Healthy controls exhibited higher expression of the interleukin-8 precursor CXCL8 compared to other groups. Conversely, the interleukin-8 receptor, CXCR1, had higher expression in mild and moderate IgG− subjects. Lower expression of chemokine CCL3 and cytokine IL1B 12 , and higher expression of the myeloid cell plasticity regulator KLF6 13 were observed in both IgG− and IgG+ COVID-19 subjects.
Chromatin remodeling prior to seroconversion can distinguish mild from moderate COVID-19. To test whether chromatin accessibility can distinguish COVID-19 disease severity in individuals early in their disease course (prior to seroconversion) an Assay for Transposase-Accessible Chromatin using sequencing (ATAC-seq) 14 was used to quantify differential chromatin accessibility in three groups: subjects with mild symptoms, subjects with moderate symptoms and pre-pandemic healthy controls. A set of 455 differentially accessible chromatin regions (DARs) was found to distinguish IgG− subjects with mild symptoms from those with moderate symptoms at p ≤ 0.05 ( Fig. 2A). Of these markers, 73 regions shared differential accessibility in mild and moderate IgG− subjects, demonstrating both severity-associated evidence of chromatin remodeling and a smaller conserved epigenetic response. Additionally, a set of 375 regions of differentially accessible chromatin was identified in IgG+ subjects (Fig. 2B). The effect sizes of each DAR were compared across healthy controls, mild, and moderate IgG− subjects to determine the overlap features associated with either cohort. We determined that many of the markers that distinguish mild subjects from healthy controls also distinguished moderate subjects from healthy controls with little difference in effect size. In contrast, the markers that distinguished mild and moderate subjects from each other had higher accessibility in healthy controls and the lowest accessibility in moderate subjects. DARs specific to the comparisons in 2A-B as heat map depictions (Fig. 2C,D) for prior and post seroconversion shows overlap between mild and moderate symptom subjects compared to healthy controls. Functional enrichment was then performed for each set of differentially accessible regions to identify an association with transcription factor motifs and downstream pathways. In IgG− subjects with mild symptoms, we found an enrichment of transcription factors and pathways related to viral infection (Fig. 2E). Fewer enriched pathways were observed in the IgG− subjects with moderate symptoms. Transcription factors found to be enriched in accessible chromatin from IgG+ subjects with mild symptoms were associated with dysregulation of myeloid development via HOXA3, C/EBP-B and C/EBP-D 15 (Fig. 2F).
CD14+ monocytes are the primary cell type differentially activated in IgG− subjects. The chromatin accessibility signature detected in the bulk ATAC-seq datasets identified candidate epigenetic biomarkers associated with COVID-19 severity and showed that PBMCs undergo extensive chromatin remodeling in response to SARS-CoV-2 infection. To understand how each cell type contributed to this signature and to track the evolution of the chromatin landscape of each cell type during seroconversion, we performed single-cell (sc) ATAC-seq profiling of the PBMCs isolated from subjects with mild or moderate symptoms and healthy controls ( Table 1 and Supplemental Table S2). The relative abundance of cell types collected from each subject cohort were similar between IgG− and IgG+ timepoints (  Table S6). Analysis of transcription factor motif accessibility was applied to PBMCs from IgG− mild or moderate subjects to identify shared regulatory mechanisms. Transcription factors with highly accessible motifs and low gene expression were identified, consistent with epigenetically primed chromatin remodeling without corresponding transcriptional activity (Fig. 3B). These motifs, including those of the CREB and ETS protein families, were characterized by increased binding activity as measured by footprint depth and relatively low average per-cell gene expression. Inflammatory regulators including JUN/ FOS, proteins from the KLF family, and IRF1 showed elevated gene expression. A complete list of transcription factors associated with accessible regions and low gene expression is found in Supplemental Table S7.
Comparison of PBMCs collected from subjects with mild or moderate symptoms identified specific peaks enriched on the day of study enrollment (IgG−) versus a timepoint 2-4 weeks later (IgG+) for all cell types (Fig. 3C). A total of 5988 peaks were identified as uniquely accessible at one of the three longitudinal (IgG−, www.nature.com/scientificreports/ mixed serology, and IgG+) timepoints. A similar experimental approach was applied to each major cell type individually to determine the extent of chromatin remodeling prior to seroconversion. The cells with enrichment of significant peaks at the early, IgG− timepoint were primarily CD14+ monocytes (Fig. 3D).   were identified in all cell types collected from IgG− subjects with mild or moderate symptoms by counting the number of peak-to-gene linkages for each gene in the scATAC-seq datasets 7 . A total of 1109 genes with greater than 10 peak-to-gene linkages were defined as DORC genes, and those known to be regulated by a super-enhancer were labeled (Fig. 4A). These DORCs showed cell type-specific profiles and appeared to have higher accessibility than gene expression in these cells, suggesting that these genes were being primed for activation (Fig. 4B). Accessibility at the DORCs regulated by a super-enhancer was enriched in CD14+ monocytes and other myeloid cells, including dendritic cells and CD16+ monocytes, indicated by the presence of proximal gene regulatory elements with correlated accessibility to DORC gene loci (Fig. 4C). These DORC genes regulated by super-enhancers play a role in priming active chromatin states in CD14+ monocytes, consistent with the cell fate decisions identified that distinguished subjects with variable symptom severity.
Epigenetic profiles of CD14+ monocyte subpopulations differentiate mild from moderate IgG− COVID-19. We next identified transcriptomic and epigenetic profiles of sub-populations of CD14+ monocytes that distinguished IgG− subjects with mild from moderate symptoms. Single-cell gene expression data for CD14+ monocytes were processed using the uniform manifold and projection (UMAP) algorithm. Clusters 2 and 4 identified a unique sub-population of CD14+ monocytes in subjects with mild symptoms and cluster 7 identified a sub-population of CD14+ monocytes in subjects with moderate symptoms (Fig. 5A). Transcription factors and downstream targets were identified for each cluster by measuring correlated gene expression www.nature.com/scientificreports/ with promoter-proximal transcription factor motifs (Supplemental Methods). CD14+ monocyte subpopulations from subjects with mild symptoms exhibited significant up-regulation of the regions that regulate IRF7, IRF1, and STAT1 transcription factors, whereas the CD14+ monocyte subpopulations from subjects with moderate symptoms were characterized by up-regulation of the CEBPB and KLF3 regulatory regions (Fig. 5B). Single-cell ATAC-seq data for CD14+ monocytes were processed using a similar approach to identify clusters with correlated gene expression regulatory activity (Fig. 5C). Activity at the transcription factor motifs was estimated by measuring genome-wide motif accessibility, and the correlations between scRNA-seq and scATAC-seq clusters were calculated in a pairwise manner (Fig. 5D). Differential gene expression between the mild and moderate sub-populations identified interferon-stimulated gene activation as a distinguishing feature of symptom severity. In addition to these transcriptomic markers, gene expression measurements of enriched monocyte-related pathways were calculated to quantify the average transcriptomic response in each cell subset (Fig. 5F, Supplemental Tables S8, S9). To perform this analysis, scRNA-seq clusters containing CD14+ monocytes from both mild and moderate subjects were merged were appropriate (Supplemental Methods) and two published datasets were incorporated to enable comparison with more severe COVID-19 disease phenotypes 16,17 . We identified an enrichment of interferon response activity in the mild and moderate monocytes from Mild 1, Mild 2 and Moderate 1 CD14+ monocyte subsets (scRNA-seq clusters 2, 4 and 7) which was consistent with our bulk RNA-seq data on the same cohorts. CD14+ monocytes from subjects with severe or critical disease did not share the same strong gene expression signature of interferon signaling. Furthermore, monocytes from these clusters have elevated expression of TNF-a, NF-kB, toll-like receptor (TLR) signaling, and monocytic recruitment factor CCR2 which helps explain why these cells cluster separately in Fig. 5A. Expression of these pathways is also elevated in monocyte subsets from subjects with more severe disease, including the ICU subjects treated with TLR agonists 16 . Pathways identified centered on those of the adaptive immune response. Differential regulation of pathways related to the transition from innate to adaptive immunity, marked by dectin-1 and the toll-like receptor cascade, suggested that activation of T and B cell receptor signaling were observed in these sub-populations [18][19][20] .

Discussion
This study demonstrated that the epigenomes of critical subsets of PBMCs were remodeled extensively early in SARS-CoV-2 infection and reflected disease severity in IgG− subjects with mild and moderate symptoms. Specifically, differential activity of transcription factors and chromatin accessibility prior to anti-spike IgG seroconversion were observed to distinguish disease severity, preceding the later transcriptional response. The fact that these signals were robust even in a small cohort of subjects with subtle differences in disease severity (WHO OS 1 vs 2) suggests the potential power of epigenetic approaches to classify subsets of SARS-CoV-2 infection. This study complements previous work that identified signatures of accessible chromatin in severe or convalescing COVID-19 [21][22][23][24] .
We found that transcription factor motifs enriched in peaks accessible prior to seroconversion had the highest occupancy in cells from the myeloid lineage, specifically CD14+ monocytes, dendritic cells, and CD16+ monocytes. KLF and CREB transcription factor families, which are known to regulate monocyte-macrophage polarization, had elevated motif accessibility in subjects with mild or moderate symptoms, suggesting an epigenetic priming mechanism controlling cell fate 25,26 . This distinction is important as not only are CD14+ monocytes generally more activated in early COVID-19, as has been described, but here were differentially activated across subjects with variable disease severity. These differences included modulation of IFN-γ, as reported previously. Specifically, mild symptom severity (or effective control) was marked by upregulation of classical antiviral pathways including those regulating IRF1 and IRF7. With increased severity, these antiviral signals diminished, suggesting that dysregulated and less effective responses underlied moderate disease. These early epigenetic changes occurred prior to transcriptional manifestations and in this cohort were also more stable and predictive of disease severity than gene expression alone. These findings suggest that epigenetic approaches focusing on chromatin accessibility may offer even greater promise than other host-based molecular diagnostics and predictive tools, as has been suggested for DNA methylation-based profiles 27 . Our work is consistent both with previous studies that identified DNA-methylation signatures, which included viral response and interferon signaling that will predicted SARS-CoV-2 infection and clinical outcome 27 , and with genome-wide DNA methylation signatures associated with severe COVID-19 that highlight hypermethylation in IFN-related genes, hypomethylation in inflammatory genes and increased epigenetic age 28 .
Interestingly, CD14+ monocytes underwent the most extensive chromatin remodeling over time and exhibited epigenetic profiles at early seronegative times that distinguished disease severity. In contrast, after seroconversion, motifs enriched in accessible peaks were not seen in monocytes, instead showing the highest activity in B cells and plasmablasts, consistent with the transition from inflammatory signaling to adaptive immune development in this phase. Our work is consistent with epigenetic profiles of immune cells of individuals convalescing from COVID-19 that show establishment of immunological memory 24 .
There are limitations to this study, including the relatively small sample size, nonstandard enrollment and sampling times that vary across individuals, and a lack of subjects with critical illness to examine how these epigenetic changes manifest in more severe disease. However, development of a host-response assay that leverage the highly sensitive epigenetic biomarkers established early during infection has the potential to fill a clear unmet clinical need in the care of patients with COVID-19 29,30 . For example, in addition to the well characterized respiratory system damage of COVID-19, the deleterious effects on the central nervous system (CNS) can be devastating and include headache, anosmia (loss of smell), hyposmia (loss of taste), disturbance of smell, taste or vision, epileptic seizures, Guillain-Barre syndrome and intracerebral hemorrhage 31 . Our future studies will determine association between epigenetic biomarkers and specific symptoms that reflect damage to the CNS or peripheral nervous systems. www.nature.com/scientificreports/ In summary, we found that the evolution of the chromatin landscape in circulating leukocytes during COVID-19 primes host immunological responses at early times, is mediated primarily by CD14+ monocytes and correlates with an observed divergence in disease severity. These changes temporally precede transcriptional manifestations of pathways related to adaptive immunity. Together, these observations offer novel insight into severity-associated variations in host responses to SARS-CoV-2 infection and suggest that detection of critical components of chromatin remodeling in early disease may offer promise for a new class of diagnostic tools for COVID-19. The COVID-19 pandemic has highlighted with excruciating clarity the need for prognostic biomarkers to triage patients and personalize their treatment relevant to their likelihood to experience disease progression. Our findings offer a potential pathway via translation to diagnostic platforms capable of detecting markers of chromatin accessibility. Furthermore, improved understanding of the transcriptional and epigenetic response may reveal novel therapeutic opportunities. We envision that the rapid development of nanotechnology that supported mRNA vaccines can be harnessed to support detection of peripheral blood-based epigenetic biomarkers of early pre-seroconversion COVID-19 32 .

Methods
Cohort recruitment, biological sample collection, and initial phenotyping. The study was approved by the Duke University Institutional Review Board. Protection of human subjects was in accordance with research protocols approved by the Duke University Institutional Review Board, consistent with the Declaration of Helsinki. Written informed consent was obtained from all research subjects or their legally authorized representatives. Subjects with confirmed or suspected SARS-CoV-2 infection or their close contacts were identified in the outpatient setting and enrolled into the Molecular and Epidemiological Study of Suspected Infection protocol (MESSI, IRB Pro00100241). All close contacts and subjects with mild or moderate COVID-19 were longitudinally sampled from enrollment to convalescent phase. Biological samples and demographic information were collected prospectively at first visit (Day 0) and at weekly intervals on Day 7 and Day 14. At each visit, infection with SARS-CoV-2 was confirmed using quantitative polymerase chain reaction (qPCR) of nasopharyngeal (NP) swab samples, and serology testing was performed for IgG against the SARS-CoV-2 spike domain. All subjects with mild or moderate COVID-19 progressed from seronegative (IgG−) to seropositive (IgG+). Close contacts were qPCR negative and IgG− at all time points; healthy controls were enrolled pre-pandemic and were not tested for SARS-CoV-2 or spike protein IgG. Self-reported symptom surveys were performed at each visit. To quantify symptom severity, the sum of 38 defined symptoms, each scored 0-4 (0-none, 1-mild, 2-moderate, 3-severe, 4-very severe), was determined from symptom onset through each longitudinal sample collection. SARS-CoV-2 q-PCR tests used virus RNA extracted from NP samples in 140 µL of viral transport medium (VTM) using QIAamp Viral RNA Mini Kit (QIAGEN, Cat# 52904) according to manufacturer's instructions. SARS-CoV-2 nucleocapsid (N1) and human RNase P (RPP30) RNA copies were determined using 5 µL of isolated RNA in the CDC-designed kit (CDC-006-00019, Revision: 03, Integrated DNA Technologies 2019-nCoV kit). SARS-CoV-2 IgG ELISA tests for antibody response were performed using the anti-SARS-CoV-2 spike S1 domain IgG ELISA assay (EUROIMMUN Medizinische Labordiagnostika AG, Lübeck, Germany).

Purification of PBMCs.
Peripheral blood mononuclear cells (PBMCs) were purified using the Ficoll-Hypaque density gradient method according to manufacturer's instructions. Briefly, whole blood was collected in ACD Vacutainer tubes (Becton Dickinson, Cat# 364606) and processed within 8 h by dilution 1:2 in PBS, layered onto the Ficoll-Hypaque (Sigma Aldrich, Cat# GE17-1440-02) in 50 ml conical tubes, and centrifuged at 420 × g for 25 min. Buffy coat was collected, washed twice in D-PBS (Sigma Aldrich, Cat# D8537) by centrifugation at 400 × g for 10 min to isolate PBMCs which were assessed for viability and cell count using a Vi-Cell automated cell counter (Beckman-Coulter). PBMCs were adjusted to 10 × 10 6 cells/ml in cryopreservation media (90% FBS, 10% DMSO), frozen at − 80 °C using CoolCell LX (BioCision) for 12-24 h and stored in liquid nitrogen vapor phase.
RNA extraction, total RNA-seq, and analysis. RNA was extracted from 300 K cells using the Zymo Direct-zol RNA Miniprep Kit (Zymo Research, Cat# R2051) and RNA quality assessed using the Agilent DNA Screentape assay. The RNA Integrity Number (RIN) scores for all samples were > 7.0. Total RNA libraries were generated using the NuGEN Ovation® SoLo RNA-Seq Library Preparation Kit (Tecan Life Sciences, Cat# 0500-96). Libraries were sequenced using Illumina NovaSeq 6000 instrument with S4 flow cell and 150 base pair paired-end reads (Illumina, Cat# 20012866). FASTQ files were generated from the NovaSeq BCL outputs and quality was assessed with FASTQC 33 . Differentially expressed genes were identified between subjects with different disease severity using the limma package and voom to model variance 34 . Nuclei purification, ATAC-seq and analysis. Nuclei were extracted from frozen PBMCs. Briefly, 100 K cells were spun down at 300×g for 5 min at 4 °C. The supernatant was removed, and cells were mixed with 100 µL of lysis buffer (10 mM NaCl, 3 mM MgCl2, 10 mM Tris-HCl pH7.4, 0.1% Tween-20, 0.1% NonidetTM P40) and lysed on ice for 4 min. Wash buffer (1 mL; 10 mM NaCl, 3 mM MgCl2, 10 mM Tris-HCl pH7.4, 0.1% Tween20) was added before centrifuging at 500×g for 5 min at 4 °C. ATAC-seq libraries were generated as described 14  www.nature.com/scientificreports/ ATAC-seq pipeline (https:// github. com/ ENCODE-DCC/ atac-seq-pipel ine) using the MACS2 peak-caller with all default parameters. Differential accessibility was calculated between groups of subjects with different disease severity using the csaw package 35 .
Single-cell (sc)RNA-seq and analysis. Frozen PBMCs were thawed, and count and cell viability were measured by Countess II. The cell viability exceeded 80% for all samples except PBMC samples from CC subjects, which had viability between 70 and 80%. For single cell (sc) RNA-seq, 200 K cells were aliquoted, spun down, resuspended in 30 µL PBS + 0.04%BSA + 0.2U/µL RNase inhibitor, and counted using Countess II. GEM generation, post GEMRT cleanup, cDNA amplification, and library construction were performed following 10X Genomics Single Cell 5' v1 chemistry and quality was assessed using Agilent DNA Screentape assay. Libraries were then pooled and sequenced using Illumina NovaSeq 6000 platform with the goal of reaching saturation or 20,000 unique reads per cell on average. Sequencing data were used as input to the 10× Genomics Cell Ranger pipeline to demultiplex BCL files, generate FASTQs, and generate feature counts for each library. Dimensionality reduction and cell type annotation was accomplished using gene-barcode matrices generated using CellRanger count were analyzed using Seurat 3 with the default parameters unless otherwise specified 36 . Regulatory network inference was accomplished by converting the scRNA-seq Seurat object into a SingleCellExperiment and used as input to analysis with the SCENIC package 37 .
scATAC-seq and analysis. PBMCs were thawed and nuclei were extracted as for ATAC-seq. The singlecell suspensions of scATAC-seq samples were converted to barcoded scATAC-seq libraries using the Chromium Single Cell 5′ Library, Gel Bead and Multiplex Kit, and Chip Kit (10 × Genomics). The Chromium Single Cell 5′ v2 Reagent (10× Genomics, Cat# 120237) kit was used to prepare single-cell ATAC libraries according to the manufacturer's instructions. Quality was assessed using Agilent DNA Screentape assay. Libraries were then pooled and sequenced using Illumina NovaSeq 6000 platform with the goal of reaching saturation or 25,000 unique reads per nuclei on average. Sequencing data were used as input to the 10× Genomics Cell Ranger ATAC pipeline to demultiplex BCL files, generate FASTQs, and generate feature counts for each library.scRNA-seq and scATAC-seq were integrated using fragment file outputs generated using CellRanger ATAC count were analyzed using ArchR following the standard workflow and with default parameters unless otherwise specified 38 . Feature and motif enrichment analysis (peak calling) was performed using MACS2 via the addReproduciblePeakSet() method in ArchR which uses pseudo-bulk replicates of cells grouped on a specific design variable. The correlations between chromVAR transcription factor deviation scores and scRNA-seq derived gene expression data were calculated using the ArchR method correlateMatrices() to identify activators and repressors. Topic-based clustering was performed for the CD14+ monocytes from the mild or moderate subject cohorts using the R package cisTopic 39 . Detailed methods for qPCR and antibody tests, and epigenetic and genomic profiling, including RNA extraction, sequencing, differential expression analysis and differential chromatin accessibility analysis for both bulk and single cell analytic approaches can be found in Supplemental Information.

Data availability
The sequencing datasets and related clinical metadata tables are available via the NIH/NCBI Gene Expression Omnibus (GEO) repository using accession number GSE206284.