Dissecting PGE_2-driven inhibition of T cell activation using single-cell multi-omic and inflammatory bowel disease genetic association analysis

Domain: chemistry
Published: 2026-04-07
Editorial outcome: minor_revision

Abstract
Genetic and immunologic studies implicate the interleukin (IL)-23/T helper (Th)17 pathway in
inflammatory bowel disease (IBD). IL-23 and IL-1β drive human Th17 differentiation, while
prostaglandin E2 (PGE ) and transforming growth factor-β (TGF-β) further modulate Th17
2
development and plasticity. However, how these inflammatory mediators influence human T cell
regulatory programs remains incompletely understood. We used single-cell multi-omics to profile
171,829 peripheral blood T cells from 25 healthy donors in 64 samples exposed to activation
stimuli, IL-1β and IL-23 alone or in combination with PGE , TGF-β, or both. PGE broadly
2 2
suppressed T cell activation, except in Th17 and T follicular helper cells, and markedly altered
chromatin accessibility and gene expression, particularly in Th17, Th1, and regulatory T cells,
where IBD-associated SNPs were enriched in open chromatin and connected to cell type-specific
cis-regulatory elements. Our study demonstrates the utility of single-cell multi-omics for defining
stimulus-specific effects on T cells and prioritizing disease-associated genes.

Introduction
The inflammatory bowel diseases (IBD), Crohn's disease (CD) and ulcerative colitis (UC), are
chronic, immune-mediated inflammatory disorders of the gastrointestinal tract that affect about
3.7 million people in America and Europe, and hundreds of thousands more globally1-3. The
pathogenesis of IBD is multifactorial, involving both genetic and non-genetic elements such as
environmental factors and gut microbiota2. Further research to unravel the cellular and molecular
mechanisms of IBD is needed for development of more targeted, effective therapies.

Several converging lines of genetic and immunologic evidence suggest that interleukin (IL)-23
signaling and the T helper (Th)17 immune pathway play key roles in the pathogenesis of IBD4-15.
IL-23 and IL-1β drive human Th17 differentiation, and prostaglandin E2 (PGE2) and transforming
growth factor-beta (TGF-β) have also been reported to influence Th17 development and/or
plasticity11-13,16-39. These four inflammatory mediators are also overexpressed in IBD tissues40-45,
and genes that encode their receptors (IL23R, IL1R1, and PTGER4) and/or downstream signal
transducers (e.g., TYK2, JAK2, STAT3, SMAD3) are implicated in genome-wide association studies
(GWAS) that have identified thousands of IBD-associated variants4-7.

The strongest association signals in most IBD risk loci lie in non-coding genomic regions where
causal variants cannot be distinguished from neighboring non-causal variants in linkage
disequilibrium, and their presumed cell-type and microenvironment-specific regulatory effects on
gene expression is usually not obvious. Integration of multimodal omics data in disease-relevant
cell types and contexts with epigenetic fine mapping promises to bridge knowledge gaps between
statistical genetic associations and understanding of underlying disease mechanisms46-53.

These rationales motivated us to isolate T cells from the peripheral blood of 25 healthy human
donors with ages 18-35 years, and to expose the T cells to activation stimuli, IL-1β and IL-23 alone
or in combination with PGE2, TGF-β, or both PGE2 and TGF-β to model defined inflammatory
signaling contexts relevant to IBD. By studying T cells from healthy donors with an age range
similar to the peak ages of IBD onset in adults, confounding variables such as IBD-directed
medications, comorbidities, and age-related immune changes that could influence T cell
responses were avoided. The T cells were then profiled using DOGMA-seq, a single-cell trimodal
omics technique incorporating simultaneous assays for gene expression, cell surface protein
epitope antibody-derived tags (ADTs), and the assay for transposase accessible chromatin (ATAC)
within individual cells, enabling higher resolution annotation of specific T cell subtypes and more
robust dissection of complex regulatory networks in comparison to single modal or bimodal
assays54-56.

Leveraging this multimodal omics dataset, we systematically characterized T cell heterogeneity,
quantified PGE2-associated transcriptional and epigenetic changes, and integrated these signals
with IBD GWAS single nucleotide polymorphisms (SNPs) to prioritize candidate genes and putative
cis-regulatory elements. Using IBD genetics as a test case, we demonstrated how this resource
can be applied to link disease-associated non-coding variants to cell type-specific regulatory
programs. To validate key findings, we additionally analyzed bulk ATAC-seq and RNA-seq data
from Th17- and Treg-enriched CD4+CD45RO+CD196+ T cells.

In this work, we establish a single-cell multi-omics T cell resource that advances understanding of
how inflammatory cues intersect with genetic risk for IBD. Beyond IBD, this framework provides
a broadly applicable platform for interpreting genetic associations and guiding functional studies
across immune-mediated diseases, including other disorders involving the IL-23/IL-17 pathway,
such as psoriasis and ankylosing spondylitis8,15.

Results
Simultaneous trimodal single-cell measurement and high-resolution annotation of human
peripheral blood T cells
To characterize the various immune cell states in T cells and the gene regulatory elements in their
subpopulations, we profiled 171,829 high-quality T cells from 25 healthy human donors using
DOGMA-seq (Fig. 1A, Supplementary Fig. 1A). T cells from the first seven study participants were
exposed to anti-CD3/CD28/CD2 activation stimuli, IL-1β and IL-23 (Act IL1B IL23) alone, and in
combination with either PGE2 (Act IL1B IL23 PGE2) or TGF-β (Act IL1B IL23 TGFB), or both PGE2
and TGF-β (Act IL1B IL23 PGE2 TGFB) (7 study participants x 4 conditions), but after we observed
minimal effects of TGF-β (details presented in the Cell type-specific changes upon exposure of T
cells to PGE2 section), cells from the remaining 18 study participants were exposed to just the
conditions without and with PGE2 (Act IL1B IL23 and Act IL1B IL23 PGE2) (18 study participants x
2 conditions), with distinct hashtags added to each of the 64 samples to mitigate batch effects
between samples collected from the same study participant but stimulated differently
(Supplementary Fig. 2A, Supplementary Fig. 1B-C).

We simultaneously assayed 36,601 genes, 163 ADTs, and 252,371 ATAC peaks for each profiled
cell, and this facilitated precise annotation of T cell subpopulations based on canonical markers
in all three modalities (Supplementary Fig. 2B-D). ADTs were particularly effective in
differentiating between CD4+ and CD8+ T cells, and between naïve and memory T cells.
Additionally, ATAC peaks and their associated transcription factor binding motifs57 played a vital
role in identifying T cell subtypes, such as Th1 and Th17 cells.

To address the dynamic changes in T cells upon activation58, we computed activated/resting
scores for each cell based on two lists of signature genes (Fig. 1B, Supplementary Fig. 2E). The
scores are highly correlated with CD25 and CD127 ADTs (Fig. 1C). As both up-regulation of CD25
and down-regulation of CD127 are associated with T cell activation59,60, this finding provides
independent validation of the activated/resting scores. These scores serve as a quantitative proxy
for the continuum of T cell activation states, enabling systematic characterization of activation
heterogeneity across subpopulations61. Leveraging both canonical markers and activated/resting
scores, we identified 20 subtypes of T cells (Fig. 1D, Supplementary Fig. 3), including CD4+ naïve,
CD4+ regulatory, CD4+ memory (further split into four subtypes), CD8+ naïve, CD8+ regulatory,
CD8+ memory, and gamma-delta T cells, and mucosal-associated invariant T cells (MAITs).
Importantly, each cell type is comprised of cells from all study participants, and CD4+ T cells
account for more than half of the T cell population (Supplementary Fig. 1D). These
subpopulations are particularly distinguishable based on ADT and ATAC data (Supplementary Fig.
1E), highlighting the value of multi-omics in comprehensively characterizing T cells.

To further validate the activated/resting scores, an integrative analysis was performed to
combine in-house and public DOGMA-seq data, leveraging the ground truth of T cells treated
with and without anti-CD3/CD28 beads54 (Supplementary Fig. 4A-B). We observed that the
distribution of untreated T cells matches cells with high resting scores from our data, while the
majority of treated T cells overlap with cells exhibiting high activated scores but still contain cells
in the resting state (Supplementary Fig. 4C-D), suggesting that some T cells treated with anti-
CD3/CD28 beads were not fully activated, at least at the time of profiling (16h), or may represent
the lowly active state defined in a recent study58.

Cell type-specific changes upon exposure of T cells to PGE2
When we analyzed 84,345 T cells from seven study participants that were stimulated under four
distinct conditions (Supplementary Fig. 2A, 5B), we observed that the variation between
conditions with and without PGE2 matches the resting and activated states, but we observed
little difference between conditions with and without TGF-β (Supplementary Fig. 5C-D). This
observation is consistent with analysis of data from activated and resting T cells (Supplementary
Fig. 4E). Moreover, we observed significant changes in cell type proportions between conditions
with and without PGE2 but not between conditions with and without TGF-β (Supplementary Fig.
5A). Given these minimal effects of TGF-β, we focused subsequent experiments on two
stimulation conditions, with and without PGE2 (Supplementary Fig. 2A).

Two-sided paired Wilcoxon tests were performed to compare changes in cell type proportions in
129,845 cells from the 25 donors that were stimulated under conditions with and without PGE2.
PGE2 significantly suppresses the activation of multiple T cell subtypes, including CD4+ naïve,
CD4+ regulatory, other CD4+ memory, CD8+ naïve, CD8+ memory, and MAITs (Fig. 2A), which is
consistent with previous reports25. Interestingly, we did not observe immunosuppressive effects
in Th17 and T follicular helper (Tfh) cells.

Differential chromatin accessibility (DA) and differential gene expression (DE) analyses were
performed in 50 paired pseudo-bulk samples (25 study participants x 2 conditions) in each cell
type (Supplementary Data 1-2). While PGE2 suppresses T cell activation, it also up-regulates
numerous ATAC peaks and genes in cells of the resting state, particularly in Th1, Th17, Tfh, and
MAIT cells (Fig. 2B-C). Notably, CREM, STAT4, and PDE4D are among the top up-regulated genes.
Conversely, PTGER4, which encodes the EP4 receptor for PGE2, is strongly down-regulated
(Supplementary Fig. 6-7).

Pathway enrichment analysis reveals that PGE2 exerts opposing effects based on T cell activation
state. Key signaling pathways such as Protein Kinase A (PKA) (including CREM and PDE4D), mTOR,
and JAK/STAT (including STAT4) are up-regulated in resting T cells but down-regulated in
activated T cells (e.g., Th1, other CD4+ memory, CD8+ memory, and MAIT cells), demonstrating
a dual regulatory role (Fig. 2D).

These findings support a model where PGE2 initiates a transient immunomodulatory program via
cAMP–PKA signaling, as evidenced by the up-regulation of CREM and PDE4D62;63. While the
inducible cAMP early repressor isoforms of CREM function as transcriptional repressors of
inflammatory genes, the concurrent up-regulation of PDE4D (which degrades cAMP) and down-
regulation of PTGER4 suggest a negative feedback loop that may prevent prolonged
immunosuppression by self-limiting the signal.

PGE2-induced cAMP has long been recognized for its anti-inflammatory actions64. Interestingly,
while canonical cAMP signaling typically inhibits nuclear factor-kappaB by suppressing IκB kinase
activity65, we observed a noncanonical response in resting Th17 cells, where exposure to PGE2
was associated with transcriptional up-regulation of CREB1, NFKB1, and other cAMP pathway
components, accompanied by increased chromatin accessibility at the NFKB1 locus, suggesting
active transcriptional remodeling. This suggests that cAMP effects may be highly context-
dependent66, and is consistent with murine studies showing that PGE2–EP2/EP4 signaling
cooperates with IL-23 to enhance pathogenic Th17 differentiation67.

In addition, in cells of the resting state, increased activity in the mTOR pathway may reflect
enhanced metabolic preparedness68 and EIF2 pathway up-regulation may indicate translational
reprogramming69.

In parallel, CD4+CD45RO+CD196+ T cells were sorted to enrich Th17 and Treg cells, which were
subsequently activated, stimulated, and subjected to bulk ATAC-seq (cells collected from 33
healthy study participants) and RNA-seq (cells collected from 38 healthy study participants,
Supplementary Fig. 2A). Using bulk samples, CREM, STAT4, PDE4D, and PTGER4 are highlighted
again in DA and DE analyses (Supplementary Fig. 8A, C, Supplementary Data 3-4). To validate
single-cell DA and DE results, we conducted correlative analysis to compare effects observed in
bulk analysis to those observed in comparable single-cell analysis (Supplementary Fig. 8B, D), and
we found high directional concordance. Despite reduced granularity compared to single-cell
analysis (Fig. 2A), bulk analysis still reveals the immunosuppressive effects of PGE2
(Supplementary Fig. 8E).

Allele-specific expression (ASE) and open chromatin (ASoC) analysis reveals stimulation
dependent regulatory dynamics and enrichment in IBD loci
To investigate the impact of stimulation on cis-regulatory variation, we performed ASE and ASoC
analyses in samples from the 25 donors with DOGMA-seq data. ASE refers to the imbalance in
gene expression between two alleles at heterozygous loci, while ASoC reflects allele-specific
differences in chromatin accessibility. High-confidence genotype data for each donor was
obtained directly from the multi-modal single-cell sequencing data by performing germline
variant calling separately from the gene expression and chromatin accessibility modalities of
DOGMA-seq. After discordant variants between RNA- and ATAC-derived SNP calls, flipped alleles
and duplicates, and variants with high missingness were filtered, a total of 60,962 high-
confidence SNPs remained (see Methods for details).

We identified 8,315 significant ASE sites and 18,524 significant ASoC sites (FDR-adjusted P < 0.05)
(Supplementary Data 5-6) Their distribution across each T cell subtype is shown in Supplementary
Fig. 9A. These allele-specific signals are broadly distributed but vary with stimulation condition
and activation state. In the Act IL1B IL23 stimulation condition, activated T cells show a higher
number of significant ASE and ASoC sites than their resting counterparts (paired Wilcoxon rank-
sum P: 2.82x10-3 [ASE] and 3.91x10-3 [ASoC]). However, in the Act IL1B IL23 PGE2 stimulation
condition, this pattern is reversed: resting T cells exhibit significantly more ASoC sites (paired
Wilcoxon rank-sum P: 1.37x10-2), while ASE site numbers are not significantly different (paired
Wilcoxon rank-sum P: 0.213).

To further dissect the impact of PGE2 stimulation, we compared the numbers of significant ASE
and ASoC sites between conditions with and without PGE2 within each T cell subtype, under
otherwise identical stimulation with IL-1β and IL-23. We observe a striking state-dependent
effect. In activated T cells, PGE2 stimulation is associated with a significant reduction in both ASE
and ASoC signals (paired Wilcoxon rank-sum P = 9.77x10⁻³ for both), suggesting dampened allelic
imbalance and potentially reduced regulatory activity under inflammatory conditions. In contrast,
in resting T cells, PGE2 stimulation is associated with a significant increase in the number of both
ASE and ASoC sites (paired Wilcoxon rank-sum P = 1.95x10⁻³), indicating enhanced allele-specific
regulation in the absence of prior activation. These findings further support a model in which
PGE2 exerts inhibitory effects on activated T cells, while promoting greater regulatory
heterogeneity or plasticity in resting T cells. This stimulation-state interaction highlights context-
specific regulatory dynamics of T cells in response to external signals.

While many of the ASE signals are stimulation-specific, a substantial proportion are shared across
all four stimulation conditions, reflecting a stable core of allele-specific regulatory activity.
Specifically, 3,339 ASE variants and 6,198 ASoC variants are consistently detected across all
conditions (Supplementary Fig. 9B–C). In contrast, several hundred sites are unique to individual
conditions, indicating stimulus-specific regulatory effects. Notably, 2,100 variants exhibit
consistent allele-specific signals at both the chromatin and gene expression levels, suggesting
coordinated effects on regulatory architecture.

To evaluate the potential disease relevance of these regulatory variants, we assessed whether
stimulation-responsive ASE and ASoC sites are enriched in genomic regions identified in IBD
GWAS6. Both ASE and ASoC sites show strong enrichment in IBD GWAS loci, highlighting their
potential regulatory involvement in disease susceptibility. This enrichment is most pronounced
for shared sites across all stimulation conditions, but is also evident in stimulation-specific sites.
Specifically, ASE sites show pronounced enrichment with Wilcoxon rank-sum test observed vs.
expected P values of 2.23x10-6 for Act IL1B IL23 PGE2-specific sites, 1.54x10-3 for Act IL1B IL23
TGFB-specific sites, 1.63x10-2 for Act IL1B IL23 PGE2 TGFB-specific sites, 1.78x10-6 for Act IL1B
IL23-specific sites, and 1.53x10-113 for shared sites (Supplementary Fig. 9D). Similarly, ASoC sites
display significant enrichment with P values of 1.62x10-4 for Act IL1B IL23 PGE2-specific sites, 2.90
x10-4 for Act IL1B IL23 TGFB-specific sites, 2.34x10-2 for Act IL1B IL23 PGE2 TGFB-specific sites,
5.62x10-12 for Act IL1B IL23-specific sites, and 2.00x10-199 for shared sites (Supplementary Fig. 9E).
These findings indicate a robust association between ASE and ASoC sites and IBD genetic risk loci,
underscoring the potential functional relevance of these regulatory elements in the pathogenesis
of IBD.

Cell type-specific enrichment of IBD GWAS SNPs in ATAC peaks
Our comprehensively characterized T cell atlas with stratification into activated and resting states
provides an unprecedented opportunity to investigate the T cell subtype specificity of IBD GWAS
SNPs across T cell subpopulations. We focused on 4,329 SNPs from 151 IBD credible sets, which
have a greater than 95% probability of containing the causal variant for each fine-mapped IBD
locus6,70. These variants, hereafter referred to as IBD credible set SNPs, were prioritized based on
Bayesian fine-mapping approaches that account for linkage disequilibrium and effect sizes. The
term 'all variants' refers to the 4,329 IBD SNPs contained within all of the IBD credible sets, while
'lead variants' refers to the 151 IBD SNPs with the most significant GWAS P-value within each IBD
credible set.

