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.


Supplemental
SARS-CoV-2 quantification by qRT-PCR Nasal swab Viral Transport Medium (VTM) was aliquoted and cryopreserved from study subjects to determine SARS-CoV-2 N1 gene copy number by RT-PCR to stratify subjects as COVID PCR positive or negative. Viral RNA was extracted from 140 uL of VTM according to manufacturer's instructions (QiaAmp Viral RNA minikit).
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). Standard quantitative RT-PCR (TaqPath 1-step RT qPCR Master Mix, Thermofisher) was run with test RNA and gene-specific standard curves (2e5 copy/mL -20 copy/mL). Regression analysis was used to determine gene copy number and corrected to report copies/mL of VTM. Samples with a Ct value less than 35 were scored as COVID PCR negative and samples greater than or equal to 35 were scored COVID PCR positive.

Purification of PBMCs
PBMCs were purified using the Ficoll-Hypaque density gradient method according to manufacturer's instructions. Briefly, whole blood was collected in ACD Vacutainer tubes (BD) and processed within 8 hours by dilution 1:2 in PBS, layered onto the Ficoll-Hypaque (Sigma Aldrich) in 50 ml conical tubes, and centrifuged at 420 x g for 25 minutes. Buffy coat was collected, washed twice in D-PBS by centrifugation at 400 x g for 10 minutes to isolate peripheral blood mononuclear cells (PBMCs) which were assessed for viability and cell count using a Vi-Cell automated cell counter (Beckman-Coulter). PBMCs were adjusted to 10x10 6  Differential chromatin accessibility analysis Differential accessibility was calculated between groups of subjects with different disease severity using the csaw v1 chemistry and quality was assessed using Agilent DNA tape screen 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 10x Genomics Cell Ranger pipeline to demultiplex BCL files, generate FASTQs, and generate feature counts for each library. For dimensionality reduction and cell type annotation, gene-barcode matrices generated using CellRanger count were analyzed using Seurat 3 with the default parameters unless otherwise specified 10 . Cells with > 5% of reads mapping to the mitochondrial genome or > 2500 genes detected were removed from the analysis. Counts were log-normalized, and the top 2000 variable features were identified. Principal component analysis was performed using these variable genes, and the top 20 principal components were used for downstream analysis. UMAP dimensionality reduction was performed using the top 20 principal components identified using the Harmony package 11 . Graph-based clustering was performed with resolution = 0.5. Cell types were inferred by using the DatabaseImmuneCellExpressionData() function from the SingleR package 12 . Labels were confirmed by identification of differentially expressed genes using the FindAllMarkers() function from Seurat 10 . For regulatory network inference, the scRNA-seq Seurat object was converted into a SingleCellExperiment and used as input to analysis with the SCENIC package 13 . Cells from seronegative mild and moderate subjects were re-clustered using Monocle 3, and the top 100 marker genes were computed for each cell partition 14 . The standard workflow for running the SCENIC analysis was then performed using the count matrix for these marker genes as input (https://github.com/aertslab/SCENIC). Briefly, GENIE3 was used to identify regulons of transcription factors and their downstream regulatory targets with correlated coexpression, and AUCell was then used to score the activity of these regulons in each cluster. The 'top10perTarget' co-expression parameter value was used to prune the list of scored regulons. Gene set and pathway enrichment was measured by using the Seurat method AddModuleScore() for lists of genes retrieved from the Molecular Signatures Database v7.5.1 (https://www.gsea-msigdb.org/gsea/msigdb/) as documented in Supplemental Table   11. scRNA-seq clusters were merged as follows: Mild/Mod. 1 = clusters 9, 10; Mild/Mod. 2 = cluster 8; Mild/Mod. 3 = cluster 11; Mild/Mod. 4 = clusters 1, 3, 5, 6. The average pathway score was plotted on a heatmap for each cell subset, and hierarchical clustering was applied to group similar samples within each dataset. Oneway analysis of variance (AOV) and Tukey testing was performed comparing the average pathway score for each cell subset. The AOV adjusted p-value and the percentage of cell subset contrasts that were significantly associated with the average pathway score were recorded in Supplemental Table 12.
Single-cell (sc)ATAC-seq, data processing and analysis PBMCs were thawed and nuclei were extracted as for ATAC-seq. For scRNA-seq and scATAC-seq integration, fragment file outputs generated using CellRanger ATAC count were analyzed using ArchR following the standard workflow and with default parameters unless otherwise specified 15 . Cells with a transcription start site enrichment score < 4, cells with fewer than 1000 detected fragments, and putative doublets were removed from downstream analysis. Dimensionality reduction was computed using iterative latent semantic indexing (LSI), and batch effect correction was applied using Harmony.
Graph-based clustering was performed using the FindClusters() method from Seurat 3 with resolution = 0.8.
UMAP embeddings were calculated with the top 30 principal components from either LSI or Harmony.
Constrained integration was performed using the addGeneIntegrationMatrix() method and scRNA-seq cell type annotations were used to label the identify of scATAC-seq clusters. Pseudo-scRNA-seq gene expression profiles were estimated for each cell type using the addGeneIntegrationMatrix() function from ArchR during the constrained label transfer procedure.
For 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. Differentially accessible peaks were identified between two groups and visualized using the ArchR methods getMarkerFeatures() and markerPlot(), respectively. Significance was defined as FDR <= 0.1 and absolute log2 fold change >= 0.5 unless otherwise specified. The 'cisbp' motif set was imported from TFBSTools using the ArchR addMotifAnnotations() method, and motif enrichment in differentially accessible peaks was performed using the peakAnnoEnrichment() method. Accessibility estimates were compared to the average gene expression over all cells for the genes that encode each transcription factor. The gene integration matrix was constructed using the .getGroupMatrix() function from ArchR and specifying which subject cohorts to compare.
For integrative analysis with scRNA-seq, the correlations between chromVAR transcription factor deviation scores and gene expression data were calculated using the ArchR method correlateMatrices() to identify activators and repressors. Peak-to-gene linkages were calculated using the addPeak2GeneLinks() method in ArchR using a correlation cutoff of 0.5 and resolution = 1. This approach uses low-overlapping cell aggregates to reduce noise that arises from doing correlative analyses with sparse scATAC-seq datasets. Peak-to-gene linkages were plotted against peak accessibility at gene bodies within DORCs for each cell type. DORC genes were defined as gene loci with > 10 peak-to-gene linkages, and these sites were used as input to the web tool Seanalysis to identify regulation by a known super-enhancer in peripheral blood cells 16,17 . The activity of the top-ranked transcription factor regulators that were correlated with scRNA-seq clusters was estimated for each scATAC-seq cluster by measuring the average genome-wide motif accessibility. These activities were used to calculate Pearson correlation coefficients between scATAC-seq clusters (C1-C8) and scRNA-seq clusters (C1-C11) were calculated to identify scATAC-seq clusters with similar regulatory network activity.