Dynamic activity in cis-regulatory elements of leukocytes identifies transcription factor activation and stratifies COVID-19 severity in ICU patients

Abstract

Transcription factor programs mediating the immune response to coronavirus disease 2019 (COVID-19) are not fully understood. Capturing active transcription initiation from cis-regulatory elements such as enhancers and promoters by capped small RNA sequencing (csRNA-seq), in contrast to capturing steady-state transcripts by conventional RNA-seq, allows unbiased identification of the underlying transcription factor activity and regulatory pathways. Here, we profile transcription initiation in critically ill COVID-19 patients, identifying transcription factor motifs that correlate with clinical lung injury and disease severity. Unbiased clustering reveals distinct subsets of cis-regulatory elements that delineate the cell type, pathway-specific, and combinatorial transcription factor activity. We find evidence of critical roles of regulatory networks, showing that STAT/BCL6 and E2F/MYB regulatory programs from myeloid cell populations are activated in patients with poor disease outcomes and associated with COVID-19 susceptibility genetic variants. More broadly, we demonstrate how capturing acute, disease-mediated changes in transcription initiation can provide insight into the underlying molecular mechanisms and stratify patient disease severity.

Keywords: COVID-19; active cistrome; acute respiratory distress syndrome; biomarkers; critical care; disease stratification; endotyping; enhancer RNA; transcription factor activity; transcriptional regulation.

Graphical abstract

Figure 1

Differential cistrome activation in peripheral leukocytes of hospitalized COVID-19 patients (A and B) Longitudinal study design for sampled plasma and peripheral leukocytes of patients across different stages of lung injury. Ninety-seven samples were included for active cistrome analysis using capped short RNA-seq (csRNA-seq). For fatal cases, the collection ended after patients transitioned to comfort care. One patient declined resuscitation or intubation (DNR/DNI). (C) csRNA-seq captures short 5′-capped RNA from active cis-regulatory elements (gene promoters, enhancers), collectively termed transcription start regions (TSRs). (D) csRNA-seq identifies transcriptional activity (red) in enhancers (eRNAs) at the STAT5B (left), LITAF (center), and MX1 (right) loci. The STAT5B enhancer resides in a neutrophil-specific active chromatin region marked by acetylated histone H3K27 (H3K27ac+) and open chromatin (ATAC-seq+). The LITAF enhancer resides in a CD4 T cell-active chromatin region. MX1 eRNA correlates with MX1 gene expression over time. (E and F) Unbiased clustering of 93,465 TSRs by activity similarity across 97 samples using uniform mani-fold approximation projection (UMAP). Positively associated TSRs (Pearson’s correlation; e.g., TSR A and TSR D) clustered closer together (F). The lines represent simple linear regression with 95% confidence intervals. (G) First-day enrollment leukocytes show globally distinct cis-regulatory activity between COVID-19 patients with fast (n = 5) and prolonged recovery (n = 9). Asymptomatic individuals (n = 5) serve as controls. See also Figures S1 and S2.

Figure 2

Lung injury-associated immune cistrome activation (A and B) The enhancer activities in the STAT5B intron (TSR A) and LITAF intron (TSR B) are correlated with the modified Murray score for lung injury (Spearmen correlation, n = 92 samples with clinical scores). The lines represent simple linear regression and the 95% confidence intervals. (C) A correlation coefficient was computed for each TSR (n = 93,465) based on its activity with lung injury scores. (D) Each TSR in the UMAP is shaded based on the lung injury correlation coefficient (purple, positively correlated; gold, negatively correlated). (E) Open-chromatin ATAC-seq enrichment in neutrophils, plasmablasts, monocytes, and CD 8 T cells for each TSR visualized on the immune cistrome UMAP. Part of a neutrophil’s active cistrome is associated with lower lung injury scores (asterisks, D and E). Red delineates high ATAC-seq enrichment relative to other hematopoietic cell types based on Calderon et al. and Perez et al. (F) Conceptual schematic using logistic regression analysis to identify transcription factor (TF) motifs associated with high (purple) or low (gold) lung injury scores (MEIRLOP). (G) Each dot represents the enrichment coefficient of TF motifs in TSRs with activity correlated or anticorrelated with the lung injury score. Error bars represent the 95% confidence intervals. The enrichments of all motifs, except for T1ISRE, are statistically significant (adj. p < 1 × 10⁻⁴). See also Figures S2–S4.

Figure 3