To assess cell type-specific enrichment, we used Chromatin Element Enrichment Ranking by
Specificity (CHEERS)71, a statistical method that detects enrichment of disease-associated SNPs
in regulatory regions by accounting for the quantitative, cell-state-specific activity of chromatin
elements. By leveraging the multi-sample design of our study and addressing the sparsity of
single-cell data, we generated 500 pseudo-bulk samples (20 cell types x 25 study participants)
and conducted SNP enrichment analysis using CHEERS, which yielded enrichment P-values for
each pseudo-bulk sample (Supplementary Fig. 10A ). Individual P-values for each cell type were
subsequently combined and adjusted across distinct cell types using the most conservative
methods to narrow down the cell types most worthy of post-GWAS functional investigation (Fig.
3A).

Both lead variants and all variants in the IBD credible set SNPs are strikingly enriched in ATAC
peaks of activated Th17 (Bonferroni-P = 1.36x10-34 [lead] and 3.05x10-30 [all]), Th1 (Bonferroni-P
= 6.64x10-27 [lead] and 3.97x10-19 [all]), and CD4+ Treg cells (Bonferroni-P = 7.75x10-22 [lead] and
1.32x10-7 [all]), reflecting the known dysregulation of the balance between these cell types in IBD
pathogenesis72. In addition, our results also implicate involvement of CD8+ T cells, given the
notable enrichment of lead variants in ATAC peaks of MAITs (Bonferroni-P = 2.78x10-16 [activated]
and 7.76x10-14 [resting]), CD8+ memory (Bonferroni-P = 2.75x10-11 [activated]), and gamma-delta
T cells (Bonferroni-P = 1.10x10-10), which have gained attention in recent studies73,74. Notably,
CD-associated variants show more significant enrichment in activated states, especially in Th17
cells (Bonferroni-P = 1.12x10⁻²⁴ [CD] vs. 9.70x10⁻³ [UC]), suggesting that activated Th17 cells may
play an especially prominent role in CD pathogenesis.

To independently validate these findings, we applied Stratified Linkage Disequilibrium Score
Regression (S-LDSC), a heritability-based method that tests whether GWAS signal is enriched in
specific annotations across these cell types75,76, using IBD, CD, and UC summary statistics6. The
majority of our findings are consistent with those obtained through CHEERS analysis, with
consistent heritability enrichment in activated Th1, Th17, Tfh, CD4+ memory, and Treg cells in
the context of IBD (Fig. 3A).

Notably, a pattern of pronounced heritability enrichment in the activated states of these cells
compared to their resting counterparts emerged from our analysis. This trend is particularly
evident in activated Th1 (Bonferroni-P = 2.74x10-9 [IBD], 2.36x10-6 [CD], 9.20x10-7 [UC]), Th17
(Bonferroni-P = 1.66x10-10 [IBD], 1.37x10-7 [CD], 1.02x10-6 [UC]), Tfh (Bonferroni-P = 4.09x10-5
[IBD], 2.65x10-3 [CD], 1.74x10-3 [UC]), other CD4+ memory (Bonferroni-P = 1.53x10-9 [IBD],
2.03x10-8 [CD], 3.31x10-6 [UC]), and CD4+ Treg cells (Bonferroni-P = 1.41x10-5 [IBD], 2.70x10-6 [CD],
8.15x10-4 [UC]), where significant enrichment is observed across the IBD, CD, and UC phenotypes.

Prostaglandins are involved in diverse biological processes including IBD-related inflammation
and immune responses77, so we next examined whether PGE2 alters enrichment patterns of IBD
risk variants. We generated 480 pseudo-bulk samples (20 cell types x 2 conditions x 12 randomly
selected study participants, Methods) and conducted CHEERS analysis to compare the
enrichment of IBD-associated SNPs in samples stimulated with and without PGE2 (Supplementary
Fig. 10B), followed by combining and adjusting individual P-values (Fig. 3B). The enrichment of
IBD-associated SNPs in cell types is similar to the previous analysis without stratification by
stimulation condition. Specifically, lead variants in the IBD credible set SNPs are significantly
enriched in ATAC peaks of activated Th17 (Bonferroni-P = 5.87x10-27 [w/o PGE2] and 1.34x10-18
[w/ PGE2]), Th1 (Bonferroni-P = 2.21x10-13 [w/o PGE2] and 1.46x10-18 [w/ PGE2]), and CD4+ Treg
cells (Bonferroni-P = 2.56x10-18 [w/o PGE2] and 2.07x10-13 [w/ PGE2]), both with and without
PGE2 stimulation. Interestingly, IBD-associated SNPs show more significant enrichment in resting
Th17 (Bonferroni-P = 1.83x10-15 [w/o PGE2] and 1.05x10-7 [w/ PGE2], difference-P = 6.13x10-3),
other CD4+ memory (Bonferroni-P = 3.62x10-5 [w/o PGE2] and 1 [w/ PGE2], difference-P =
1.89x10-4), and CD8+ memory T cells (Bonferroni-P = 8.64x10-6 [w/o PGE2] and 1 [w/ PGE2],
difference-P = 8.64x10-4) without PGE2 stimulation, despite their lower cell type proportions in
the absence of PGE2 stimulation. In contrast, IBD-associated SNPs show more significant
enrichment in activated CD8+ memory T cells (Bonferroni-P = 9.67x10-9 [w/o PGE2] and 4.84x10-
17 [w/ PGE2], difference-P = 1.53x10-2) upon PGE2 stimulation.

Differences in enrichment between conditions with and without PGE2 are relatively modest for
CD-associated SNPs, with significance detected only in resting MAITs (difference-P = 7.69x10-3
[lead] and 1.58x10-2 [all]), compared with those for UC-associated SNPs, suggesting that PGE2
differentially modulates chromatin accessibility at CD- and UC-associated risk loci. Notably, this
aligns with clinical observations that non-steroidal anti-inflammatory drug (NSAID) use, which is
known to reduce endogenous PGE2 levels, is associated with increased relapse risk in CD but not
UC78. Our findings suggest a possible mechanism: the reduced regulatory flexibility at CD risk loci
could increase vulnerability to PGE2 depletion during NSAID exposure, whereas the
comparatively greater plasticity of UC loci may help buffer against such effects.

Cell type-specific peak-to-gene linkages facilitating annotation of ATAC peaks and
identification of plausible cis-regulatory elements (CREs)
To better understand the regulatory function of non-coding chromatin regions, we identified
peak-to-gene linkages, which are correlations between chromatin accessibility peaks and
expression levels of nearby genes, thereby mapping candidate CREs to their putative gene
targets79. We initially used ArchR, a computational framework that enables scalable inference of
peak-to-gene regulatory relationships from single-cell data to identify peak-to-gene linkages
within a distance of 1 Mb for T cell subpopulations80. To improve robustness, we refined the
initial linkages from ArchR by validating them using pseudo-bulk correlation analyses, generating
64 pseudo-bulk samples for each cell type.

This analysis yielded an average of 101,972 plausible peak-to-gene linkages for each
subpopulation (Fig. 3C), constituting approximately 6.5% of peak-to-gene linkages initially
discovered using ArchR (Methods). On average, 76.1% (ranging from 60.6% to 90.9%) of plausible
peak-to-gene linkages were confirmed as positive associations. Notably, these results reveal that
regulatory interactions are highly cell type-specific and often extend beyond nearest-gene
assumptions. Only 8.5% (ranging from 4.4% to 14.0%) of the top linked genes were found to be
the nearest genes to their corresponding ATAC peaks (Fig. 3D), emphasizing the need to go
beyond proximity-based annotations when assigning regulatory elements to target genes.
Interestingly, we found that both CD4+ and CD8+ naïve T cells have fewer linked peaks and genes
(Fig. 3E-F), suggesting a more transcriptionally quiescent or less differentiated regulatory state
despite their high abundance. Linked peaks and genes shared between different cell types are
summarized in Supplementary Fig. 11A-B.

Together, these regulatory maps offer valuable insights into the functional annotation of non-
coding chromatin regions. By linking T cell subtype-specific chromatin accessibility peaks that
overlap disease-associated SNPs with downstream gene targets, they help clarify cellular
contexts in which disease-associated SNPs may exert pathogenic effects.

SCARlink reveals putative enhancers and supports regulatory interpretations of peak-to-gene
linkages
We applied SCARlink, a computational method that predicts enhancer-gene relationships using
joint single-cell gene expression and chromatin accessibility profiles, to augment understanding
of cis-regulatory mechanisms and annotation of enhancer elements relevant to IBD in T cells81.
SCARlink fits a regularized Poisson regression to predict single-cell gene expression based on
accessibility of fixed-size genomic tiles and calculates cell type-specific Shapley values, which
quantify the relative importance of each tile in explaining the expression of a given gene.
Genomic regions with high Shapley values are prioritized as candidate enhancers in a cell-type-
specific manner.

SCARlink was applied to 5,000 highly variable genes (HVGs) identified in DOGMA-seq data. Genes
with sparse expression (defined as having more than 90% zero values) were excluded, resulting
in a final set of 1,256 genes for analysis. Across all T-cell subpopulations, we identified an average
of 38,819 significant gene-linked tiles per subpopulation (ranging from 9,235 to 100,313) linked
to 1,247 genes on average (ranging from 1,209 to 1,254) after FDR adjustment. Importantly,
SCARlink also provided independent support for 4,916 peak-to-gene linkages per subtype (range:
1,130–11,075), based on the colocalization of predicted enhancer tiles within or overlapping
previously linked ATAC peaks. These results suggest that a substantial proportion of peak-to-gene
relationships may be mechanistically driven by cell type-specific enhancer activity.

The full set of SCARlink predictions provided in Supplementary Data 7 serve as a high-resolution,
cell-type-resolved map of putative enhancers in human T cells and provide mechanistic insights
into how noncoding variation at ATAC peaks may influence gene expression.

Cell type-specific colocalization of IBD-associated SNPs in plausible CREs
Next, we conducted a colocalization analysis to identify ATAC peaks that are both associated with
cis-gene expression (plausible CREs) and overlap all or lead IBD credible set SNPs, thereby
revealing the enrichment of IBD credible set SNPs in specific loci and prioritizing cell types crucial
to disease pathogenesis.

Focusing on the cell types prioritized by CHEERS and S-LDSC, all variants were found to be
strikingly enriched in peaks linked to PTGER4 (Fig. 4A). For instance, the ATAC peak chr5-
40490646-40491146, covering the SNP rs10440635 (Fig. 4B), displays significant positive
association with PTGER4 across seven distinct cell types, including activated and resting CD4+
naïve T, resting CD4+ Treg, activated Th1, resting Th17 and Tfh, and activated CD8+ naïve T cells
(Supplementary Fig. 12A). SCARlink analysis reveals that the predicted PTGER4 gene expression
shows a Spearman correlation R² of 0.244 (Supplementary Fig. 13), well above the genome-wide
average of 0.116, indicating that chromatin accessibility is a major determinant of PTGER4
expression in T cells and may mediate genetic risk for IBD. We identified a significant PTGER4-
linked tile (chr5-40490500-40491000) through Shapley value analysis and expression prediction
testing that overlaps both rs10440635 and the peak chr5-40490646-40491146 in five T cell
subtypes including resting and activated CD4⁺ naïve, activated Th1, resting Th17, and activated
CD8⁺ naïve T cells, supporting the presence of a cell type-specific enhancer that may regulate
PTGER4 and contribute to IBD pathogenesis (Supplementary Fig. 13 and Supplementary Data 8).
Of note, rs10440635 not only resides within the IBD credible set spanning from 40486896 to
40623346 on chr570 but it has also been implicated in ankylosing spondylitis, a condition closely
linked with IBD, and is presumed to regulate the expression of PTGER4 in synovial biopsies82.
Although not reported among 216,597 immune cell type-specific cis-expression quantitative trait
loci (eQTL) identified in single-cell eQTL analysis of peripheral blood mononuclear cells83,
rs10440635 was a PTGER4 eQTL in eQTLGen data84 from blood (P = 2.02x10-8) and achieved
nominal significance in GTEx data85 from EBV-transformed lymphocytes (Supplementary Fig. 12B),
thereby highlighting the need for eQTL analysis with finer granularity to verify this plausible
association.

Similarly, another 12 IBD-associated SNPs are covered by an additional 10 ATAC peaks that are
correlated with PTGER4 across 13 subtypes of T cells (Supplementary Fig. 12A), including
rs9292781 and rs1122433 that map within chr5-40578427-40578927 (Fig. 4B), and were both
PTGER4 eQTLs in eQTLGen data84 from blood (P = 4.77x10-25 and 4.42x10-30). The tile chr5-
40578500-40579000, which includes rs1122433, is significantly linked to PTGER4 and overlaps
the peak-to-gene linkage between PTGER4 and chr5-40578427-40578927 in activated and
resting CD4+ naïve, activated CD8+ naïve, activated Th1, and other activated CD4+ memory T
cells, suggesting an additional enhancer region that may contribute to PTGER4 regulation and
IBD risk (Supplementary Data 8). These findings converge on the PTGER4 locus as a genetically
and epigenetically regulated target in T cells, strongly implicating PGE2–EP4 signaling in IBD
pathogenesis.

Other colocalized SNP–CRE relationships were also identified. Around LNPEP, a gene that has
received little attention in the context of IBD despite evidence indicating a decrease in its serum
protein abundance in CD patients86, 35 IBD-associated SNPs fall within 14 peaks correlated with
LNPEP expression across 16 cell types (Supplementary Fig. 14A, 15). Interestingly, the ATAC peak
chr5-96896251-96896751 covers eight IBD-associated SNPs and is a plausible CRE for LNPEP in
eight cell types (Supplementary Fig. 14B). Analysis of all variants also reveals enrichment in linked
peaks for TNFRSF14 (e.g., rs2227313 in resting CD8+ naïve T cells and six other cell types, which
was also revealed as a single-cell eQTL in multiple T cell subpopulations83, Supplementary Fig. 16-
17), and IL2RA (e.g., rs61839660, which has a posterior probability of 0.999 to be a causal variant
for CD70 in activated CD4+ and CD8+ naïve T cells, Supplementary Fig. 18-19). From SCARlink
analysis, we found that IL2RA gene expression is strongly predicted by chromatin accessibility
(Spearman R² = 0.511) and is associated with the predicted enhancer tile chr10:6052500–
6053000, which covers rs61839660 across four cell types (Supplementary Fig. 20 and
Supplementary Data 9).

In comparison, analysis of lead variants uncovered the highest number of overlaps between IBD-
associated SNPs and CREs governing CCL20 and CXCR5. For instance, rs1811711, with a posterior
probability of 0.914 for being a causal variant for UC but not directly linked to CCL20 in a previous
fine mapping study6, was identified within a plausible CRE correlated with CCL20 expression in
activated Tfh cells and seven other cell types (Supplementary Fig. 21A-B). The association
between rs1811711 and CCL20 was further supported by eQTL data from whole blood
(Supplementary Fig. 21C). Similarly, rs630923, with a posterior probability of 0.820 for being a
causal variant for CD but not directly linked to CXCR5 in a previous fine mapping study70, was
identified within a plausible CRE correlated with CXCR5 expression in activated other CD4+
memory T cells and another five cell types (Supplementary Fig. 22A-B), although the association
between rs630923 and CXCR5 hasn’t been revealed by published bulk and single-cell eQTL
analyses.

These colocalization analyses provide not only specific loci but also cell-type specificity of IBD-
associated SNPs, compared to the general enrichment pattern revealed by CHEERS and S-LDSC.
Moreover, the colocalization of IBD-associated SNPs within cell type-specific CREs provides
compelling evidence for the association between genes regulated by CREs and SNPs located
within those CREs. The integration of SCARlink results further strengthens this analysis by
providing detailed predictions of gene-linked tiles and identifying potential enhancer elements
that regulate gene expression. This relationship is not fully elucidated by existing bulk and single-
cell eQTL analyses, highlighting the importance of in-depth, cell type-specific investigations with
greater precision. Together, these findings offer mechanistic insights into how genetic variants
may influence cell type-specific regulatory programs that contribute to IBD pathogenesis.

Cell type-specific colocalization of IBD-associated SNPs, DA peaks, peak-to-gene linkages, and
DE genes
To investigate the impact of both genetic and inflammatory factors (PGE2) on immune
dysregulation in IBD87, a colocalization analysis was performed to connect IBD-associated SNPs
with DA peaks and DE genes between T cell stimulation conditions with and without PGE2
(Supplementary Data 10). We hypothesized that ATAC peaks correlated with gene expression and
covering SNPs may be CREs that are potentially impacted by genetic susceptibility. Furthermore,
if these peaks and their linked genes are also altered upon PGE2 stimulation, this could suggest
a pathway that is potentially involved in IBD pathogenesis.

Out of the 151 IBD credible sets, 72 contain 633 IBD-associated SNPs involved in colocalization
events (Supplementary Data 11-12), marking a subset of these credible sets as essential for
exploring immune dysregulation in IBD influenced by both genetic and inflammatory factors.
Moreover, these events implicate a total of 325 genes, shedding light on critical pathways,
including the IL-2 signaling pathway, Th1 cytokine production, and NK cell differentiation
(Supplementary Data 13).

The greatest number of colocalization events cluster around PTGER4, highlighting 17 events
across seven T cell subsets. For example, rs10440635 is located within the ATAC peak chr5-
40490646-40491146, which is altered by PGE2 and correlated with the DE gene PTGER4. Such
associations are consistent across multiple T cell subsets, including activated CD4+ naïve, Th1,
CD8+ naïve, and resting CD4+ naïve T cells (Fig. 5A). These observations are further supported by
bulk data analysis from cells enriched for Th17 and Treg cells (Fig. 5B). As mentioned in the former
section, rs10440635 is located within a plausible CRE regulating PTGER4 expression across seven
distinct cell types. The additional findings that both the CRE's activity and PTGER4 expression are
modulated by PGE2 in a subset of these cells— four out of the seven T cell types— further
emphasize their plausible role in IBD pathogenesis. Similarly, we observed that rs9292781 and
rs1122433 are colocalized with a DA peak and the DE gene PTGER4 across four T cell subtypes
(Fig. 5C), with validation in bulk data (Fig. 5D). Collectively, our findings strikingly highlight
rs10440635, rs9292781, and rs1122433 for their pivotal roles in modulating PTGER4 expression
through the CREs where they reside across activated CD4+ naïve, Th1, CD8+ naïve, and resting
CD4+ naïve T cells, and this is further supported by validation with Hi-C data from GM12878 cell
lines88,89 (Supplementary Fig. 23A-C). This regulation, in conjunction with their responses to PGE2
stimulation, is particularly compelling, given PGE2's capacity to directly bind to the EP4 receptor
encoded by PTGER4, thereby delineating a mechanistic link to IBD pathogenesis.

Our colocalization analysis also highlights SNPs around SATB1 (10 events across four
subpopulations, Supplementary Fig. 24A-B), which regulate the pathogenic function of Th17 cells
in the experimental autoimmune encephalomyelitis model for multiple sclerosis90, and intestinal
epithelial function91. Additionally, SNPs around GPR183 (Supplementary Fig. 24C-D) and
RASGRP1 (Supplementary Fig. 24E) are highlighted. GPR183 expression is increased in the
presence of the CC allele of another IBD-associated SNP, rs9557195, and on lymphocytes
expressing chemokine receptors CCR6 or CCR9 and Th17 cells92. On the other hand, decreased
RASGRP1 in CD4+ T cells is correlated with active inflammation, while deletion of RASGRP1 causes
defective T lymphocyte development93. These findings are consistent with findings in our bulk
data. Interestingly, colocalization events for SNPs around CREM (a top up-regulated gene upon
PGE2 stimulation across various T cell subpopulations) are identifiable only in our single-cell
analyses of specific T cell subtypes but not bulk data (Supplementary Fig. 24F). These prioritized
IBD-associated SNPs around SATB1, GPR183, RASGRP1, and CREM also showed significant
associations with lymphocyte or eosinophil counts94, which could be substantially influenced by
PGE2 stimulation. This convergence of genetic and inflammatory factors underlines the complex
interplay in IBD pathogenesis.

Valuable resource for investigating other immune-mediated inflammatory diseases
Beyond prioritizing IBD credible sets and T cell subpopulations, our comprehensively
characterized T cell atlas can be generalized as a versatile resource to investigate other immune-
mediated inflammatory diseases. Similar to the enrichment of IBD GWAS SNPs in ATAC peaks,
we
examined 4,950 fine-mapped SNPs associated with 21 immune-mediated diseases, nearly all of
which are preferentially mapped to enhancers and promoters active in CD4+ T cells51, using
CHEERS. This analysis yielded enrichment P-values for each of the 500 pseudo-bulk sample
(Supplementary Fig. 25A) that were subsequently combined and adjusted as above
(Supplementary Fig. 25B).

We found that CD8+ memory T cells, both in resting and activated states, show striking
enrichment for disease-associated SNPs across multiple disorders, particularly psoriasis, Behçet’s
disease, and systemic lupus erythematosus (SLE). For example, we found a striking enrichment
of psoriasis-associated SNPs in ATAC peaks of CD8+ memory cells, with remarkably low
Bonferroni-adjusted P-values (Bonferroni-P) of 2.17x10-49 (resting state) and 1.37x10-45 (activated
state). This finding was not reported in the earlier fine mapping study51. CD8+ T cells have
subsequently been recognized as key effector cells and the predominant intra-epidermal T cells
in psoriatic plaques, and neutralization of CD8+ T cells prevents the development of psoriasis in
a xenograft model95. We also observed a notable enrichment of Behçet's disease-associated SNPs
in resting (Bonferroni-P = 4.18x10-43) and activated (Bonferroni-P = 7.86x10-39) CD8+ memory T
cell ATAC peaks. The frequency of CD8+ T cells is increased in patients with active Behçet's
disease and HLA-B51, and knock-out of ERAP1 in the presence of HLA-B51 alters CD8 T cell
immunity96.

In addition, we found that alopecia areata-associated SNPs are enriched in activated CD4+
regulatory T cells (Treg, Bonferroni-P = 1.33x10-37) ATAC peaks. Treg in hair follicles can facilitate
hair follicle regeneration by promoting hair follicle stem cell proliferation and differentiation97,
and Treg augmentation using low-dose IL-2 is an effective treatment98. We found that ankylosing
spondylitis-associated SNPs are enriched in resting (Bonferroni-P = 3.85x10-35) and activated
(Bonferroni-P = 7.27x10-30) CD8+ memory T cells, which can be primed by microbial peptides
presented by HLA-B27 and subsequently interact with self-peptides that HLA-B27 can bind99.
Furthermore, SLE-associated SNPs are enriched in activated CD8+ memory T cell ATAC peaks
(Bonferroni-P = 2.53x10-31). Defective effector CD8+ T cell function in SLE leads to failure to
remove autoreactive B cells that induce autoimmunity100.

Colocalization analysis of immune-mediated disease SNPs within cell type-specific CREs also
prioritize genetic loci and cell types crucial to disease pathogenesis. Focusing on the ten cell types
prioritized by CHEERS, psoriasis-associated SNPs were found to be enriched in peaks linked to
IFNLR1 (Supplementary Fig. 26A), as well as RUNX3 (Supplementary Fig. 27A) expression. In
psoriasis, Th17-secreted IL-29 mediates cutaneous anti-viral competence by activating IFNLR1101,
and we found that eight psoriasis-associated SNPs around IFNLR1 fall within five peaks correlated
with IFNLR1 expression across nine cell types (Supplementary Fig. 26B), involving both proximal
(e.g., rs7552167 in resting Th17 cells) and distal (e.g., rs11249213, rs4265380, and rs7414934 in
gamma-delta T cells) regulation of IFNLR1 expression. RUNX3 is important for CD8+ T cell
development and the differentiation of Th1 and Th17 cells102, and we found that ten psoriasis-
associated SNPs around RUNX3 fall within four peaks correlated with RUNX3 expression across
nine cell types (Supplementary Fig. 27B), involving both proximal (e.g., rs11249213, rs4265380,
and rs7414934 in resting Th1 cells) and distal (e.g., rs4090311 and rs78312791 in resting Th17
cells) regulation of RUNX3 expression. We also observed that Behcet’s disease-associated SNPs
are enriched in peaks correlated with KLRC4 expression (e.g., rs1841957 and rs1841958 in resting
CD8+ memory T cells, Supplementary Fig. 28A-B). Other notable associations include enrichment
of alopecia areata SNPs near IL21 in Th1 cells, ankylosing spondylitis SNPs near GPR65 in MAITs,
and SLE SNPs near BLK and CD44 in CD8+ memory T cells (Supplementary Fig. 29–32).

These results collectively highlight the broader applicability of our T cell atlas for elucidating
regulatory mechanisms across diverse immune-mediated diseases. In particular, the
convergence of genetic enrichment in CD8+ memory T cells across multiple conditions suggests
shared effector pathways that may underlie common immunopathology and provides a
foundation for therapeutic hypothesis generation across a spectrum of inflammatory disorders.

Methods
Enrichment and stimulation of peripheral blood T Cells
25 healthy young adult human study participants (ages 18-35 years) were enrolled with informed
consent in a research protocol compliant with all relevant ethical regulations, which was
approved by the University of Pittsburgh Institutional Review Board. T cells were enriched from
15 ml of whole blood following the MACSxpress Whole Blood Pan T Cell Isolation Kit, human
(Supplementary Data 14) protocol. T cells from the first seven study participants were cultured
immediately following isolation, while those from the remaining 18 study participants were
cryopreserved in CryoStor CS10 (BioLife Solutions) at a concentration of up to 10 million cells per
ml for subsequent processing.

For T cells from the first seven study participants (the first series of DOGMA-seq experiments),
four aliquots of one million cells each per study participant were cultured in 24-well cell culture
plates (Falcon) and stimulated overnight for 12 hours in 1 ml X-VIVO 15 Serum-free
Hematopoietic Cell Medium (Lonza) with 10 ul ImmunoCult Human CD3/CD28/CD2 T Cell
Activator (Supplementary Data 14), 50 ng/ml IL-1B (R&D Systems), and 50 ng/ml IL-23 (R&D
Systems), with either 1 uM Prostaglandin E2 (PGE2) (Sigma) or 3 ng/ml TGF-beta1 (R&D Systems),
or both, for a total of four stimulation conditions per participant. After the overnight incubation,
cells were collected into 15 ml conical tubes through a 30-um strainer and washed with up to 14
ml of PBS/0.2% BSA, centrifuged at 300 g for 10 minutes at 4°C, and resuspended in 50 ul
Pharmingen Stain Buffer (BSA) (BD Bioscience). An aliquot of each cell suspension was then
stained with ViaStain AOPI Staining Solution (Nexcelom Bioscience) and viability counts were
obtained using a Cellometer Auto 2000 Cell Viability Counter (Nexcelom Bioscience).

For T cells from the remaining 18 study participants (the second series of DOGMA-seq
experiments), cells from four participants at a time were thawed and two aliquots of one million
cells each per study participant were cultured in 24-well cell culture plates and stimulated
overnight for 12 hours in 1 ml X-VIVO 15 Serum-free Hematopoietic Cell Medium with 10 ul
ImmunoCult Human CD3/CD28/CD2 T Cell Activator, 50 ng/ml IL-1B, and 50 ng/ml IL-23, with or
without 1 uM PGE2, for a total of two stimulation conditions per participant. After overnight
incubation, cells were collected into 15 ml conical tubes through a 30-um strainer and washed
with up to 14 ml of PBS/0.2% BSA, centrifuged at 300 g for 10 minutes at 4°C, and resuspended
in 50 ul Pharmingen Stain Buffer. An aliquot of each cell suspension was then stained with
ViaStain AOPI Staining Solution and viability counts were obtained using a Cellometer Auto 2000
Cell Viability Counter.

Cell staining with barcoded antibodies
For the first series of DOGMA-seq experiments, after a 10-minute incubation with Human
TruStain FcX (BioLegend) on ice, up to 700,000 cells from each condition were stained with 1.4 ul
of a unique TotalSeq-A anti-human Hashtag antibody (Supplementary Data 14) in 70 ul for 30
minutes at 4°C on a Laminar Wash 16-well strip (Curiox Biosystems), for a total of 8 samples. Cells
were washed in the Laminar Wash Mini System (Curiox Biosystems) with 25 cycles at a flow rate
of 10 ul/s. Following cell collection and counting, up to 1 million cells from the different
hashtagged samples were pooled, centrifuged at 350g for 5 minutes at 4°C and resuspended in
50 ul. Fifty ul of TotalSeq-A Human Universal Cocktail, V1.0 (Supplementary Data 14-15), at a 0.5
dilution, were added to the cells and the total volume was split between 2 wells of a Laminar
Wash 16-well strip for a 30-minute incubation at 4°C. Cells were washed as previously described,
pooled into a 300 ul volume of PBS/0.2% BSA, filtered through a 40 um Scienceware Flowmi Cell
Strainer (SP Bel-Art) and counted.

For the second series of DOGMA-seq experiments, after a 10-minute incubation with Human
TruStain FcX on ice, approximately 300,000 cells from each condition were stained with 1 ul of a
unique TotalSeq-A anti-human Hashtag antibody in 50 ul for 30 minutes at 4°C on a Laminar Wash
16-well strip, for a total of 8 samples. Cells were washed in the Laminar Wash Mini System (Curiox
Biosystems) with 25 cycles at a flow rate of 10 ul/s. Following collection and cell count, up to 1.5
million cells from the different hashtagged samples were pooled, centrifuged at 350g for 5
minutes at 4°C and resuspended in 50 ul. Fifty ul of TotalSeq-A Human Universal Cocktail, V1.0,
at a 0.5 dilution, were added to the cells and the total volume was split between 2 wells of a
Laminar Wash 16-well strip for a 30-minute incubation at 4°C. Cells were washed as previously
described, pooled into a 300 ul volume of PBS/0.2% BSA, filtered through a 40 um Scienceware
Flowmi Cell Strainer and counted.

Single cell preparation for DOGMA-seq
The two series of DOGMA-seq experiments were conducted following the same protocol.
Specifically, cells were processed according to the CG000338 Chromium Next GEM Multiome
ATAC + Gene Expression Rev. D protocol (10X Genomics) for Gene Expression and ATAC library
construction using the Chromium Next GEM Single Cell Multiome ATAC + Gene Expression 4 rxn
kit (1000285, 10X Genomics) while modifications for DOGMA-seq were done following the “Cell
permeabilization with Digitonin (CDIG)” section of the DOGMA-seq protocol (NYGC Innovation
Lab).

Up to 1 million cells were centrifuged at 350g for 5 minutes at 4°C and permeabilized in 100 ul
chilled DIG lysis buffer (20 mM Tris-HCl pH 7.4, 150 mM NaCl, 3 mM MgCl , 0.0075% Digitonin
2
(Invitrogen) and 2 U/ul Protector RNase inhibitor (Roche Diagnostics)) for 1 min on ice, followed
by the addition of 1 ml of chilled DIG wash buffer (20 mM Tris-HCl pH 7.4, 150 mM NaCl, 3 mM
MgCl and 1 U/ul RNase inhibitor) and pipette mixing, before centrifugation at 500g for 5 minutes
2
at 4°C. The supernatant was discarded, the pellet was resuspended in a very small volume of
about 17 ul of diluted Nuclei Buffer (10X Genomics) and 5 ul were used for cell counting.
Following the CG000338 Chromium Next GEM Multiome ATAC + Gene Expression Rev. D protocol,
30,000 cells were processed for each of two transposition reactions, then added to a Barcoding
master mix and loaded into two wells of a Chip J and run on the Chromium Controller (10X
Genomics) for GEM generation.

During sample pre-amplification (step 4.2), ADT and hashtag oligonucleotide (HTO) additive
primers were added to a final concentration of 0.2 uM. After cleanup and elution in 100 ul, 30 ul
of the eluate were set aside for ADT and 10 ul were set aside for HTO i7 index PCR reactions that
were performed using Illumina TruSeq™ Small RNA RPIX (with custom indexes) and Illumina
TruSeq D70X_short primers, respectively, and the SI-PCR primer.

When generating the Gene Expression library in the Chromium Next GEM Multiome protocol, we
chose to use the Single Index Kit T Set A (1000213, 10X Genomics) instead of the Dual Index Kit
TT Set A (1000215, 10X Genomics).

Initial cDNA concentrations and ATAC library average sizes were determined by the TapeStation
High Sensitivity D5000 assay using HS D5000 Reagents (5067-5593, Agilent Technologies) and HS
D5000 ScreenTape (5067-5592, Agilent Technologies), and Gene Expression, ADT and HTO library
average sizes were determined by the TapeStation High Sensitivity D1000 assay using HS D1000
Reagents (5067-5585, Agilent Technologies) and HS D1000 ScreenTape (5067-5584, Agilent
Technologies) on the 2200 TapeStation System (Agilent Technologies). Library concentrations
were assessed by Qubit dsDNA HS assays (Q32851, Invitrogen). The Qubit concentration together
with the average size of fragments from the TapeStation assay were used to calculate the
estimated molarity of each library. The estimated molarity determined the dilutions necessary to
perform library quantification using the KAPA Library Quantification Kit, Universal (07960140001,
Roche Diagnostics) for each sample in order to obtain a more accurate molarity concentration of
the library pool for sequencing.

DOGMA-seq library pooling for sequencing
For the first series of DOGMA-seq experiments, libraries from DOGMA-seq experiments were
pooled and sequenced with libraries from other projects. The DOGMA-seq RNA libraries were
represented at either 0.81 nM or 2.14 nM, DOGMA-seq ADT libraries were represented at either
0.39 or 1.05 nM and HTO libraries at either 0.05 or 0.13 nM in two separate 10 nM pools of
libraries that were submitted to the UPMC Genome Center and sequenced using a NovaSeq 6000
S4 Reagent Kit v1.5 (200 cycles) (20028313, Illumina) on a NovaSeq 6000 sequencer (Illumina)
with a configuration 28/8/0/T151 (187 total cycles). DOGMA-seq ATAC libraries were represented
at approximately 2.5 nM and 5 nM in two separate 10 nM pools that were sequenced using using a
NovaSeq 6000 S1 Reagent Kit v1.5 (100 cycles) (20028319, Illumina) with a configuration
53/8/24/53 (138 total cycles).

For the second series of DOGMA-seq experiments, libraries from the DOGMA-seq experiments
were pooled and sequenced with libraries from other projects. The DOGMA-seq RNA libraries
were represented at approximately 1 nM, ADT libraries were represented at approximately 0.5
nM and HTO libraries at approximately 0.06 nM in two separate pools of libraries that were
submitted to the UPMC Genome Center and sequenced using a NovaSeq 6000 S4 Reagent Kit
v1.5 (200 cycles) on a NovaSeq 6000 sequencer with a configuration of either 28/10/10/151 (199
total cycles) or 28/8/0/151 (187 total cycles). DOGMA-seq ATAC libraries were represented at
either 1 nM or 2 nM in two separate pools that were sequenced using either a NovaSeq 6000 S4
Reagent Kit v1.5 (200 cycles) on a NovaSeq 6000 sequencer with a configuration of 100/8/24/100
(232 total cycles) or a NovaSeq 6000 S2 Reagent kit v1.5 (100 cycles) (20028316, Illumina) on a
NovaSeq 6000 sequencer with a configuration of 53/8/24/53 (138 total cycles).