Distinct cistrome clusters identify co-enrichment of TF motifs suggestive of co-regulatory mechanisms (A) Cluster-specific motif analysis shows enrichment of the NF-κB binding sequence upstream of transcription start sites (TSSs) in clusters A and B but not C (green histograms). The IRF motif is enriched in cluster C (yellow). The incidence of NF-κB or the IRF motif in the entire active cistrome was used as a control (n = 93,465, gray). Clusters A and C have increased enrichment for CEBP (red) and PU.1 motifs (blue), lineage-determining TFs (LDTFs) for myeloid cells. (B) Discrete TSR clusters enriched for the representative TF motifs. (C) Motif analysis depicting co-enrichment of signal-dependent, LDTF, and promoter TF motifs. Red indicates motif enrichment in the TSR cluster relative to all TSRs (log2 ratio) and blue depletion. The dot size represents the Fisher’s exact test p value. Functional enrichment/gene ontology (GO) analysis identifies top pathways from genes associated with each TSR cluster. See also Figures S5–S7.

Figure 4

The kinetics of transcriptional factor activity during the hospitalization course of COVID-19 (A) cis-regulatory activity from three COVID-19 ARDS patients on enrollment days 1 and 9. (B) Full-time-course activities for T1ISRE/STAT and STAT/BCL6 are represented by the median log2 TSR signals and the 95% confidence intervals (CIs). Dotted boxes represent day 1 and 9. The medians for healthy controls (HCs) are shown as green dashes. (C) Heatmap of individual survivors’ time courses for T1ISRE/STAT and STAT/BCL6 activity. Red, high; blue, low; gray, no data. (D) Aggregate activity for T1ISRE/STAT and STAT/BCL6 TSRs (n = 9 patients, total 63 samples). Each point represents the median log2 of TSR clusters. The color of each point indicates the lung injury score (purple, severe lung injury; gold, less severe). The line and shaded region correspond to the smooth conditional mean and 95% CIs, respectively. (E) A correlation coefficient was computed for each TSR (n = 93,465) based on its activity with the hospital time of the prolonged survivors. (F) Each TSR in the UMAP is shaded based on the hospital time correlation coefficient (blue, early hospital course; yellow, late). (G) Gene pathways enriched in the early and late TF programs. A ridge plot shows the time-TSR activity correlation coefficient of genes in the respective pathways. See also Figures S8–S10.

Figure 5

Distinct TSR clusters exhibit significant enrichment of single-nucleotide polymorphisms (SNPs) associated with COVID-19 clinical outcomes (A) The LZTFL1 locus in chromosome 3p21.31 exemplifies the overlap of COVID-19-associated SNPs (p < 5 × 10⁻⁸) with cis-regulatory elements detected by csRNA-seq. (B) SNP A (rs34460587, −log p > 15 for hospitalization versus non-hospitalized COVID-19 cases) lies within −300 to +100 bp of the TSS located in the intergenic region between the CCR1 and XCR1 genes. (C) UMAP showing the global distribution of TSRs overlapping COVID-19-associated SNPs. (D) TSRs associated with higher lung injury or TF clusters (T1ISRE/STAT, STAT/BCL6, E2F/MYB, and ARE/SMAD/AP1) have significant overlaps with disease-associated SNPs. The color represents the adjusted p values, determined by regulatory element locus intersection (RELI), that account for the underlying genetic structure.

Figure 6

Dysregulated E2F/MYB, STAT/BCL6, and T1ISRE/STAT programs are associated with severe COVID-19 (A) csRNA-seq and total RNA-seq performed on matched RNA from peripheral leukocytes of COVID-19 patients (n = 55). To represent TF activity, genes specific to the TF network with expression highly correlated with the TSR cluster activity (csRNA-seq) were selected. (B) The median TF-specific gene expression positively correlated with the median cistromic activity for the E2F/MYB, STAT/BCL6, and T1ISRE/STAT clusters, respectively. Pairwise Pearson’s correlation analysis also shows specificity in TF activity. (C–F) Combinatorial TF activity signatures are associated with COVID-19 outcomes using a published validation cohort (n = 100). TF activities for STAT/BCL6, E2F/MYB, and T1ISRE/STAT were computed and ranked with bulk leukocyte RNA-seq (C). Hospital-free day 45 (HFD45) delineates clinical severity. TF activity signatures of (D) high STAT/BCL6 and high E2F/MYB or (E) high STAT/BCL6 and low T1ISRE/STAT correlate with the highest clinical severity. Combination of T1ISRE/STAT with E2F/MYB does not correlate with severity outcome (F). Error bars represent 95% CIs from the mean (horizontal lines). One-way ANOVA with adjusted p values; ∗p <0.05, ∗∗p < 0.005, ∗∗∗p < 0.0005. See also Figures S11–S14.

Figure 7

Summary (A) The active cistromes from patients’ leukocytes were profiled using csRNA-seq. Distinct clusters of cis-regulatory elements with co-enrichment of TF motifs, or “TF networks,” track clinical disease severity. (B) The activities of the TF networks in early or late hospital progression prognosticated outcomes.