All sequencing was performed on distinct samples.

Activation and stimulation of CD4+CD45RO+CD196+ peripheral blood T cells
39 healthy young adult human study participants (ages 18-35 years) were enrolled with informed
consent in a research protocol, which was approved by the University of Pittsburgh Institutional
Review Board. CD4+ T cells were enriched from 50 mL of peripheral blood using the RosetteSep
Human CD4+ T Cell Enrichment Cocktail (Supplementary Data 14) and sorted for
CD4+CD45RO+CD196+ cells using fluorescence-activated cell sorting. The sorted cells were then
incubated overnight at 1×106 cells/mL in serum-free X-VIVO-20 medium (Lonza) at 37°C, 5% CO2.
The following morning, the cells were stimulated with T Cell Activation/Expansion Kit, human
(Supplementary Data 14) at a ratio of 1 bead per 5 cells, along with IL-23 (50 ng/mL) and IL-1β
(50 ng/mL), with or without the addition of PGE2 (1 μM) for 6 hours.

Bulk ATAC and RNA sequencing
ATAC-seq libraries were prepared from 50,000 cells from each of the two stimulation conditions
per study participant according to published protocols with minor modifications107,108. To avoid
difficulties in estimating molarity for optimal cluster density during sequencing, size selection of
ATAC-seq libraries was performed using Pippin Prep with a 2% agarose gel cassette (Sage
Sciences), targeting fragments between 140 and 650 bp. Library quantification was assessed
using Bioanalyzer High Sensitivity DNA chip (Agilent Technologies) and KAPA Library
Quantification qPCR assays. Both ATAC-seq libraries from 33 study participants passed quality
control and were subsequently sequenced. HiSeq PE Rapid Cluster Kit v2 (PE–402–4002, Illumina)
and HiSeq Rapid SBS Kit v2 (200 Cycle) (FC–402–4021, Illumina) were used for sequencing of
library pools performed on HiSeq 2500 sequencers (Illumina) in Rapid Run Mode with 50 bp
paired-end sequencing (15 pairs of libraries). Later sequencing was performed on NextSeq 500
sequencers (Illumina) using High Output v2 (150 cycles) kits (15058247, Illumina), generating 75
bp paired-end sequencing, which was adopted for cost-effectiveness upon availability.

Following ATAC-seq library preparation, remaining cells from each treatment condition were
stored in QIAzol (QIAGEN) at -80°C for subsequent batch-extraction of total RNA. The quality and
concentration of extracted RNA were assessed using a Bioanalyzer and Quant-iT RiboGreen RNA
Assay Kit (Q33140, Invitrogen). For RNA-seq library preparation, the TruSeq Stranded Total RNA
Library Prep Gold (48 Samples) (20020598, Illumina) was used according to the manufacturer’s
protocol. Library pools containing up to 4 paired RNA-seq samples underwent 75 bp paired end
sequencing on NextSeq 500 sequencers using High Output v2 kits, resulting in a total of 38 pairs
of RNA samples being sequenced.

All sequencing was performed on distinct samples.

Raw DOGMA-seq sequencing data processing
The sequenced RNA and ATAC libraries were processed and aligned to GRCh38 human reference
genome using Cell Ranger ARC software (version 2.0.1). This software provided UMI counts, peak
counts, and fragments for each barcode. Joint cell calling was conducted to distinguish cells from
the background based on RNA UMIs per barcode and ATAC transposition events in peaks per
barcode, resulting in the generation of filtered UMI count and peak count matrices. The
sequenced ADT and HTO libraries were aligned using the kallisto | bustools workflow (kallisto
version 0.46.1, bustools version 0.39.3)109, resulting in the generation of raw UMI count matrices,
which were subsequently filtered according to the barcodes in the filtered UMI count and peak
count matrices for RNA and ATAC libraries. Cell demultiplexing was carried out based on the
filtered UMI count matrices for the HTO libraries using the function HTODemux() from the R
package Seurat (version 4.1.0)110, and only singlets remained for downstream analyses.

Peak calling, quality control, and integration of DOGMA-seq data
To obtain a consistent set of ATAC peaks for cells sequenced in different libraries, peaking calling
was performed for all singlets using the function CallPeaks() from the R package Signac (version
1.5.0)111, which invokes MACS2 (version 2.2.7.1)112 to call peaks based on fragments from all
libraries. Subsequently, counts in each called peak were quantified for all singlets using the
FeatureMatrix() function from Signac. A new peak count matrix based on the consistent set of
ATAC peaks across all libraries was generated to enable comparative analysis of chromatin
accessibility across libraries and mitigates batch effects between different libraries.

Filtered UMI count matrices for all RNA and ADT libraries were separately merged, resulting in
two UMI count matrices for RNA and ADT data containing only singlets. To ensure high-quality
data, quality control was conducted based on the three matrices containing only singlets.
Specifically, singlets were filtered with the following criteria applied: a minimum of 200 genes
detected, a minimum of 1,000 total RNA UMIs observed, a maximum of 20% mtRNA reads, a
maximum of 50 ADT UMIs from nine isotype controls observed, a minimum of 1,000 ATAC
fragments observed, and a minimum TSS enrichment score of 4. The remaining high-quality
singlets were subjected to downstream analyses.

For RNA data, regularized negative binomial regression (SCTransform)113 was used to normalize
the UMI count matrix, with glmGamPoi114 invoked to improve the speed and the variable
percent.mt regressed out. 3,000 highly variable genes were identified and used in principal
component analysis (PCA) dimensionality reduction. To mitigate batch effects, Harmony115 was
run on the first 50 components of PCA reduction, with distinct study participants serving as the
covariate and the maximum number of rounds set as 20. Stimulation conditions were not used
as a covariate to run Harmony due to our paired study design. Samples from the same study
participant but treated under different conditions were sequenced in the same batch, and
therefore, any observed variation between samples from the same study participant was
expected to be primarily due to biological differences rather than batch effects. The first 30
Harmony-adjusted components were used in Uniform Manifold Approximation and Projection
(UMAP)116 dimensionality reduction.

For ADT data, centered log ratio (CLR) transformation was used to normalize the UMI count
matrix. All 163 protein tags were set as variable features and used in PCA dimensionality
reduction. To mitigate batch effects, Harmony was run on the first 50 components of PCA
reduction, with distinct study participants and whether cryopreserved serving as the covariates
and the maximum number of rounds set as 20. Significant batch effects were observed at the
ADT level between fresh frozen and cryopreserved cells but not at the RNA or ATAC levels.
Therefore, whether or not cryopreserved served as an additional covariate to be adjusted only in
ADT data. The first 30 Harmony-adjusted components were used in UMAP dimensionality
reduction.

For ATAC data, Latent Semantic Indexing (LSI)117 was used to normalize and reduce the
dimensionality of the peak count matrix, with all the peaks selected as variable features. To
mitigate batch effects, Harmony was run on the first 50 components of LSI reduction, with
distinct study participants serving as the covariate and the maximum number of rounds set as 20.
The first 30 (excluding the first component due to correlation with sequencing depth) Harmony-
adjusted components were used in UMAP dimensionality reduction.

The Seurat function FindMultiModalNeighbors() was executed using the Harmony-adjusted
components from all three modalities as described above, followed by UMAP dimensionality
reduction based on the constructed weighted nearest neighbor (WNN) graph. In addition, a
similar process was performed, but only using components from the ADT and ATAC data.
Preliminary clustering was conducted using the SLM algorithm with a resolution of 0.2, based on
the WNN graphs or shared nearest neighbor (SNN) graphs for each modality.

Preliminary cell-type annotation was carried out using a reference-mapping approach. Anchors
at the RNA level were identified between DOGMA-seq data and Azimuth PBMC reference110.
Subsequently, cell type labels were transferred from the reference to DOGMA-seq data, and
clusters mainly consisting of non-T cells were filtered out, ensuring that only T cells were included
in the subsequent analyses. As a result of these quality control and filtering steps, 171,829 high-
quality T cells remained.

Peak calling was repeated for the remaining T cells, resulting in a peak count matrix containing
235,469 ATAC peaks. Additionally, the RNA UMI count matrix contained 36,601 genes, while the
ADT UMI count matrix contained 163 protein tags. The normalization, dimensionality reduction,
batch correction, and preliminary clustering steps described above were also repeated at each
modality for the remaining T cells. chromVAR (version 1.16.0)118 was run using the JASPAR2020
motif set to calculate the enrichment of chromatin accessibility at different TF motif sequences
in each cell.

Definition of activated and resting scores
To address the dynamic changes in T cells upon activation, we computed activated/resting scores
for each cell based on two lists of signature genes, which were derived from bulk RNA-seq data
profiling CD4+ T cells cultured with and without anti-CD3/CD28 beads. Cano-Gamez et al. isolated
naïve (CD25-CD45RA+CD45RO-) and memory (CD25-CD45RA-CD45RO+) CD4+ T cells from PBMCs,
cultured them separately with (activated) or without (resting) anti-CD3/anti-CD28 coated beads
for 16 hours or five days61. Bulk RNA-seq was performed on the cultured cells
(EGAD00001005291), and four lists of differentially expressed (DE) genes were identified for each
cell type and time point, comparing the activated and resting states.

Signature genes for the activated state were defined as DE genes with log2-fold change (LFC) > 1
and FDR-adjusted P-value (FDR) < 0.05, while signature genes for the resting state were defined
as those with LFC < -1 and FDR < 0.05. To ensure stringency, only the overlapping signature genes
between naïve and memory T cells across two time points were defined as the final lists of
signature genes for activated (containing 925 genes) and resting states (containing 1,761 genes),
which were used to calculate activated and resting scores based on RNA data using the Seurat
function AddModuleScore().

Cell type annotation of DOGMA-seq data
To overcome the computational challenge raised by the large DOGMA-seq dataset, cell type
annotation was conducted using a two-stage approach. First, a representative subset of 20,000
T cells was randomly down-sampled from the original 171,829 high-quality T cells. The
normalization, dimensionality reduction, and batch correction steps described above were
repeated at each modality for the representative subset. Clustering was conducted using the SLM
algorithm with a resolution ranging from 0.02 to 0.2, based on the WNN graphs or SNN graphs
for each modality. The R package clustree119 was used to determine the minimum resolution
needed to identify biologically meaningful clusters.

Based on the WNN graph combining ADT and ATAC data, five biologically meaningful clusters can
be easily identified. To identify subclusters within each cluster, iterative normalization,
dimensionality reduction, batch correction, and clustering steps were executed. Canonical
markers at RNA, ADT, and ATAC levels (listed in Supplementary Fig. 2A, C, D) were examined in
each subcluster, as well as the activated/resting scores. The subclusters within the five clusters
were subsequently annotated as distinct subtypes of T cells.

After identifying the 20 subtypes of T cells in the representative subset, the reference-mapping
approach was used to transfer cell type labels to the original 171,829 high-quality T cells. Label
transfer was conducted individually for each library (171,829 cells split by libraries), with
predicted labels merged at the end. Anchors were identified at RNA, ADT, and ATAC levels using
the annotated subset as a reference, resulting in three matrices of prediction scores for each cell
type based on the three modalities. To integrate the prediction scores across the three modalities,
a consensus strategy was utilized, where cells were assigned the labels with the highest
prediction scores across the three modalities.

To support future annotation efforts of Th1 and Th17 cells, we have developed epigenetic and
transcriptomic marker panels specific to these subsets, comparing them to all other CD4+
memory T cells (Supplementary Data 16-19).

Integrative analysis of in-house and public DOGMA-seq data
Mimitou et al. cultured PBMCs in complete medium with (STIM) or without (CTRL) anti-CD3/anti-
CD28 beads for 16 hours54, and the cells were subjected to DOGMA-seq (GSE156478). Following
the same procedure as above, 21,057 T cells were identified for integrative analysis. Counts in
the same set of 235,469 ATAC peaks were quantified for all cells based on the assumption that
the set of ATAC peaks in public DOGMA-seq data is expected to be a subset of the 235,469 ATAC
peaks.

A subset of 20,000 T cells was randomly down-sampled from the annotated in-house DOGMA-
seq data, remaining the previously calculated activated/resting scores. Raw RNA and ATAC data
from in-house and public DOGMA-seq data were integrated, with batch effects between the two
sources corrected by Harmony. Notably, ADT data were not used in the integration process, as
distinct panels of protein tags were added in the two experiments. The UMAP plots shown in
Supplementary Fig. 4 were based on ATAC data since those based on both RNA and ATAC data
or RNA data only were less distinct.

Analysis of DOGMA-seq data under different stimulation conditions
To compare the effects of the four distinct stimulation conditions, 84,345 cells from the first
DOGMA-seq experiment (7 study participants x 4 conditions) were extracted. The same
normalization, dimensionality reduction, and batch correction (using only distinct study
participants as the covariate) steps as described above were repeated at each modality, followed
by UMAP dimensionality reduction using components of both ADT and ATAC data. Next, the
activated/resting scores were recalculated for the extracted cells. Two-sided paired Wilcoxon
tests were also performed to compare the changes in cell type proportions upon different
stimulations. P

To focus on the effects of PGE2, 129,845 cells stimulated under conditions with and without PGE2
were extracted from the 25 study participants (25 study participants x 2 conditions), which were
subjected to a similar analysis as described above.

Cell type-specific differential accessibility (DA) and differential expression (DE) analysis
To compare the changes in chromatin accessibility, gene expression, and protein abundance
upon PGE2 stimulation (Act_IL1B_IL23_PGE2 vs. Act_IL1B_IL23), DA and DE analyses were
performed in each cell type using the pseudo-bulk method, which was expected to outperform
generic and specialized single-cell DE methods120. The R package Libra (version 1.0.0) was
modified to allow for paired samples when employing DESeq2 (version 1.34.0)’s likelihood ratio
test (LRT). Specifically, for each cell type, pseudo-bulk samples were generated by summing peak
counts and gene counts of all cells from the specific cell type and one of the 50 samples (25 study
participants x 2 conditions). The resulting P-values were adjusted for multiple testing using the
FDR method, and the significance level was defined as FDR-adjusted P-values < 0.05. For DA
analysis, ATAC peaks were annotated based on their top linked genes, or in cases where no linked
genes were identified, the nearest genes were used for annotation.

To summarize significant DE genes identified in each cell type, Ingenuity Pathway Analysis (IPA)121
was performed based on all available canonical pathways curated by IPA. For each cell type,
genes with FDR-adjusted P-values < 0.05 were used as inputs. Since the bulk data analysis
identified more than 8,000 genes with FDR-adjusted P-values < 0.05, we used an FDR-adjusted P-
value threshold of < 10-7 to select genes as inputs, in order to comply with IPA’s limitation of a
maximum of 8,000 genes.

To compare the results obtained from DOGMA-seq data with those from bulk data enriched for
Th17 and Treg cells, pseudo-bulk samples were generated by summing peak counts and gene
counts of all cells from the comparable cell types (activated and resting Th17 and CD4+ regulatory
T cells) and one of the 50 samples. The resulting pseudo-bulk samples merging four cell types
were subjected to a similar analysis as described above.

Raw bulk sequencing data processing
The sequenced ATAC libraries were trimmed and aligned to GRCh38 human reference genome
using Bowtie2 (version 2.4.5)122 in very sensitive mode. Reads with MAPQ scores less than 20,
aligned to chrM, and located in genomic blacklist regions annotated by ENCODE were filtered out,
as well as duplicated reads. Counts in the same set of 252,371 ATAC peaks as DOGMA-seq data
were quantified for all samples based on the assumption that the set of ATAC peaks in bulk data
is expected to be a subset of the 252,371 ATAC peaks.

The sequenced RNA libraries were trimmed, aligned to GRCh38 human reference genome using
STAR (version 2.7.6a)123, and quantified using RSEM (version 1.3.1)124 through the nf-core
(version 2.7.2)125 RNA-seq pipeline (version 3.4). The resulting matrices of gene and transcript
counted for 29,972 genes and 87,725 transcripts.

Analysis of bulk ATAC and RNA data
DA analysis was performed on ATAC data using DESeq2’s likelihood ratio test (LRT) in 66 samples
(33 study participants x 2 conditions) to compare the effects of PGE2 stimulation while allowing
for paired samples. Similarly, DE analysis was performed on RNA data using DESeq2’s likelihood
ratio test (LRT) in 76 samples (38 study participants x 2 conditions) to assess the effects of PGE2
stimulation while allowing for paired samples. The resulting P-values were adjusted for multiple
testing using the FDR method, and the significance level was defined as FDR-adjusted P-values <
0.05.

The correlations between effects obtained from bulk analysis and comparable single-cell analysis
were examined, and signals with consistent directions and FDR-adjusted P-values < 0.05 in both
bulk and sing-cell analyses were defined as replicated signals.

Peak-to-gene linkages were also identified in bulk data based on the 64 samples (32 study
participants x 2 conditions) subjected to both ATAC-seq and RNA-seq. Specifically, both ATAC and
RNA data were VST-normalized, and correlations between peaks and genes within a distance of
1M bps were assessed using the R package MatrixEQTL (version 2.3) with a linear model and no
covariates.

Germline variant calling from single-cell multimodal omics data
We used Monopogen126 to accurately detect individual genotypes at each genomic site using
single-cell DOGMA-seq data. Separate BAM files were generated for ATAC-seq and RNA-seq data
for each sample, enabling independent SNP calling from the two modalities. Genotype data
derived from ATAC-seq and RNA-seq were subsequently merged, and discordant calls between
the two datasets were excluded to maintain consistency. To ensure data integrity, duplicate
variants and flipped alleles were removed. Additionally, genotypes with sequencing depth below
a threshold of 10 reads were filtered to improve reliability.

To facilitate integration across samples, we consolidated genotype data from all individuals,
enabling SNP calling directly from the combined dataset to capture variants shared among the
samples. To further enhance the quality of the final genotype dataset, we applied a missingness
filter with a threshold of 0.5, retaining genotypes with sufficient data coverage across the
samples. Identified SNPs were then compared against the SNPs from 1000 Genomes Project
(1KGP)127 to identify shared variants.

After resolving discordances, the two datasets were merged. Duplicate variants and flipped
alleles were removed, and SNPs with missingness greater than 0.5 were filtered out, resulting in
a final set of 60,962 SNPs across 25 samples. Among the identified SNPs, the majority (68.0%)
were called from chromatin accessibility data, 13.9% were derived from gene expression data,
and 18.1% were called from both data modalities. A substantial overlap was observed between
the identified SNPs and those cataloged in the 1KGP127 , with overlap proportions of
approximately 93.3% for chromatin accessibility-derived SNPs and 95.6% for gene expression-
derived SNPs.

Allele-specific gene expression and chromatin accessibility
For the called genotype data, phasing of all biallelic common variants was performed using
shapeit4 (version 4.1.3)128 to distinguish SNPs mapping to maternal versus paternal alleles.
Selection of heterozygous sites was carried out using SelectVariants from GATK (version 4.5.0.0),
and allelic gene expression and peak reads at each heterozygous site were counted using GATK
ASEReadCounter129,130. To detect maximal allelic imbalance, the allelic peak reads for each
heterozygous SNP in each sample were summed. Only biallelic SNP sites with a total of at least
10 reads and a minimum of 2 reads for each allele were retained131,132. Finally, a binomial test
was employed to ascertain significant differential accessibility between alleles, considering the
alternate read counts as success trials and total read counts as trials. The threshold of FDR-
adjusted p < 0.05 was considered as significant.

To evaluate the concentration of IBD-related signals within ASoC/ASE SNPs, we compared these
with the P-value distributions from IBD GWAS summary statistics6. A quantile-quantile (QQ) plot
was then generated to visually evaluate the P-value distribution, aiding in the detection of
potential deviations from the expected uniform distribution. This was achieved by calculating the
expected and observed -log10 transformed P-values for each group, allowing us to visually
inspect for signals of true association as opposed to mere noise. To quantify the differences
between the observed and expected P-value distributions, we employed the Wilcoxon log-rank
sum test.

Enrichment of GWAS SNPs in ATAC peaks
The analysis assessed two sets of fine-mapped GWAS SNPs: the first combines 4,329 SNPs for IBD
from 151 credible sets identified in two studies6,70, and the second includes 4,950 SNPs for 21
autoimmune diseases51. IBD SNPs were further divided into six groups based on their specificity
to CD and UC, as well as whether they are lead variants in credible sets of SNPs. Genome
coordinates of two sets of GWAS SNPs were converted from GRCh37 to GRCh38 for consistency.
Variants in strong linkage disequilibrium (LD) with R2> 0.8 were then identified for each of the 27
lists of SNPs (21 autoimmune diseases + 6 groups of IBD SNPs) using CHEERS (master branch)71,
based on the European-ancestry GRCh38 LD files provided by CHEERS.

The peak count matrix containing 252,371 ATAC peaks was used to generate 500 pseudo-bulk
samples (20 cell types x 25 study participants) by summing peak counts of all cells from a specific
cell type and study participant. The resulting pseudo-bulk samples were then normalized using
CHEERS. For each list of SNPs, CHEERS computed the enrichment of SNPs in the 252,371 ATAC
peaks across the 500 pseudo-bulk samples, resulting in enrichment P-values per sample. Next,
the most conservative methods were used to combine enrichment P-values at the cell type level.
Specifically, individual P-values for each cell type were combined using Fisher’s method, resulting
in 20 combined P-values, which were further adjusted using Bonferroni’s method across the 20
cell types.

When considering stimulation conditions, a similar analysis was performed but confined to cells
under conditions with and without PGE2. In this analysis, 12 study participants were randomly
selected from the 25 study participants, and 480 pseudo-bulk samples (20 cell types x 2
conditions x 12 study participants) were generated by summing peak counts of all cells from a
specific cell type and one of the 12 selected study participants. Down-sampling was necessary
due to computational challenges, as the running time of CHEERS on 500 pseudo-bulk samples
had already reached the 6-day wall time of the computational cluster. To compare the
differences in enrichment between conditions with and without PGE2, individual P-values were
converted to Z-scores, which were compared between conditions using two-sided paired t-tests.

Heritability enrichment analysis
Heritability enrichment analysis through the Stratified Linkage Disequilibrium Score Regression
(S-LDSC) method allows for a nuanced understanding of how genetic variants contribute to
disease heritability133,134. S-LDSC operates by partitioning heritability across different genomic
annotations, enabling researchers to assess whether certain genetic markers or regions are more
significantly associated with a trait than would be expected by chance. This method relies on the
principle that genetic variants in linkage disequilibrium (LD) with causal variants will collectively
contribute to the trait's heritability. By correlating the degree of LD (measured as LD score) with
trait heritability, S-LDSC can estimate the per-SNP heritability attributable to specific genomic
features.

To further validate the enrichment of specific T cell subtypes in IBD GWAS, we applied the S-LDSC
method to partition heritability for IBD, CD and UC6,70 with 53 baseline annotations alongside
binary annotations specific to each T cell subtype75,76. In our analysis, we employed Wilcoxon
rank sum test within Seurat110 to identify the differentially accessible chromatin regions for each
cell type within the samples without PGE2 simulation. From this, we selected the top 20% of DA
peaks by P value for each cell type, extending the region of interest by adding 100-kb windows
to both sides of these peaks. SNPs located within these expanded regions were annotated as '1',
whereas all other SNPs received a '0' annotation. The 503 European individuals from the 1KGP
were included as our LD reference panel127,135. We set a rigorous Bonferroni correction threshold
at the P-value of 8.33 x 10-4, accounting for adjustments across 3 traits and 20 cell types.

Cell type-specific peak-to-gene linkages
To reduce the substantial noise in the correlative analysis of sparse single-cell data, the ArchR
(version 1.0.1) function addPeak2GeneLinks() was used to generate low-overlapping aggregates
of cells for each cell type, which were used to identify peak-to-gene linkages within a distance of
1M bps with FDR-adjusted P-values < 0.05. On average, the analysis identified 1,556,942 peak-
to-gene linkages (ranging from 391,031 to 3,393,196) for the 20 cell types.

To obtain more conservative results, cell type-specific peak-to-gene linkages were replicated in
pseudo-bulk analysis. Specifically, for each cell type, pseudo-bulk samples were generated by
summing peak counts and gene counts of all cells from the specific cell type and one of the 64
samples (7 study participants x 4 conditions + 18 study participants x 2 conditions). For each cell
type, pseudo-bulk ATAC and RNA data were normalized across the 64 samples using variance
stabilizing transformation (VST) from DESeq2136. Then, correlations between peaks and genes
within a distance of 1M bps were assessed using the R package MatrixEQTL137 with a linear model
and no covariates.

Cell type-specific peak-to-gene linkages identified using ArchR were retained as plausible linkages
for downstream analysis if they were replicated in pseudo-bulk analysis with consistent direction
and raw P-value < 0.05. On average, 101,972 plausible peak-to-gene linkages (ranging from
42,724 to 197,245) were retained for the 20 cell types, which correspond to approximately 6.5%
of peak-to-gene linkages initially discovered using ArchR.

The nearest genes for ATAC peaks were annotated using HOMER (version 4.11)138 based on the
hg38 (GRCh38) human reference genome. The top linked genes with the highest significance
were then assessed to determine whether they were also the nearest genes for specific peaks.

SCARlink analysis
To further investigate the effects of cis-regulatory elements and validate the results of our
previous ArchR peak-to-gene linking analysis using single-cell multi-omics data, we applied a
single-gene level model, SCARlink81, to identify potential cell type-specific enhancers of genes.
This method begins by dividing chromatin accessibility data into 500 bp tiles. For each gene,
regions spanning 250 kb upstream to 250 kb downstream of the gene body were selected to
extract tile-level chromatin accessibility features from scATAC-seq data.

SCARlink employs regularized Poisson regression to predict single-cell gene expression based on
single-cell chromatin accessibility. From the regression coefficients, cell type-specific Shapley
values are calculated for each tile to determine their relative importance. Tiles with high Shapley
values are often associated with enhancers, suggesting the presence of potential regulatory
elements within or near those regions. This approach thus enables the identification of putative
cell type-specific enhancers for individual genes.

For this analysis, we randomly selected 50,000 cells from our DOGMA-seq dataset and focused
on 5,000 HVGs. Genes with sparse expression (defined as having more than 90% zero values) were
excluded. For each retained gene, SCARlink models were trained, and Shapley values for
tiles within each gene region were computed and subsequently standardized using z-score
transformation. The significance of the difference in gene expression prediction with and without
a specific tile was assessed on held-out test cells using a paired Wilcoxon signed-rank test, with
FDR adjustment applied to account for multiple hypothesis testing. These computations were
performed across all 20 T cell subtypes to systematically evaluate the enhancer potential of
specific tiles.

Colocalization of GWAS SNPs and plausible peak-to-gene linkages
CHEERS provides a general assessment of the enrichment pattern of GWAS SNPs, narrowing
down the cell types worth exploring, but it may lack the granularity for each locus. A
complimentary analysis was performed to map the two sets of fine-mapped GWAS SNPs
described above to the 252,371 ATAC peaks, especially those linked with nearby genes. The set
of SNPs associated with autoimmune diseases was mapped to 668 ATAC peaks, for which a total
of 7,587 plausible cell type-specific peak-to-gene linkages were identified across the 20 cell types.
Similarly, the set of SNPs associated with IBD was mapped to 460 ATAC peaks, for which a total
of 4,877 plausible cell type-specific peak-to-gene linkages were identified across the 20 cell types.

To explore a specific disease of interest, the corresponding list of SNPs was mapped to peaks
linked with nearby genes, with a focus on the cell types prioritized by CHEERS. Specifically,
plausible cell type-specific peak-to-gene linkages in the prioritized cell types were extracted if the
peaks of these linkages were mapped by the disease-associated SNPs. The results were then
summarized in a gene-centric manner by counting the number of SNPs within linked peaks for a
specific gene across all prioritized cell types. The locus around the gene with the highest number
of SNPs mapping to its linked peaks across all prioritized cell types was selected for visualization.
Cell type-specific peak-to-gene linkages mapped by SNPs around the locus were also visualized
at the pseudo-bulk level. We also evaluated the overlap between SCARlink gene-linked tiles,
colocalized IBD SNPs, and plausible peak-to-gene linkages to explore their potential regulation
via enhancer elements.

Cell type-specific colocalization of IBD SNPs, DA peaks, peak-to-gene linkages, and DE genes
To investigate the potential relevance of ATAC peaks mapped by IBD SNPs upon PGE2 stimulation,
colocalization analysis was performed to establish a connection between IBD SNPs, DA peaks,
peak-to-gene linkages, and DE genes for each cell type.

The analysis began by identifying all linked genes for each DA peak in each cell type and examining
whether any of these linked genes were also DE genes. By doing so, linked DE genes for each DA
peak were identified through peak-to-gene linkages, and the consistency among the three types
of signals was checked. For example, if both a DA peak and its linked DE genes were up-regulated
upon PGE2 stimulation, a positive peak-to-gene linkage was expected for a consistent signal. DA
peaks consistently linked with at least one DE gene were hypothesized to be crucial upon PGE2
stimulation through the regulation of the linked DE genes. Next, the overlap between ATAC peaks
mapped by IBD SNPs and DA peaks consistently linked with at least one DE gene was extracted
and defined as the colocalization of IBD SNPs, DA peaks, peak-to-gene linkages, and DE genes.

The same colocalization analysis was also performed in bulk data and used as a replication for
some of the results identified in the single-cell analysis of the 20 cell types. Selected colocalization
signals were subsequently visualized using scatter plots to depict peak-to-gene linkages and
boxplots to depict DE genes and DA peaks, as well as a highlight for IBD SNPs falling within the
DA peaks.

Following the colocalization analysis, 72 out of the 151 IBD credible sets, encompassing 633 IBD-
associated SNPs participating in colocalization events, were prioritized, implicating a total of 325
genes. Subsequently, a Gene Ontology (GO) enrichment analysis139,140 was performed on these
325 genes to identify overrepresented pathways.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Review of the Manuscript\nComments on the Abstract\nThe abstract addresses important themes in immunology and IBD research, but it currently lacks clarity, coherence, and\nfocus. Several elements are missing or introduced too abruptly. Key points:\n• The research goal and motivation are not clearly stated. The rationale for introducing PGE2 is missing and leaves the\nreader wondering about its relevance.\n• The structure is somewhat disjointed. Findings are listed without a logical progression from background → research\nquestion → methods → main findings → broader implications.\n• Some findings are described too vaguely (e.g., \"PGE2 inhibits activation\" — which T cells? By what measure?) and others\ntoo densely.\n• The final sentence is overly promotional; it would be better to focus on insights gained rather than dataset availability.\n• Stylistically, long sentences and repetitive phrasing (e.g., \"we found,\" \"we identified,\" \"moreover\") make it harder to read\nsmoothly.\nSuggestions:\n• Clarify the goal and background, particularly around PGE2 and its relevance to IBD.\n• Follow a clear logical arc in the abstract.\n• Use more specific and interpretable language for the findings.\n• Tone down promotional language.\nComments on the Introduction\nThe introduction is largely appropriate until the paragraph beginning with \"Thus, to further elucidate.\" Some issues arise at\nthis point:\n• It is not clearly stated what types of T cells were used, nor how they were prepared.\n• DOGMA-seq is introduced abruptly without even a short definition. A brief one-sentence explanation could help (e.g.,\n\"DOGMA-seq enables simultaneous profiling of chromatin, transcriptomics, and surface proteins in single cells.\").\n• There is no transition explaining why single-cell multi-omics is necessary for this study. One linking sentence could greatly\nimprove clarity.\n• PGE2 is introduced without sufficient background on its relevance to T cell biology or IBD.\nSuggestions:\n• Add 1–2 bridging sentences clarifying why DOGMA-seq was necessary and why PGE2 was studied.\n• Soften overly enthusiastic language about the methods; let the value of the approach emerge naturally from the results.\nResults: Simultaneous trimodal measurement and annotation of human peripheral blood T cells\nThis section contains important information, but there are areas where the writing could be made much clearer:\n• The rationale for PGE2 stimulation remains unclear. Even a single sentence explaining that PGE2 is elevated in IBD and\nmodulates immune responses would help set expectations.\n• The switch from four stimulation conditions to two is confusing and should be explained more simply (\"Minimal effects were\nobserved with TGF-β, so subsequent experiments focused on PGE2 alone.\").\n• There is too much technical detail about DOGMA-seq in the Results section; this could be streamlined to keep the focus on\nbiological findings.\n• The relevance of the activated/resting scores could be better framed to explain biological meaning, not just statistical\nvalidation.\n• Some apparent paradoxes (e.g., overall suppression of activation markers yet upregulation of pathways like mTOR and\nSTAT4) are not explained.\nSuggestions:\n• Briefly frame PGE2 relevance upfront.\n• Explain design choices (4 vs 2 stimulation conditions) cleanly.\n• Streamline technical descriptions to avoid distracting from main findings.\n• Clarify the biological relevance of the activation/resting comparisons.\nResults: Allele-specific expression (ASE) and open chromatin (ASoC) analyses\nThe ASE/ASoC analysis is very valuable but quite difficult to follow:\n• The section opens with technical steps (germline variant calling, discordant variants) without stating why this matters — a\nshort framing sentence would help.\n• ASE and ASoC abbreviations are heavily used without reminding the reader what they mean.\n• Wilcoxon P-values are reported, but their biological meaning is not interpreted.\n• Excessive enumeration of SNP counts and condition-specific sites makes it hard to retain the main point.\n• The final, very important finding — that ASE and ASoC sites are enriched in IBD loci — is somewhat buried.\nSuggestions:\n• Briefly explain the goal of the ASE/ASoC quality control at the beginning.\n• Re-introduce ASE and ASoC definitions clearly.\n• Interpret the biological meaning of increased/decreased ASE and ASoC after PGE2 stimulation.\n• Emphasize the enrichment in IBD loci more clearly and earlier in the paragraph.\nResults: Enrichment of IBD GWAS SNPs in ATAC peaks\nThe general structure of this section is logical, but it can be improved:\n• \"Credible SNPs\" should be defined briefly when first mentioned.\n• CHEERS and S-LDSC are introduced without explaining their purposes in simple terms.\n• P-values are listed in dense groups without summarizing the overall patterns or trends.\n• Some key biological implications, particularly about activated T cell subtypes mediating IBD risk, are understated.\n• The relevance of PGE2 shifting enrichment patterns is interesting but could be better linked to clinical observations (e.g.,\nNSAID effects in Crohn's vs UC).\nSuggestions:\n• Add short, plain-language explanations of credible SNPs, CHEERS, and S-LDSC.\n• Summarize biological conclusions first before presenting detailed statistics.\n• Conclude the section with a clear statement about the clinical relevance of these findings.\nResults: Cell type-specific peak-to-gene linkages\nThis section is important but could be made easier to follow:\n• \"Peak-to-gene linkage\" should be briefly explained for readers less familiar with single-cell analysis.\n• The use of ArchR is introduced without explaining its role.\n• Again, the biological interpretation of the findings should come before listing numbers.\n• A short closing statement explaining why this matters (e.g., better annotation of non-coding disease variants) would\nstrengthen the section.\nSuggestions:\n• Add short intuitive definitions and interpretations before presenting numbers.\n• End the section by framing why these regulatory maps are important for understanding IBD\ngenetics.\nResults: SCARlink predictions\nWhile SCARlink is a strong addition to the study, it is not well introduced:\n• Shapley values are mentioned without explaining what they mean in this context.\n• The SCARlink method workflow is not clearly outlined.\n• The biological relevance to enhancer prediction for IBD is not made explicit.\n• Too many raw numbers are listed before any interpretation.\nSuggestions:\n• Introduce Shapley values with one sentence explaining the concept intuitively.\n• Briefly describe the goal of SCARlink (predicting enhancers linked to gene expression).\n• State the biological relevance upfront before listing numbers.\nResults: Colocalization of IBD-associated SNPs with CREs and DE genes\nThis section is better structured, but could still be made more efficient:\n• Lists of SNPs and genes (e.g., PTGER4, SATB1, GPR183) are too long and could be grouped more compactly.\n• The biological meaning and clinical relevance of the PTGER4 findings should be more clearly emphasized.\n• Redundancy about PTGER4 could be trimmed slightly.\nSuggestions:\n• Highlight the major discoveries first (especially PTGER4), and summarize other loci briefly.\n• Link the findings more explicitly to mechanisms of IBD pathogenesis.\nResults: Application to Other Immune-Mediated Diseases\nThe extension to other diseases is an important point but could be framed more cleanly:\n• Again, the lists of SNPs and genes are too long and detailed for the main text.\n• A short summary of the main diseases and key findings at the start would help.\n• The clinical relevance of connecting CD8+ memory T cells to psoriasis, Behçet’s, and SLE could be emphasized more\nclearly.\nSuggestions:\n• Summarize the major disease associations upfront.\n• Group minor findings together rather than listing them separately.\n• Mention broader translational relevance.\nDiscussion\nThe Discussion undersells the major scientific contributions of this study:\n• The opening is too focused on techniques (DOGMA-seq) rather than biological discoveries.\n• The important findings regarding PTGER4, Th17/Treg dynamics, and PGE2 in IBD pathogenesis are not highlighted as\nthey should be.\n• There is too much repetition of the Results without synthesis or proposing new models.\n• The potential clinical relevance of the findings is understated.\n• The final paragraph is too generic; it should more confidently state what new understanding the study provides.\nSuggestions:\n• Start with the major biological discoveries, not the techniques.\n• Frame a model linking genetic variants, chromatin accessibility, PGE2 signaling, and T cell state transitions in IBD.\n• Highlight translational implications, including PTGER4 as a potential therapeutic target.\n• End with a confident, forward-looking conclusion emphasizing the study's impact.", "reviewer_2": "The authors conducted an elegant multiomic analysis to explore the connection between inflammatory bowel disease (IBD)\nand genetic predisposition, applying a sophisticated computational workflow that, to the best of my knowledge, has not been\npreviously implemented. Their approach is undoubtedly innovative and contributes valuable insights.\nHowever, the manuscript's structure and flow present challenges, making it difficult for a broad readership to follow. The\nextensive use of technical terminology and abbreviations complicates the overall clarity, potentially obscuring the\nmanuscript's key message.\nDespite the novelty of the computational approach and the compelling nature of the study, the work lacks mechanistic\nvalidation of its findings in IBD, which appears to be the primary aim as suggested in the title. Given the complexity of IBD\npathogenesis, attributing T-cell inhibition solely to PGE2 seems overly simplistic. A broader perspective incorporating\nadditional mechanisms underlying T-cell-mediated immunity in IBD would enhance the work’s depth.\nI strongly recommend that the authors refine and streamline the text to improve accessibility while emphasizing the\nrelevance of their findings to IBD pathogenesis. Expanding the discussion beyond genetic associations to highlight key\ninflammatory and immune-mediated pathways—fundamental aspects of IBD—would significantly strengthen the\nmanuscript's clinical impact.", "reviewer_3": "As a reviewer coming from a clinical/transational background, I am not in a position to fully evaluate the technical aspects of\nthe DOGMA-seq single-cell multiomics platform. However, I find the characterization and depth of analysis regarding the\nnumber and subsets of T cells quite impressive and commendable.\nThat said, I believe the claim that this constitutes a model of IBD inflammation—based solely on in vitro stimulation of T cells\nwith IL-1β, IL-23, and either PGE2, TGF-β, or both—appears overly simplistic. While the approach clearly identifies key\nresponse elements to these specif-ic immune mediators, it does not recapitulate the complex immunopathology of IBD as\nseen in patients.\nIn line with this, I would have appreciated data derived from the DOGMA-seq platform on mucosal T cells obtained from IBD\npatients, both during remission and active disease, or from other immune-mediated inflammatory diseases (IMIDs), to\nstrengthen translational relevance.\nFinally, I would encourage the authors to reflect on whether their DOGMA-seq platform data provide any insights into the\npotential relationship between PGE2 signaling and non-response to anti-TNF therapy, an important and clinically relevant\narea that may benefit from the mechanistic resolution their platform enables."}, "author_response": "POINT-BY-POINT RESPONSE TO REVIEWER COMMENTS\nReviewer #1:\nRemarks to the Author:\nComments on the Abstract\nThe abstract addresses important themes in immunology and IBD research, but it currently lacks clarity,\ncoherence, and focus. Several elements are missing or introduced too abruptly. Key points:\n• The research goal and motivation are not clearly stated. The rationale for introducing PGE2 is\nmissing and leaves the reader wondering about its relevance.\n• The structure is somewhat disjointed. Findings are listed without a logical progression from\nbackground → research question → methods → main findings → broader implications.\n• Some findings are described too vaguely (e.g., \"PGE2 inhibits activation\" — which T cells? By\nwhat measure?) and others too densely.\n• The final sentence is overly promotional; it would be better to focus on insights gained rather\nthan dataset availability.\n• Stylistically, long sentences and repetitive phrasing (e.g., \"we found,\" \"we identified,\"\n\"moreover\") make it harder to read smoothly.\nSuggestions:\n• Clarify the goal and background, particularly around PGE2 and its relevance to IBD.\n• Follow a clear logical arc in the abstract.\n• Use more specific and interpretable language for the findings.\n• Tone down promotional language.\nR: We thank the reviewer for the insightful comments about our abstract. After consulting with the editor,\nwe have expanded the abstract beyond the 150 words recommended in the “Brief guide for submission to\nNature Communications” document to clarify the study goals and rationale and reorganize the structure\nso that it follows a more coherent progression from background and motivation to methods, main findings,\nand broader implications. The revised abstract explicitly highlights the relevance of PGE2 and the other\ninflammatory mediators in the T cell stimulation cocktails to IBD and T cell biology, and it provides a\nmore precise overview of methods and major results, including specific T cell subsets that are most\naffected by PGE2. We have removed overly promotional language and focused on presenting key\nbiological insights. We believe that the revised abstract better conveys the scope and significance of our\nwork.\nComments on the Introduction\nThe introduction is largely appropriate until the paragraph beginning with \"Thus, to further elucidate.\"\nSome issues arise at this point:\n• It is not clearly stated what types of T cells were used, nor how they were prepared.\n• DOGMA-seq is introduced abruptly without even a short definition. A brief one-sentence\nexplanation could help (e.g., \"DOGMA-seq enables simultaneous profiling of chromatin, transcriptomics,\nand surface proteins in single cells.\").\n• There is no transition explaining why single-cell multi-omics is necessary for this study. One\nlinking sentence could greatly improve clarity.\n• PGE2 is introduced without sufficient background on its relevance to T cell biology or IBD.\nSuggestions:\n• Add 1–2 bridging sentences clarifying why DOGMA-seq was necessary and why PGE2 was\nstudied.\n• Soften overly enthusiastic language about the methods; let the value of the approach emerge\nnaturally from the results.\nR: In the revised Introduction, we explicitly state that we studied T cells isolated from the peripheral blood\nof healthy human donors with ages 18-35 years, an age range similar to the peak ages of IBD onset in\nadults. The relevance of PGE2 and the other inflammatory mediators in the T cell stimulation cocktails to\nT cell biology and IBD is discussed. DOGMA-seq is introduced as a single-cell trimodal omics technique\nincorporating simultaneous assays for gene expression, cell surface protein epitope antibody-derived tags,\nand the assay for transposase accessible chromatin within individual cells, enabling higher resolution\nannotation of specific T cell subtypes and more robust dissection of complex regulatory networks in\ncomparison to single modal or bimodal assays. Overly enthusiastic language about the methods has been\ndeleted.\nResults: Simultaneous trimodal measurement and annotation of human peripheral blood T cells\nThis section contains important information, but there are areas where the writing could be made much\nclearer:\n• The rationale for PGE2 stimulation remains unclear. Even a single sentence explaining that PGE2\nis elevated in IBD and modulates immune responses would help set expectations.\n• The switch from four stimulation conditions to two is confusing and should be explained more\nsimply (\"Minimal effects were observed with TGF-β, so subsequent experiments focused on PGE2\nalone.\").\n• There is too much technical detail about DOGMA-seq in the Results section; this could be\nstreamlined to keep the focus on biological findings.\n• The relevance of the activated/resting scores could be better framed to explain biological meaning,\nnot just statistical validation.\n• Some apparent paradoxes (e.g., overall suppression of activation markers yet upregulation of\npathways like mTOR and STAT4) are not explained.\nSuggestions:\n• Briefly frame PGE2 relevance upfront.\n• Explain design choices (4 vs 2 stimulation conditions) cleanly.\n• Streamline technical descriptions to avoid distracting from main findings.\n• Clarify the biological relevance of the activation/resting comparisons.\nR: We thank the reviewer for the thoughtful comments and suggestions regarding the Results section. We\nhave carefully revised the text to improve clarity and to better explain the rationale, design decisions, and\nbiological interpretation of our findings.\n1. The relevance of PGE2 and the other inflammatory mediators in the T cell stimulation cocktails\nto T cell biology and IBD is discussed in the Introduction.\n2. We have clarified the stimulation design. Specifically, we now explain that we initially\nevaluated four stimulation conditions but observed minimal effects of TGF-β, so we focused\nsubsequent experiments on PGE2 by comparing two conditions with or without PGE2.\n3. To maintain focus on biological insights, we have streamlined the technical details in the Results\nsection. Methodological aspects such as the calculation of activated/resting scores, pseudo-bulk\nconstruction, and DOGMA-seq data processing are now described in the Methods section only.\nThis allows the Results to more directly emphasize key biological findings.\n4. In response to the suggestion for better framing of the activated/resting comparisons, we now\ndescribe the biological relevance of this dichotomy. Specifically, we show that while PGE2\nbroadly suppresses T cell activation across multiple subtypes, including CD4+ naïve, CD4+\nregulatory, CD4+ memory, CD8+ naïve, CD8+ memory, and MAIT cells, it also induces\nupregulation of numerous genes in resting T cells, including CREM, STAT4, and PDE4D, along\nwith activation of Protein Kinase A (PKA), mTOR, EIF2, and IL-23–related signaling pathways.\nIn contrast, these same pathways are down-regulated in activated cells, particularly in Th1, CD4+\nmemory, CD8+ memory, and MAIT cells.\nWe now provide a mechanistic interpretation of this pattern. PGE2 may initiate a transient\nimmunomodulatory program via the cAMP–PKA axis, while also engaging a negative feedback\nloop through upregulation of PDE4D (which degrades cAMP) and downregulation of PTGER4\nand the EP4 receptor for PGE2 that PTGER4 encodes, preventing prolonged immunosuppression.\nSimultaneously, we observe upregulation of STAT4 related pathways in resting T cells, suggesting\nthat resting T cells may be transcriptionally primed for cytokine-driven differentiation.\nThe suppression of T cell activation and upregulation of pathways such as mTOR and STAT4\nreflect dual and context-specific roles of PGE2. In resting cells, these pathways may promote\nmetabolic readiness and selective transcriptional priming, whereas in activated cells, they are\ndown-regulated, possibly due to exhaustion, desensitization, or dominance of TCR-driven\nprograms.\nAdditionally, we now include new observations related to Th17 cells, where we detect a\nnoncanonical transcriptional response. Upon PGE2 stimulation, CREB1, NFKB1, and other\ncAMP-related genes are upregulated in resting Th17 cells. Chromatin accessibility also increases\nnear the NFKB1 locus, suggesting transcriptional remodeling. These findings are consistent with\nprior murine studies showing that PGE2–EP2/EP4 signaling can synergize with IL-23 to enhance\nTh17 pathogenic programs via CREB1, STAT3, and NF-κB1. We interpret this as further evidence\nof the flexible, context-dependent reprogramming induced by PGE2.\nResults: Allele-specific expression (ASE) and open chromatin (ASoC) analyses\nThe ASE/ASoC analysis is very valuable but quite difficult to follow:\n• The section opens with technical steps (germline variant calling, discordant variants) without\nstating why this matters — a short framing sentence would help.\n• ASE and ASoC abbreviations are heavily used without reminding the reader what they mean.\n• Wilcoxon P-values are reported, but their biological meaning is not interpreted.\n• Excessive enumeration of SNP counts and condition-specific sites makes it hard to retain the\nmain point.\n• The final, very important finding — that ASE and ASoC sites are enriched in IBD loci — is\nsomewhat buried.\nSuggestions:\n• Briefly explain the goal of the ASE/ASoC quality control at the beginning.\n• Re-introduce ASE and ASoC definitions clearly.\n• Interpret the biological meaning of increased/decreased ASE and ASoC after PGE2 stimulation.\n• Emphasize the enrichment in IBD loci more clearly and earlier in the paragraph.\nR: We thank the reviewer for the constructive comments about the ASE/ASoC section. We have revised\nthe text substantially to improve clarity, structure, and biological interpretation.\n1. We now begin the section with a concise statement outlining the purpose of the ASE and ASoC\nanalyses and the rationale for initial variant-level quality control. We explain that we performed\ngermline variant calling separately from the gene expression and chromatin accessibility\nmodalities of DOGMA-seq, allowing us to infer sample-matched genotypes without requiring\nexternal genotyping arrays, and then filtered discordant variants between RNA- and ATAC-\nderived SNP calls, flipped alleles and duplicates, and variants with high missingness to arrive at\nhigh confidence genotype calls for accurate and allele-consistent integration of chromatin\naccessibility and gene expression data. This provides biological context for the technical\ndescription that follows. We have moved all details regarding proportions and numbers of\nidentified SNPs except the total number of high-confidence SNPs that remained to the Methods\nsection.\n2. To improve readability, after defining abbreviations for allele-specific expression (ASE) and\nopen chromatin (ASoC) when first mentioned in the Introduction, we have reintroduced these\ndefinitions again in the title of the Results subsection where ASE and ASoC analyses are presented.\n(We have also re-introduced abbreviations for differential chromatin accessibility (DA) and gene\nexpression (DE) a few times throughout the manuscript, and we have provided an additional table\nfor abbreviations in the supplementary files.)\n3. We have added biological interpretation of the Wilcoxon test results. Specifically, we highlight\nthat the increased number of significant ASE and ASoC sites in resting versus activated T cells\nafter exposure to PGE2 suggests enhanced allele-specific regulation in the resting state, consistent\nwith our previous findings that PGE2 suppresses activation programs while maintaining or even\nenhancing transcriptional and epigenetic plasticity in resting T cells. Conversely, the reduction of\nASE and ASoC sites in activated cells under PGE2 stimulation suggests an inhibitory effect on\ntranscriptional complexity or regulatory heterogeneity in an already activated state. These\nbiological interpretations are now stated clearly alongside the statistical comparisons.\n4. In response to the concern about excessive enumeration, we have reduced the number of exact\nSNP counts presented in the main text and moved detailed breakdowns to the Supplementary\nTables and Figures. The main paragraph now focuses on the key patterns and biological\nimplications of stimulation-specific versus shared allele-specific regulation.\n5. Finally, and importantly, we have revised the structure of the section so that the enrichment of\nASE and ASoC sites in IBD GWAS loci, a central result of this analysis, is presented earlier and\nmore prominently. We now emphasize this finding as a major conclusion of the section,\nhighlighting its relevance for understanding the genetic architecture of IBD and the stimulation-\ndependent regulatory landscape of human T cells.\nWe hope these changes improve the clarity and impact of this section in line with the reviewer’s helpful\nsuggestions.\nResults: Enrichment of IBD GWAS SNPs in ATAC peaks\nThe general structure of this section is logical, but it can be improved:\n• \"Credible SNPs\" should be defined briefly when first mentioned.\n• CHEERS and S-LDSC are introduced without explaining their purposes in simple terms.\n• P-values are listed in dense groups without summarizing the overall patterns or trends.\n• Some key biological implications, particularly about activated T cell subtypes mediating IBD\nrisk, are understated.\n• The relevance of PGE2 shifting enrichment patterns is interesting but could be better linked to\nclinical observations (e.g., NSAID effects in Crohn's vs UC).\nSuggestions:\n• Add short, plain-language explanations of credible SNPs, CHEERS, and S-LDSC.\n• Summarize biological conclusions first before presenting detailed statistics.\n• Conclude the section with a clear statement about the clinical relevance of these findings.\nR: We thank the reviewer for the constructive feedback. In the revised manuscript, we have made several\nchanges to improve clarity, interpretation, and clinical relevance in this section:\n1. We now define “IBD credible set SNPs” at the point they are first mentioned. We explain that\nthese SNPs were prioritized based on Bayesian fine-mapping approaches that account for linkage\ndisequilibrium and effect size, and are comprised of the SNPs in each of 151 fine-mapped IBD\nloci which have a greater than 95% probability of containing the causal variant.\n2. To aid understanding of the analytical tools used, we have added concise explanations for both\nCHEERS and S-LDSC. We describe CHEERS as a method that detects enrichment of disease-\nassociated SNPs in regulatory regions by accounting for the quantitative, cell-state-specific\nactivity of chromatin elements. We describe our application of S-LDSC as a heritability-based\nmethod that tests whether GWAS signal is enriched in specific annotations across cell types.\n3. We have restructured the paragraph to present biological conclusions first, followed by\nsupporting statistics. For example, we now summarize that IBD SNPs are most strongly enriched\nin activated Th17, Th1, and CD4+ Treg cells, and that CD-associated SNPs show stronger\nenrichment in activated states than UC-associated SNPs, before detailing the exact P-values.\n4. We have expanded our interpretation of biological trends, emphasizing the prominent role of\nactivated CD4+ T cell subsets (Th1, Th17, CD4+ Treg) across both CHEERS and S-LDSC results,\nwhich supports their key involvement in IBD pathogenesis.\n5. To address the reviewer’s final point, we have revised the last paragraph of the section to\nexplicitly link PGE2-dependent enrichment pattern shifts to clinical observations that NSAID use,\nwhich is known to reduce endogenous PGE2 levels, is associated with increased relapse risk in\nCD but not UC. We state that reduced regulatory flexibility at CD risk loci could increase\nvulnerability to PGE2 depletion during NSAID exposure, whereas the comparatively greater\nplasticity of UC loci may help buffer against such effects.\nWe believe these revisions significantly enhance the clarity, biological insight, and translational relevance\nof the section and hope they adequately address the reviewer’s concerns.\nResults: Cell type-specific peak-to-gene linkages\nThis section is important but could be made easier to follow:\n• \"Peak-to-gene linkage\" should be briefly explained for readers less familiar with single-cell\nanalysis.\n• The use of ArchR is introduced without explaining its role.\n• Again, the biological interpretation of the findings should come before listing numbers.\n• A short closing statement explaining why this matters (e.g., better annotation of non-coding\ndisease variants) would strengthen the section.\nSuggestions:\n• Add short intuitive definitions and interpretations before presenting numbers.\n• End the section by framing why these regulatory maps are important for understanding IBD\ngenetics.\nR: We appreciate the reviewer’s thoughtful comments and suggestions regarding the clarity and\ninterpretability of this important section. In the revised manuscript, we have made several changes to\naddress these points:\n1. We have added a brief, intuitive definition of peak-to-gene linkages early in the section,\nexplaining that they are correlations between chromatin accessibility peaks and expression levels\nof nearby genes, thereby mapping candidate cis-regulatory elements (CREs) to their putative gene\ntargets.\n2. We clarified the role of ArchR as a computational framework for inferring peak-to-gene linkages\nfrom single-cell multi-omic data.\n3. To improve readability, we restructured the paragraph so that biological insights appear before\nthe detailed numerical results. Specifically, we first highlight that peak-to-gene linkages often\nreveal regulatory relationships beyond nearest-gene assignments, and that naïve T cells tend to\nhave fewer such linkages despite their abundance, indicating lower regulatory complexity in their\nresting states.\n4. Finally, we added a closing sentence emphasizing the significance of these regulatory maps for\nIBD genetics. We state that identifying cell type-specific peak-to-gene linkages enables more\naccurate annotation of non-coding disease variants, and provides mechanistic hypotheses linking\nGWAS loci to gene regulation in relevant T cell subsets.\nWe believe these revisions improve the accessibility and biological interpretation of the section and thank\nthe reviewer again for helping to strengthen this part of the manuscript.\nResults: SCARlink predictions\nWhile SCARlink is a strong addition to the study, it is not well introduced:\n• Shapley values are mentioned without explaining what they mean in this context.\n• The SCARlink method workflow is not clearly outlined.\n• The biological relevance to enhancer prediction for IBD is not made explicit.\n• Too many raw numbers are listed before any interpretation.\nSuggestions:\n• Introduce Shapley values with one sentence explaining the concept intuitively.\n• Briefly describe the goal of SCARlink (predicting enhancers linked to gene expression).\n• State the biological relevance upfront before listing numbers.\nR: We thank the reviewer for the constructive comments. In the revised manuscript, we clarify the goal of\nSCARlink at the beginning of the section by explicitly stating that it predicts enhancer-gene relationships\nusing single-cell gene expression and chromatin accessibility data. We also introduce Shapley values as a\nmeasure of the relative importance of fixed-size genomic tiles in explaining the expression of a given gene.\nTo improve readability, we mention the biological rationale for applying SCARlink, to augment\nunderstanding of cis-regulatory mechanisms and annotation of enhancer elements relevant to IBD in T\ncells, earlier in the section. Finally, we have restructured the presentation of the numerical results to\nemphasize interpretation first, followed by supporting statistics. These revisions collectively improve the\nclarity, accessibility, and biological framing of the SCARlink analysis.\nResults: Colocalization of IBD-associated SNPs with CREs and DE genes\nThis section is better structured, but could still be made more efficient:\n• Lists of SNPs and genes (e.g., PTGER4, SATB1, GPR183) are too long and could be grouped\nmore compactly.\n• The biological meaning and clinical relevance of the PTGER4 findings should be more clearly\nemphasized.\n• Redundancy about PTGER4 could be trimmed slightly.\nSuggestions:\n• Highlight the major discoveries first (especially PTGER4), and summarize other loci briefly.\n• Link the findings more explicitly to mechanisms of IBD pathogenesis.\nR: We thank the reviewer for the helpful feedback. In response:\n1. We have streamlined the lists of SNPs and genes throughout this section by consolidating\noverlapping observations and presenting the PTGER4 locus and associated variants as a single\ncoherent paragraph. Other loci such as LNPEP, TNFRSF14, IL2RA, CCL20, and CXCR5 are now\ndescribed more succinctly in grouped summaries to maintain readability while preserving key\nbiological insights, such as predicted enhancer activity, peak-to-gene linkages, and supporting\nevidence from eQTL and fine-mapping studies.\n2. We agree that the biological and clinical relevance of PTGER4 deserved clearer emphasis. We\nhave revised the text to make our findings more interpretable and mechanistically meaningful.\nSpecifically, we now highlight that:\n2.1 PTGER4 expression is strongly predicted by chromatin accessibility (Spearman R² =\n0.244 from SCARlink analysis);\n2.2 Multiple IBD SNPs colocalize with enhancer regions linked to PTGER4 across diverse\nT cell subsets;\n2.3 These findings converge on the PTGER4 locus as a genetically and epigenetically\nregulated target in T cells, strongly implicating PGE2–EP4 signaling in IBD pathogenesis;\n2.4 The SNP rs10440635 is located in an IBD credible set, was previously implicated in\nankylosing spondylitis, and regulates PTGER4 in synovial biopsies, adding cross-disease\nrelevance.\n3. To address redundancy, we have removed repetitive mentions of colocalization results and\nmerged previously separated pieces of evidence (e.g. SCARlink predictions) into integrated\nsummaries for each key region.\n4. Finally, we have added a concluding sentence that explicitly links these colocalization results\nto mechanistic insights, emphasizing how genetic variants may influence cell type-specific\nregulatory programs that contribute to IBD pathogenesis.\nResults: Application to Other Immune-Mediated Diseases\nThe extension to other diseases is an important point but could be framed more cleanly:\n• Again, the lists of SNPs and genes are too long and detailed for the main text.\n• A short summary of the main diseases and key findings at the start would help.\n• The clinical relevance of connecting CD8+ memory T cells to psoriasis, Behçet’s, and SLE could\nbe emphasized more clearly.\nSuggestions:\n• Summarize the major disease associations upfront.\n• Group minor findings together rather than listing them separately.\n• Mention broader translational relevance.\nR: We thank the reviewer for the helpful suggestions to improve clarity and structure. In the revised\nmanuscript, we have substantially reorganized this section to highlight the major findings more clearly.\nSpecifically:\n1. We now open the section with a concise summary of the main disease associations, emphasizing\nthe prominent enrichment of psoriasis-, Behçet’s-, and SLE-associated SNPs in ATAC peaks of\nCD8⁺ memory T cells.\n2. To enhance readability, we have grouped secondary associations (e.g., those for alopecia areata,\nankylosing spondylitis) into a more compact paragraph, while retaining all corresponding\nSupplementary Figures.\n3. We have clarified the clinical relevance of the CD8⁺ memory T cell findings by explicitly linking\nthe genetic enrichment patterns to known effector functions and disease mechanisms in psoriasis,\nBehçet’s, and SLE.\n4. Finally, we have revised the concluding paragraph to more directly emphasize the broader\ntranslational value of our atlas in illuminating shared regulatory mechanisms and informing future\ntherapeutic hypotheses across immune-mediated diseases.\nImportantly, we preserved all Supplementary Figures to retain the full detail for interested readers, while\nimproving the narrative focus in the main text as recommended. We hope these revisions address the\nreviewer’s concerns.\nDiscussion\nThe Discussion undersells the major scientific contributions of this study:\n• The opening is too focused on techniques (DOGMA-seq) rather than biological discoveries.\n• The important findings regarding PTGER4, Th17/Treg dynamics, and PGE2 in IBD pathogenesis\nare not highlighted as they should be.\n• There is too much repetition of the Results without synthesis or proposing new models.\n• The potential clinical relevance of the findings is understated.\n• The final paragraph is too generic; it should more confidently state what new understanding the\nstudy provides.\nSuggestions:\n• Start with the major biological discoveries, not the techniques.\n• Frame a model linking genetic variants, chromatin accessibility, PGE2 signaling, and T cell state\ntransitions in IBD.\n• Highlight translational implications, including PTGER4 as a potential therapeutic target.\n• End with a confident, forward-looking conclusion emphasizing the study's impact.\nR: We thank the reviewer for the thoughtful and constructive feedback. In response, we have substantially\nrevised the Discussion section to better highlight the key biological discoveries and enhance the narrative\nsynthesis of our findings. Specifically:\n1. We have revised the opening paragraph to center the discussion on the biological context and\nmotivation, rather than beginning with a technical description of DOGMA-seq.\n2. We more prominently feature PTGER4 and its dual role as both a dynamically regulated\ntranscriptional target and a genetically implicated node in the PGE2 signaling pathway. We\nemphasize its potential as a translationally actionable therapeutic target in IBD, aligning with the\nreviewer’s request to highlight clinical relevance.\n3. In the revised Discussion, we propose a model in which PGE2 modulates T cell fate transitions\nin an activation-state-dependent manner, intersecting with inherited susceptibility through\nchromatin remodeling and gene regulation. We further describe how PGE2 selectively represses\neffector programs in activated cells while priming resting cells, offering a novel perspective on\nhow IBD-relevant microenvironment and genetic factors interact in shaping T cell responses\nduring disease pathogenesis.\n4. To reduce repetition, we have streamlined descriptive results and refocused on conceptual\ninsights, including the regulatory convergence between responses to PGE2, chromatin\naccessibility, and GWAS SNP enrichment in Th17, Th1, and Treg cells.\n5. Finally, we have strengthened the concluding paragraph to more confidently articulate the\nmechanistic model proposed by this study. Rather than ending with a generic summary, we now\nexplicitly state how our multi-omic atlas enables identification of disease-relevant regulatory\nmechanisms and offers a path forward for precision immunomodulation in IBD and other immune-\nmediated disorders.\nWe believe these revisions significantly improve the clarity and impact of the Discussion. We are grateful\nfor the reviewer’s insightful suggestions, which helped us better articulate the biological significance and\ntranslational implications of our findings.\nReviewer #2:\nRemarks to the Author:\nThe authors conducted an elegant multiomic analysis to explore the connection between inflammatory\nbowel disease (IBD) and genetic predisposition, applying a sophisticated computational workflow that,\nto the best of my knowledge, has not been previously implemented. Their approach is undoubtedly innovative\nand contributes valuable insights.\nHowever, the manuscript's structure and flow present challenges, making it difficult for a broad readership\nto follow. The extensive use of technical terminology and abbreviations complicates the overall clarity,\npotentially obscuring the manuscript's key message.\nDespite the novelty of the computational approach and the compelling nature of the study, the work lacks\nmechanistic validation of its findings in IBD, which appears to be the primary aim as suggested in the title.\nGiven the complexity of IBD pathogenesis, attributing T-cell inhibition solely to PGE2 seems overly\nsimplistic. A broader perspective incorporating additional mechanisms underlying T-cell-mediated\nimmunity in IBD would enhance the work’s depth.\nI strongly recommend that the authors refine and streamline the text to improve accessibility while\nemphasizing the relevance of their findings to IBD pathogenesis. Expanding the discussion beyond genetic\nassociations to highlight key inflammatory and immune-mediated pathways—fundamental aspects of\nIBD—would significantly strengthen the manuscript's clinical impact.\nR: We sincerely thank the reviewer for their thoughtful and constructive feedback. In response to these\nvaluable suggestions, we have made the following revisions to improve the clarity and impact of our\nmanuscript:\n1. We have revised the Introduction and Discussion sections to better highlight the key biological\nquestions and major discoveries of the study. Specifically, we have restructured the Discussion to\nhighlight our central findings: the activation state-specific effects of PGE2, its differential\nreprogramming of Th17, Th1, and Treg cells, and the mechanistic insights into PTGER4 regulation.\nWe have also incorporated a conceptual model describing how PGE2 may intersect with genetic\nrisk variants and regulatory elements to modulate T cell phenotypes relevant to IBD. In doing so,\nwe have shifted technical descriptions of the DOGMA-seq assay and computational pipelines to\nthe Methods section, thereby improving accessibility for a broader readership and reducing the\ntechnical density of the main text.\n2. To further enhance clarity, we have moved several technical descriptions related to the activation\nscore to the Methods section. In addition, we have added a table of abbreviations in the\nSupplementary Table 18 to assist readers in following the manuscript.\n3. In response to the reviewer’s concern about oversimplification of attributing T cell inhibition\nsolely to PGE2, our study did not aim to fully recapitulate the complex mucosal microenvironment\nof IBD. Rather, we used a reductionist approach to study the effects of PGE2, TGF-β, or PGE2\nand TGF-β in T cells that were also exposed to activation stimuli, IL-1β, and IL-23 based on\nrationales described in the Introduction and the Discussion. We intentionally studied T cells from\nhealthy, young adult donors to avoid confounding variables such as effects of IBD-directed\nmedications, comorbidities, and aging on the immune responses that we aimed to study, and to\nalign with the typical age of IBD onset in adults. Within this well-controlled experimental system,\nthe observed inhibition and reprogramming effects were indeed primarily attributable to PGE2.\nWe now more explicitly communicate the motivation for our experimental design and\nacknowledge its limitations in the Discussion. Given the current lack of high-resolution single-cell\nmulti-omic studies specifically investigating the effects of PGE2 in human T cells that are exposed\nto activation stimuli, IL-23 and IL-1β, which play key roles in Th17 biology and IBD pathogenesis,\nwe believe that our dataset lays important groundwork for future patient-based studies.\n4. Finally, we have expanded the Discussion to better establish relevance of our findings within\nthe broader framework of IBD pathogenesis and immune-mediated inflammatory pathways. In\nparticular, we emphasize how PGE2 signaling interacts with T cell activation states and genetic\narchitecture, including loci like PTGER4, to modulate transcriptional and epigenetic programs.\nWe also highlight the translational relevance of these insights, especially the potential of the EP4\nreceptor for PGE2, encoded by PTGER4, as both a genetically regulated effector and druggable\ntarget. In addition, we point to the broader utility of our single-cell multi-omic T cell atlas as a\nresource for exploring shared effector circuits across multiple immune-mediated disorders. While\nwe acknowledge that functional validation remains an important next step, our study establishes a\nmechanistically grounded framework and a set of actionable hypotheses for future investigations,\nranging from CRISPR perturbation studies to patient-derived T cells and in vivo disease models.\nWe are grateful for the reviewer’s comments, which have helped us improve the clarity, focus, and clinical\nrelevance of our manuscript. We hope that the revised version more effectively communicates both the\nconceptual novelty and the broader translational potential of our findings.\nReviewer #3:\nRemarks to the Author:\nAs a reviewer coming from a clinical/transational background, I am not in a position to fully evaluate the\ntechnical aspects of the DOGMA-seq single-cell multiomics platform. However, I find the\ncharacterization and depth of analysis regarding the number and subsets of T cells quite impressive and\ncommendable.\nThat said, I believe the claim that this constitutes a model of IBD inflammation—based solely on in vitro\nstimulation of T cells with IL-1β, IL-23, and either PGE2, TGF-β, or both—appears overly simplistic.\nWhile the approach clearly identifies key response elements to these specif-ic immune mediators, it does\nnot recapitulate the complex immunopathology of IBD as seen in patients.\nIn line with this, I would have appreciated data derived from the DOGMA-seq platform on mucosal T\ncells obtained from IBD patients, both during remission and active disease, or from other immune-\nmediated inflammatory diseases (IMIDs), to strengthen translational relevance.\nFinally, I would encourage the authors to reflect on whether their DOGMA-seq platform data provide any\ninsights into the potential relationship between PGE2 signaling and non-response to anti-TNF therapy,\nan important and clinically relevant area that may benefit from the mechanistic resolution their platform\nenables.\nR: We sincerely thank the reviewer for their thoughtful and constructive comments, and for recognizing\nthe depth and quality of our T cell characterization. Below, we provide point-by-point responses.\n1. We agree with the reviewer that our in vitro model cannot fully recapitulate the complex\nmucosal microenvironment of IBD. As clarified in the revised Introduction and Discussion, our\ngoal was not to simulate the full immunopathology of IBD, but rather to model key aspects of\nIBD-relevant immune signaling, focusing on the effects of PGE2, TGF-β, or PGE2 and TGF-β in\nT cells that were also exposed to activation stimuli, IL-1β, and IL-23 based on rationales described\nin the Introduction and the Discussion.\nImportantly, we intentionally used T cells from healthy, young adult donors for several reasons: a)\nTo avoid confounding effects from immunomodulatory or biologic therapies often used in IBD\nand other immune-mediated diseases, which could influence T cell responses in unpredictable\nways; b) To exclude potential confounding from immune aging, ensuring more consistent T cell\nbiology across donors; c) The age range of our study subjects aligns with the peak onset age for\nIBD in adults, improving the relevance of our findings.\nMoreover, because there are currently no published single-cell multi-omic datasets focused on\nPGE2-mediated changes in T cells from healthy human donors that are also exposed to activation\nstimuli, IL-23 and IL-1β, which play key roles in Th17 biology and IBD pathogenesis, our study\nprovides a necessary and well-controlled foundation for future work on patient samples. While\npatient-derived mucosal T cells would indeed further enhance translational relevance, such work\nis beyond the scope of the current study. However, we fully recognize its value and plan to pursue\nthis in future research.\n2. We thank the reviewer for the insightful suggestion to consider potential links between PGE2\nsignaling and non-response to anti-TNF therapy2, a clinically important and timely direction. In\nour dataset, we observed that TNF expression was significantly downregulated across nearly all\nactivated T cell subtypes (with the exception of CD4+ regulatory and Gamma Delta cells) upon\nPGE2 stimulation. This may reflect a feedback mechanism by which PGE2 suppresses TNF-\nmediated effector responses, aligning with its known anti-inflammatory roles in certain contexts.\nFurthermore, recent single-cell RNA-seq studies in IBD patients have shown that non-remission\nanti-TNF ulcerative colitis patients exhibit elevated expression of IL21, TNFRSF1B, and immune\ncheckpoint genes (LAG3, CTLA4, TNFRSF4, TNFRSF18, HAVCR2) in Th17 cells3. Using this\ngene set, we computed an anti-TNF non-remission module score and found that PGE2-stimulated\nTh17 cells consistently showed significantly higher scores in both resting and activated states\n(Figure 1A–B; Wilcoxon rank-sum test: P = 0.00014 for resting, P < 2.2×10⁻¹⁶ for activated).\nThese findings suggest that PGE2 exposure may be associated with poor anti-TNF responsiveness\nin Th17 cells.\nFigure 1. PGE2 promotes T cell states associated with anti-TNF non-response in Th17 cells.\n(A–B) Anti-TNF non-remission module scores (based on IL21, TNFRSF1B, LAG3, CTLA4,\nTNFRSF4, TNFRSF18, and HAVCR2) in resting (A) and activated (B) CD4+ memory Th17 cells\nunder IL-1β/IL-23 stimulation with or without PGE2. Module scores were significantly elevated\nin PGE2-stimulated cells in both states (Wilcoxon rank-sum test, P = 0.00014 for resting, P <\n2.2×10⁻¹⁶ for activated), indicating a shift toward T cell phenotypes observed in anti-TNF non-\nresponders.\n(C) Drug2cell analysis showing reduced predicted sensitivity to the anti-TNF agent Adalimumab\n(CHEMBL1201580) in PGE2-stimulated Th17 cells. Scores reflect enrichment of drug target gene\nexpression mapped to ADALIMUMAB.\n(D) Distribution of Th17 cells in WNN UMAP space under IL-1β/IL-23 stimulation with (orange)\nor without (blue) PGE2. Reduced Adalimumab target scores in (C) correspond to the PGE2-\nstimulated cluster in (D), supporting a potential link between PGE2 exposure and anti-TNF\nresistance phenotypes.\nTo further explore this hypothesis, we applied Drug2cell4 to score enrichment of anti-TNF target\ngenes in PGE2-stimulated vs. unstimulated Th17 cells (Figure 1C–D). Using target gene sets for\nADALIMUMAB derived from the ChEMBL database5, we found that PGE2-stimulated Th17\ncells exhibited significantly lower drug–target scores, implying reduced predicted sensitivity to\nTNF blockade (P = 3.27×10⁻¹⁵, FDR = 2.82×10⁻¹⁴, log₂FC = –1.19). Notably, prior studies have\nimplicated COX-2–PGE2 signaling in impaired epithelial regeneration and reduced\nresponsiveness to TNF inhibitors in UC6, consistent with our observations. These findings\nreinforce the concept that PGE2 may contribute to treatment resistance in IBD by reprogramming\nT cells toward non-responder states. While PGE2 can exert anti-inflammatory effects, its\nactivation of regulatory pathways in Th17 cells, already implicated in anti-TNF non-response—\nraises important considerations for therapeutic strategies.\nGiven the scope and mechanistic focus of the current study, which centers on mapping PGE2-\nmediated chromatin and transcriptional changes, we have opted not to include these findings in\nthe main text. However, we believe they offer a valuable translational perspective and plan to\nexplore them further in future patient-focused studies that directly integrate clinical treatment\noutcomes.\nWe sincerely appreciate the reviewer’s thoughtful feedback, which has helped us strengthen both the\nclarity and translational relevance of our manuscript.\nReferences:\n1. Lee J, Aoki T, Thumkeo D, Siriwach R, Yao C, Narumiya S. T cell-intrinsic prostaglandin E(2)-\nEP2/EP4 signaling is critical in pathogenic T(H)17 cell-driven inflammation. J Allergy Clin\nImmunol 143, 631-643 (2019).\n2. Genaro LM, et al. Anti-TNF therapy and immunogenicity in inflammatory bowel diseases: a\ntranslational approach. Am J Transl Res 13, 13916-13930 (2021).\n3. Thomas T, et al. A longitudinal single-cell atlas of anti-tumour necrosis factor treatment in\ninflammatory bowel disease. Nat Immunol 25, 2152-2165 (2024).\n4. Kanemaru K, et al. Spatially resolved multiomics of human cardiac niches. Nature 619, 801-810\n(2023).\n5. Santos R, et al. A comprehensive map of molecular drug targets. Nat Rev Drug Discov 16, 19-34\n(2017).\n6. Li Y, et al. COX-2-PGE(2) Signaling Impairs Intestinal Epithelial Regeneration and Associates\nwith TNF Inhibitor Responsiveness in Ulcerative Colitis. EBioMedicine 36, 497-507 (2018)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Thank you for addressing all comments so completely.\nTo a large extend the manuscript now reads a lot better.\nThere is duplication in information now, which should be removed tomorrow improve readability.\nAll in all the edits have made the goal - method - conclusion of the research much more clear.\nThe strength in the method used and championed by this manuscript lies in the richness of the information on a single cell level: having communicating layers of RNA, open chromatin and cell surface markers is powerful in itself.\nThe main issues of the conclusions drawn here are\n- based on a limited set of individuals (25) who are healthy (might not be a bad thing, stimulation is controlled which is not the case in diseased individuals, however, we are not well informed whether these stimuli equal disease stimuli)\n- specific set of stimuli on which main conclusions are based are only a subset of the total initial study\n- conclusion is quite narrow as if this were a hypothesis-p[roving study instead if a hypothesis generating study (which I think it actually is)", "reviewer_2": "The authors clearly respond to all my concerns.", "reviewer_4": "The reviewer notes that the authors have addressed the primary concerns regarding the simplification of their in vitro model and the use of healthy donor T cells. The response appropriately acknowledges that the system does not recapitulate the complex mucosal microenvironment of IBD, and the justification for using healthy donors to avoid confounding from medications, comorbidities, and immune aging is reasonable. The clarifications make it clear that the study is intended as a mechanistic, hypothesis-generating analysis, rather than a definitive model of IBD.\nRegarding the potential link between PGE2 signaling and anti-TNF non-response, the authors have provided additional analyses showing that PGE2-stimulated Th17 cells display higher anti-TNF non-remission module scores and reduced predicted sensitivity to Adalimumab using Drug2cell. These analyses, while exploratory and not included in the main text, provide a mechanistic rationale connecting their findings to a clinically relevant phenotype.\nHowever, the manuscript continues to overstate the applicability to IBD, particularly in framing the stimulation conditions as an IBD microenvironment. The data derive from healthy T cells exposed to IL-1β and IL-23 with or without PGE2 (and TGF-\nβ), which allows for identification of PGE2-specific effects within this context. Nevertheless, the manuscript does not discuss the effects of IL-1β/IL-23 stimulation alone compared to unstimulated controls, nor the combined effect of IL-1β/IL-23/PGE2 versus baseline, leaving a gap in understanding how much of the transcriptional and chromatin changes are attributable to cytokine activation versus PGE2-specific modulation.\nThe study’s strongest disease-relevant contribution comes from the integration of IBD GWAS SNPs, which serve as a “test case” to prioritize genes and putative regulatory elements. Emphasizing the creation of a large, multimodal T cell resource and its utility in systematically linking genetic risk variants to regulatory programs would better reflect the primary translational contribution, rather than over-interpreting the results as modeling the IBD mucosal microenvironment.\nOverall, the authors have satisfactorily addressed the initial concerns regarding the scope and translational interpretation of their work. The manuscript would benefit from reframing to temper disease-specific claims, discussing the baseline effects of cytokine stimulation relative to unstimulated controls (if performed - and if not; a discussion thereof), and highlighting the\nDOGMA-seq resource and GWAS integration as the key contributions."}, "author_response": "POINT-BY-POINT RESPONSE TO REVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nThank you for addressing all comments so completely.\nTo a large extend the manuscript now reads a lot better.\nThere is duplication in information now, which should be removed tomorrow improve readability.\nAll in all the edits have made the goal - method - conclusion of the research much more clear.\nThe strength in the method used and championed by this manuscript lies in the richness of the\ninformation on a single cell level: having communicating layers of RNA, open chromatin and cell\nsurface markers is powerful in itself.\nR: We thank the reviewer for this positive and constructive feedback. We have carefully revised\nthe manuscript to improve readability and removed redundant statements, particularly those\nrepeatedly emphasizing the advantages of the multi-omic approach. These edits reduce duplication\nwhile preserving the core scientific messages. We appreciate the reviewer’s suggestion, which\nhelped us further clarify the structure and narrative of the manuscript.\nThe main issues of the conclusions drawn here are\n- based on a limited set of individuals (25) who are healthy (might not be a bad thing, stimulation\nis controlled which is not the case in diseased individuals, however, we are not well informed\nwhether these stimuli equal disease stimuli)\nR: We thank the reviewer for this thoughtful comment. We agree that our study is based on a\nlimited number of healthy participants with ages 18-35 years (n = 25), which is an important\nconsideration. We deliberately chose healthy individuals to enable tightly controlled stimulation\nconditions and to avoid confounding variables such as age-related immune changes and\nmedications, disease heterogeneity, and chronic inflammation that are inherent to patient samples.\nImportantly, the stimuli used in our experiments (T cell activation, IL-23, IL-1β, and PGE2, and\ninitially, TGF-β) are well established to play central roles in IBD and other autoimmune diseases.\nThis design allowed us to systematically interrogate the effects of defined inflammatory signals in\na reproducible manner.\nAt the same time, we acknowledge that in vitro stimulation may not fully recapitulate the\ncomplexity of the in vivo disease microenvironment. We have now added this point as a limitation\nin the Discussion and emphasized that future studies using patient-derived samples and in vivo\nmodels will be essential to validate and extend our findings.\n- specific set of stimuli on which main conclusions are based are only a subset of the total initial\nstudy\nR: We thank the reviewer for this observation. In our initial experimental design, T cells were\nexposed to activation stimuli, IL-1β and IL-23 alone or in combination with PGE2, TGF-β, or both\nPGE2 and TGF-β. However, as reported in the Results (Supplementary Fig. 5), we observed\nminimal effects of TGF-β on T cell states, cell-type proportions, and regulatory programs in\nsamples from the first seven study participants. Therefore, subsequent analyses focused on PGE2\nin a data-driven manner. We have clarified this rationale more explicitly in both the Results and\nMethods sections to justify why downstream analyses centered on PGE2 stimulation.\n- conclusion is quite narrow as if this were a hypothesis-p[roving study instead if a hypothesis\ngenerating study (which I think it actually is)\nR: We fully agree with the reviewer. We have revised the Conclusion to better reflect the\nhypothesis-generating nature of our study. Specifically, we softened causal language (e.g.,\nchanging “revealed” to “suggested”) and reframed key findings to emphasize nomination rather\nthan confirmation of PTGER4 regulatory mechanisms. In addition, we explicitly highlighted the\nneed for future functional validation, including CRISPR-based perturbation, analyses of patient-\nderived samples, and in vivo models, to test the hypotheses generated by our integrative analyses.\nWe appreciate the reviewer’s suggestion, which helped us improve the accuracy and framing of\nour manuscript.\nReviewer #4 (Remarks to the Author):\nThe reviewer notes that the authors have addressed the primary concerns regarding the\nsimplification of their in vitro model and the use of healthy donor T cells. The response\nappropriately acknowledges that the system does not recapitulate the complex\nmucosal microenvironment of IBD, and the justification for using healthy donors to avoid\nconfounding from medications, comorbidities, and immune aging is reasonable. The clarifications\nmake it clear that the study is intended as a mechanistic, hypothesis-generating analysis, rather\nthan a definitive model of IBD.\nRegarding the potential link between PGE2 signaling and anti-TNF non-response, the authors have\nprovided additional analyses showing that PGE2-stimulated Th17 cells display higher anti-TNF\nnon-remission module scores and reduced predicted sensitivity to Adalimumab using Drug2cell.\nThese analyses, while exploratory and not included in the main text, provide a mechanistic\nrationale connecting their findings to a clinically relevant phenotype.\nHowever, the manuscript continues to overstate the applicability to IBD, particularly in framing\nthe stimulation conditions as an IBD microenvironment. The data derive from healthy T cells\nexposed to IL-1β and IL-23 with or without PGE2 (and TGF-β), which allows for identification of\nPGE2-specific effects within this context. Nevertheless, the manuscript does not discuss the effects\nof IL-1β/IL-23 stimulation alone compared to unstimulated controls, nor the combined effect of\nIL-1β/IL-23/PGE2 versus baseline, leaving a gap in understanding how much of the transcriptional\nand chromatin changes are attributable to cytokine activation versus PGE2-specific modulation.\nThe study’s strongest disease-relevant contribution comes from the integration of IBD GWAS\nSNPs, which serve as a “test case” to prioritize genes and putative regulatory elements.\nEmphasizing the creation of a large, multimodal T cell resource and its utility in systematically\nlinking genetic risk variants to regulatory programs would better reflect the primary translational\ncontribution, rather than over-interpreting the results as modeling the IBD mucosal\nmicroenvironment.\nOverall, the authors have satisfactorily addressed the initial concerns regarding the scope and\ntranslational interpretation of their work. The manuscript would benefit from reframing to temper\ndisease-specific claims, discussing the baseline effects of cytokine stimulation relative to\nunstimulated controls (if performed - and if not; a discussion thereof), and highlighting the\nDOGMA-seq resource and GWAS integration as the key contributions.\nR: We thank the reviewer for this thoughtful and constructive evaluation of our revised manuscript.\n1. Regarding the concern that we overstated the applicability of our in vitro system to IBD,\nwe agree and have revised the manuscript throughout to temper disease-specific language.\nSpecifically, we now more precisely describe our experimental system as healthy donor T\ncells exposed to activation stimuli, IL-1β and IL-23 alone or in combination with PGE2,\nTGF-β, or both PGE2 and TGF-β, rather than framing it as a direct surrogate for the IBD\nmucosal microenvironment. This clarification emphasizes that our design was intended to\ninterrogate PGE2-specific effects (and initially, also TGF-β-specific effects) within a\ndefined inflammatory signaling context, rather than to fully model in vivo disease\nconditions.\nWe have now explicitly added a limitations paragraph in the Discussion, noting that our\nstudy was conducted in a relatively modest number of healthy participants with ages 18-\n35 years and relied on in vitro stimulation to model inflammatory conditions. We clarify\nthat this design enabled tightly controlled experimental perturbations and avoided\nconfounding variables such as age-related immune changes, and medications, disease\nheterogeneity, and comorbidities inherent to patient samples, while also acknowledging\nthat in vitro stimulation cannot fully recapitulate the complexity and dynamics of the in\nvivo disease microenvironment. We further emphasize that future studies incorporating\npatient-derived samples, tissue-resident immune cells, and in vivo models will be essential\nto validate and extend our findings in physiologically and clinically relevant contexts.\n2. Regarding the reviewer’s comments about effects of cytokine stimulation relative to\nunstimulated controls, we agree that this is an important consideration. Given the well-\nestablished roles of IL-23 signaling and the Th17 immune pathway in IBD pathogenesis\nand of IL-1β and IL-23 in driving human T cell differentiation, all T cells that we studied\nwere exposed to activation stimuli, IL-1β and IL-23. We specifically interrogated PGE2-\nmediated modulation (and initially, also TGF-β-mediated modulation) in this context, and\nwe did not include unstimulated controls. We now explicitly acknowledge in the\nDiscussion that we cannot fully disentangle transcriptional and chromatin changes driven\nby activation stimuli, IL-1β and/or IL-23 alone from those specifically attributable to PGE2\nalone, and that future studies incorporating such baselines will be important to disentangle\ncytokine-driven effects from PGE2-specific regulatory programs.\n3. Regarding the reviewer’s comment that the strongest disease-relevant contribution of this\nstudy lies in the integration of IBD GWAS data, we agree. In response, we have revised\nboth the Introduction and Discussion to more explicitly emphasize this aspect.\nSpecifically, we added the following statements to the Introduction: “Leveraging this\nmultimodal omics dataset, we systematically characterized T cell heterogeneity, quantified\nPGE2-associated transcriptional and epigenetic changes, and integrated these signals with\nIBD GWAS variants to prioritize candidate genes and putative cis-regulatory elements.\nUsing IBD genetics as a test case, we demonstrated how this resource can be applied to\nlink disease-associated non-coding variants to cell type-specific regulatory programs.”\nThis highlights the use of IBD genetics as a test case to demonstrate how our resource can\nbe leveraged to link non-coding risk variants to cell type–specific regulatory programs.\nIn addition, we revised the Discussion to emphasize that “a major translational contribution\nof this study is the generation of a single-cell multi-omics T cell resource and its integration\nwith IBD GWAS data to systematically prioritize regulatory elements and candidate genes\nacross cell states.”\nConsistent with this reframing, we have shifted the narrative emphasis away from modeling\nthe IBD mucosal microenvironment and toward the broader utility of the DOGMA-seq\nresource and GWAS integration framework. As noted above, we also tempered disease-\nspecific claims and expanded the Discussion of baseline controls and study limitations, in\nline with the reviewer’s suggestions.\nWe appreciate the reviewer’s guidance, which helped us clarify the primary translational\ncontribution and strengthen the overall framing of the manuscript."}]

Decision letter