{'fraction_non_alphanumeric': 0.029459901800327332, 'fraction_numerical': 0.02648236152465837, 'mean_word_length': 5.229455840805798, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 3, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A single-cell transcriptional gradient in human cutaneous memory T cells restricts Th17/Tc17 identity

Summary The homeostatic mechanisms that fail to restrain chronic tissue inflammation in diseases, such as psoriasis vulgaris, remain incompletely understood. We profiled transcriptomes and epitopes of single psoriatic and normal skin-resident T cells, revealing a gradated transcriptional program of coordinately regulated inflammation-suppressive genes. This program, which is sharply suppressed in lesional skin, strikingly restricts Th17/Tc17 cytokine and other inflammatory mediators on the single-cell level. CRISPR-based deactivation of two core components of this inflammation-suppressive program, ZFP36L2 and ZFP36, replicates the interleukin-17A (IL-17A), granulocyte macrophage-colony-stimulating factor (GM-CSF), and interferon gamma (IFNγ) elevation in psoriatic memory T cells deficient in these transcripts, functionally validating their influence. Combinatoric expression analysis indicates the suppression of specific inflammatory mediators by individual program members. Finally, we find that therapeutic IL-23 blockade reduces Th17/Tc17 cell frequency in lesional skin but fails to normalize this inflammatory-suppressive program, suggesting how treated lesions may be primed for recurrence after withdrawal of treatment.

In brief Cook et al. show that Th17/Tc17 expression occurs in a gradient within psoriatic single T cells, paralleled by attenuation of an inflammation-suppressive geneexpression program termed ZIST. Persistent ZIST suppression, even in therapeutically improved psoriatic lesions, suggests it contributes to recurrence of disease after treatment withdrawal.

INTRODUCTION
Cytokine signaling drives common types of pathologic skin inflammation, such as the Th1 and Th17 populations causative for psoriasis vulgaris. 1 Over the past 10 years, targeted inhibition of these pathways has generated effective, biologic therapies for diseases such as psoriasis and revealed skin-resident T cell populations that support recurrent disease. 2 However, many mechanistic aspects of tissue inflammation remain unclear, particularly how pathogenic activation of these pathways is restrained between clinical flares. Nor is it certain how initial inflammatory activation amplifies within immune cell populations to produce symptomatic skin lesions.
The emergence of single-cell RNA sequencing (scRNA-seq) has made it possible to analyze cutaneous biological responses in ever higher resolution. 3-6 These profiles can uniquely reveal inflammatory abnormalities in immune and stromal cells, providing a finegrained portrait of how pathogenic changes reprogram complex tissues. However, these approaches have not yet closely examined how healthy, uninflamed T cells transition to a pathogenic state in chronic skin disease. In experimental settings outside of skin, T cell populations appear comprised of continuums of cellular identity, which robustly influence inflammatory activation. 7 It is yet unknown if similar gradients operate in clinical disease and how they control pathogenic immunological behavior.
To investigate T cell identity and consequence in an in vivo setting, we scRNA-seq-profiled CD45 + immune cells isolated from lesional skin in eight patients with active cases of the prototypical Th17/Tc17 skin disease psoriasis vulgaris and compared them with seven uninflamed, normal controls. In addition to searching for molecular abnormalities distinguishing inflamed from normal cells 3 , we sought to determine if such alterations correlated with each other in biologically specified programs. We also asked how pathogenic immune cell identities in psoriasis lesions, as revealed by scRNA-seq, might evolve during therapeutic treatment of interleukin-23 (IL-23) pathway overactivity.

RESULTS
scRNA-seq of T cells isolated from psoriatic lesions and uninflamed skin reveals both pathogenic Th17/Tc17 activation and suppression of inflammation-suppressive genes scRNA-seq using the 10x Genomics Chromium droplet-based methodology was performed in tandem with cellular indexing of transcriptomes and epitopes by sequencing (CITE-seq; see STAR Methods) on CD45 + immune cells isolated from lesional skin biopsies of 8 patients with psoriasis, comparable in cohort size to recent scRNA-seq studies (Tables S1 and S2). 4,8 Cutaneous CD45 + cells were also isolated and profiled from seven uninflamed donors. We performed high-resolution unsupervised clustering on CD3 + cells and examined major classes that each represented more than 5% of all T cells in our samples (Table S3). We defined these six T cell classes passing this filter based on their transcriptional and epitope (CITE-seq) markers ( Figure 1A), including central memory cells (Tcms), two tissue-resident memory (Trm1 and Trm2) populations, one exhausted and one activated cytotoxic CD8 + T lymphocyte population (CTLex and CTLac, respectively), and regulatory T cells (Tregs). Tcms were CD45RA À (epitope)/ CD62L + /CCR7 + , while the two closely related Trm classes were CD69 + (epitope)/CD103 +9, 10 and also CXCR6 +11 , distinguished from each other by only a small set of non-canonical markers (Table S3). Relative to CTLac, CTLex harbored elevated PDCD1, LAG3, and PRF1. All six classes were  Table S3; Figure S1). We next examined differential gene expression in T cell subpopulations in the 8 psoriasis samples against our 7 uninflamed, normal controls using a pseudo-bulk approach. 12,13 Elevation in psoriatic T cells of Th17/Tc17 cytokines such as IL17A, IL17F, IL26, 14 IFNG, and CXCL13 15 occurred predominantly in a population harboring skin-resident memory T cells markers (Trm1; Figure 1C; Table S4). IL22 expression was lower and thus less reliable for bulk comparison in these data. A broad array of T cell activation markers such as CTLA4, CD2, and CD82 were also elevated in these single cells expressing psoriatic cytokines, as was suppression of a series of inflammation-suppressive transcripts not previously reported as a major feature of T cells in psoriasis or other inflammatory skin diseases. Such transcripts include ZFP36L2 and ZFP36, two tristetraprolin family members previously shown to be capable of repressing inflammatory potential in murine and human T and B cells. [16][17][18][19] In fact, among the 140 transcripts repressed >1.5 average log2 fold change (FC) in the Th17/Tc17 cytokine-expressing Trm1 cell class with an adjusted p value <0.01, we also observed multiple other genes known to suppress inflammatory cellular identity, including inhibitors of JAK/STAT (SOCS1 20 and SOCS3 21 ) and Th17 signaling (CD69 22 ), negative regulators of CD8 + cytotoxic function (LGALS1 23 ), and promoters of Treg activity (CCL22 24 ). Downregulation of these inflammation-suppressive genes was also observed in psoriatic samples in other T cell classes, most strongly Tcm and Trm2 (Table S4), although T cells expressing canonical psoriatic cytokines classified primarily in the Trm1 category.
Because Th17 and Tc17 cells have each been reported as arising from skin-resident memory cells, 2,25,26 we were curious if both populations could be identified in our skin-resident memory T cells, which harbored virtually all expression of psoriatic cytokines. Psoriatic T cells in the Trm1 subpopulation were subclassified as CD4 + and CD8 + based on expression of CD4, CD8A, or CD8B (STAR Methods; Table S3) as in other recent studies. 4 As seen in Figure 1D, both the psoriatic CD4 + and CD8 + subpopulations substantively expressed IL17A, IL17F, IL26, and IFNG; we hereafter refer to these classes as Th17 and Tc17, respectively.
A ZFP36L2-centric cell identity gradient defines pathogenic cytokine and inflammatory transcript expression in psoriatic skin-resident memory T cells on the single-cell level The detection of these transcriptional differences allowed us to evaluate two differing hypotheses regarding T cell identity in psoriasis. First, upregulation of psoriatic inflammatory transcripts and downregulation of inflammation-suppressive individual gene transcripts might simply occur in an unrelated pattern within Th17/Tc17 single cells. Alternatively, expression of these two gene sets might be coordinately linked in specific, related patterns on the single-cell level. To distinguish these possibilities, we examined correlations between IL17A and IL17F and other highly expressed transcripts in skin-resident memory T cells. To eliminate artifacts, gene values were batch corrected at the sample level using the CPCA method in the R package iCellR (STAR Methods). Intriguingly, not only was expression of IL17A and IL17F each highly correlated with one another, and also with IFNG, but the greatest levels of anti-correlation were observed against many of the inflammation-suppressive transcripts we described above as downregulated in this class. In fact, by this measure, ZFP36L2 was the single most anti-correlated transcript out of 16,348 assessed genes versus both IL17A (d = À0.58, p < 2 3 10 À5 ) and IL17F (d = À0.65, p < 2 3 10 À5 ) expression in Th17/Tc17 cells (Figures 2A and S2). A similar, high degree of anti-correlation was also observed between psoriatic cytokine transcripts and ZFP36, whose expression closely tracks that of ZFP36L2 in single T cells (Table S5).
The correlation graphs in Figure 2B show that ZFP36L2 and inflammatory transcript expression opposed each other in single Th17/Tc17 cells as a continuous gradient rather than, for example, in discrete classes. This gradient was relatively muted in uninflamed Trms and most evident in the transition from Trms in healthy, uninflamed tissue to Th17/Tc17 cells in lesional psoriatic skin. This phenomenon presented visually as the blue, normal T cells observed at the upper left end of each gradient (with high ZFP36L2 expression and low cytokine expression), in contrast to the dark red lesional T cells extending toward the lower right. Substantial expression of key Th17 cytokines such as IL17A, IL17F, and IL26 was restricted sharply to psoriatic single Th17/Tc17 cells whose ZFP36L2 level was suppressed below that of Trms in normal, uninflamed skin.
It was not clear whether such gradient patterns arose primarily out of either the Th17 or Tc17 subpopulations. We therefore repeated this analysis after compartmentalizing CD4 + and CD8 + subpopulations as described above. As illustrated in Figure S3, the gradients involving ZFP36L2 were robustly detected in both subpopulations, confirming a recurrent patterning of pathogenic cytokines within both Th17 and Tc17 identities. Thus ZFP36L2 deficiency, paralleled by depressed ZFP36 levels, appears to specify a pathogenic inflammatory state in psoriatic T cells, a transcriptional identity largely absent in the corresponding resident-memory T cells isolated from non-inflamed skin ( Figure 2B, blue data points).
The prominent anti-correlation between tristetraprolin family members and Th17/Tc17 cytokines also led us to evaluate their relationship beyond human skin T cells. In mice, thymically derived murine gd T cells produce IL17 isoforms and drive psoriasiform inflammation. These gd T cells express cognate Trm markers, proliferate in situ in response to local tissue cytokines, and do not freely recirculate out of the skin and are thus frequently regarded as skin-resident analogs in psoriasis models. In re-analyzing single-cell data from imiquimod-provoked gd T cells in mice, 31 we again observed a steep anti-correlation between levels of Zfp36l2 and the Th17/Tc17 cytokines Il17a, Il17f, and Il22 ( Figure 2C). In data from stimulated human peripheral T cells, 7 we also found that IL17F was expressed almost exclusively in low ZFP36L2-expressing single cells ( Figure S4), indicating that this relationship extends to psoriatic cytokine expression outside of tissue-resident settings.
While we primarily detected and modeled significant psoriatic cytokine expression in cells falling in the Trm1 subpopulation, anti-correlation of other inflammatory markers against ZFP36L2/ZFP36 levels was also observed in our Tcm, Treg, and CD8 + classes ( Figure 2D; Table S5), suggesting that tristetraprolin family member expression pervasively defines inflammatory cell identity across skin T cell classes.  Figure S2. (B) ZFP36L2 suppression predicts induction of numerous other pro-inflammatory mediator transcripts from Th17/Tc17 in lesional psoriatic skin from this study (red, all 8 samples pooled) versus healthy controls (blue, all 7 samples pooled). The y axis shows ZFP36L2 expression, the x axis shows imputed cytokine transcript levels, and each point represents a single T cell. (C) scRNA-seq from prior studies of CD45 + cutaneous immune cells isolated from imiquimod-provoked (red, 3 samples pooled) versus control treated mice (blue, 3 samples pooled) 31 show maximal expression of IL17A, IL17F, and IL22 with highest ZFP36L2 suppression. The y axis shows ZFP36L2 expression, the x axis shows imputed cytokine transcript levels, and each point represents a single gd T cell. For both axes, the standard imputed expression value (iCellR) has been normalized and log-transformed (Seurat; STAR Methods). (D) Relationship between ZFP36L2 loss and inflammatory mediator expression seen in all major skin T cell classes but strongest in skin-resident memory T cells. ZFP36L2 or ZFP36 knockout in primary CD4 + and CD8 + T cells elevates protein and transcript expression of 3 0 ARE-containing cytokines We were curious if marked suppression of ZFP36L2 and ZFP36 mRNA abundance simply correlated with increased expression of key cytokines or, instead, functionally elevated their levels. Both ZFP36L2 and ZFP36 bind AU-rich elements (AREs) in the 3 0 untranslated region (UTR) of mRNA transcripts, directly inhibiting ribosomal access and translation. Tristetraprolin family members have been shown to suppress ARE-harboring cytokine transcript levels in mouse T cell populations 32 and in human Tregs. 33 We utilized CRISPR-Cas9 ribonucleoprotein complexes to genetically disrupt both loci in primary peripheral T cells because of the technical challenges in obtaining human cutaneous T cells in sufficient numbers for these assays. As the relationship between ZFP36L2 and Th17/Tc17 cytokine production was broadly conserved in both tissue-resident and other skin T cell classes ( Figures 2B and 2D), we targeted knockout to bulk CD4 + and CD8 + T cells. Two distinct guide RNAs (gRNAs) were designed for both ZFP36L2 and ZFP36; greater than 81% allele knockout was confirmed for each gRNA, using both Sanger-based (Synthego ICE 34 ) and next-generation sequencing approaches (knockout percentages in Table S6; exemplar Sanger traces shown in Figure S5). A negative control gRNA for the human AAVS1 safe-harbor integration locus was tested to control for non-specific effects of Cas9-mediated double-strand breaks.
Using biological triplicates of our CRISPR-Cas9 ribonucleoprotein-based knockouts of ZFP36L2 and ZFP36 ( Figure S5), we assessed levels of 3 0 ARE-containing cytokines by flow cytometry. Tumor necrosis factor alpha (TNF-a) and granulocyte macrophage-colony-stimulating factor (GM-CSF) are known to be elevated in psoriatic T cells, 35 although usually below detection limits for 10x Genomics 3 0 scRNA-seq analysis. 4 As shown in Figure 3A, knockout of ZFP36L2 in both CD4 + and CD8 + T cells significantly induced interferon gamma (IFNg), TNF-a, and IL-17A. By comparison, ZFP36 knockout induced IFNg and TNFa in both populations but did not appreciably affect IL-17A expression ( Figure 3B). GM-CSF was assayed specifically for the ZFP36 knockout, based on a previously reported regulatory relationship. 36 The observed cytokine induction under these conditions was substantial. ZFP36L2 knockout resulted in as much as a 1.6fold increase in IFNg-positive CD8 + T cells, from 29% to 47% (guide L2g2, p < 0.01; Figure 3A, second row). ZFP36L2 knockout also induced TNF-a 2.4-fold in CD8 + T cells, from 4.2% to 10% positivity (guide L2g2, p < 0.01), and IL-17A 2.1fold (p < 0.0001; Figures 3A and S5). Assessments for TNF-a and IFNg were repeated in CRISPR knockout CD4 + and CD8 + T cells obtained from a different donor for one guide ( Figure S5).
As the binding of tristetraprolin family members to cytokine mRNA also targets them for deadenylation and degradation, 37 we also assessed ZFP36L2 influence on transcript levels. We bulk RNA-seq-profiled CD4 + primary human T cells after ZFP36L2 disruption performed with three distinct donors (Table S7; STAR Methods). Applying the differential expression analysis package sleuth to these data, 38 among transcripts for secreted cytokines, we detected that ZFP36L2 knockout significantly increased levels of IL17A and IFNG, which are elevated in ZFP36L2-deficient T cells (p < 0.01; Figure 2A). However, other pro-inflammatory genes elevated in ZFP36L2-deficient T cells but lacking 3 0 AREs, such as CXCL13 and CD82, were not induced by its disruption, suggesting that there are additional mechanisms responsible for their upregulation. To further distinguish direct targets of tristetraprolin proteins in T cells, we examined data from a recent study employing a binding-based biochemical assay to identify Zfp36 binding partners in wildtype CD4 + cells. 19 Of 382 transcripts showing high-magnitude suppression versus ZFP36 level in our study (r < À0.5, p < 0.01), 199 were also identified in this screen (Table S5).

ZFP36L2 defines programmatic attenuation of other inflammation-suppressive transcripts in rash lesional T cells
Our findings support a model in which repressed ZFP36L2 and ZFP36 expression, in psoriatic Th17/Tc17 cells, elevates the abundance of 3 0 ARE-containing cytokines. However, many other pro-inflammatory transcripts that were elevated in ZFP36L2/ZFP36-deficient T cells lack such 3 0 regulatory elements (e.g., CD2, CD82, and T cell receptor components), suggesting a different regulatory cause. We speculated that other inflammation-suppressive genes might be downregulated in concert with ZFP36L2 in single Th17/Tc17 cells, more broadly activating an inflammatory state. We thus examined the 50 transcripts most correlated with ZFP36L2 in Trms/Th17/Tc17 cells (r > 0.4, <0.2% of assessed transcripts, p < 6.1 3 10 À364 , excluding ribosomal or mitochondrial transcripts; Table S5) for prior biological evidence of immune pathway restriction. Strikingly, many inflammation-suppressing genes that we noted earlier as downregulated in psoriatic T cells were found in this set, including ZFP36, SOCS1, SOCS3, CCL22, and CD69, as well as the global quiescence enforcers BTG1 39 and TXNIP 40 ( Figure 4A). We thus termed these genes, which appear coherently suppressed in the highest cytokine-expressing single Th17/Tc17 cells, the ZFP36L2 inflammation-suppressive transcriptional program, or ZIST. ZFP36L2 correlates with other ZIST genes in both skin-resident CD4 + and CD8 + subfractions ( Figure S3). However, its statistical correlation is weaker in the Treg subpopulation and CD8 + cells outside of the Trm1 cluster ( Figure 4B). ZIST transcripts were not identified as significantly altered in T cells in a study of collagenase artifacts in human tissue processing. 41 We hypothesized that the elevated expression of multiple inflammation-suppressive genes (i.e., ZIST) normally restricts T cell inflammatory tone. Exogenous stimulation then downregulates ZIST components, an effect observed for both ZFP36L2 and ZFP36 in peripheral T cells, 7 potentiating inflammation. While expression of ZIST members appears strongly correlated within Trms/Th17/Tc17 cells, we were curious if their individualistic effects on specific inflammatory targets might still be distinguishable on a single-cell level. Such associations would implicate individual ZIST members as regulators of distinct aspects of psoriatic T cell activation. To test this hypothesis, we implemented a single-cell, combinatoric approach to test correlations between (1) transcript levels of all possible, non-repeating sets of predictor (i.e., ZIST) genes and (2) expression of individual inflammatory mediator genes described above.
By testing all possible ZIST sets in their prediction of inflammatory transcript expression, this approach is capable of testing correlative effects of groups of genes but also quantitatively discerning if components of sets have outsize effects on targets. Briefly, the batch-corrected matrices described above were passed to Monocle3, and all possible combinations of predictor (ZIST) genes were then summed using the aggregate_gene function, in order to calculate Spearman rank correlation with respect to inflammatory mediator target genes (Table S8; STAR Methods).
These analyses reveal two intriguing features of psoriatic inflammatory activation in resident T cell populations. First, while downregulated ZFP36L2 and ZFP36 rank as the best individual predictors of inflammatory gene expression, combinations of suppressed ZIST program genes consistently perform even more strongly ( Figure 4C; Table S5). Therefore, we hypothesize that a yet unidentified global regulator coordinately suppresses these ZIST program genes, producing a single, major inflammatory gradient in T cells from psoriatic skin. Secondly, even within this dominant axis, our combinatorial analysis can still distinguish gene-specific patterns. For example, ZFP36L2 suppression is more predictive than that of ZFP36 for IL17A (r = À0.513 versus À0.416) and IL17F (r = À0.625 versus À0.486) expression. In contrast, combined ZFP36L2 and ZFP36 suppression better predict expression of IFNG (r = À0.481 versus À0.418) and CXCL13 (À0.578 versus À0.490) than ZFP36L2 alone.
Notably, these findings reflect our CRISPR-based experimental data that ZFP36L2, but not ZFP36, represses IL-17A in human CD4 + and CD8 + T cells, while both suppress IFNg (Figure 3). Figure 4D also shows an influence of BTG1 suppression on CD82 expression not observed for cytokine targets; this prediction agrees with the experimental induction of CD82 by BTG1 knockout, 39 suggesting that our model may predict specific regulatory relationships for additional validation. While these analyses are intended for examination of broad trends rather than ascribing significance to each individual ZIST combination, we also implemented a stepwise model comparison algorithm using R/cocor, establishing significance between combinatoric pairs (STAR Methods; Table S8).
ZIST program expression defines a trajectory reflecting an inherent molecular gradient from normal to psoriatic skin-resident memory cells Given the coherence of the ZIST program in Th17/Tc17 single cells, we thought that its magnitude might reveal an inherent trajectory defining molecular evolution from normal to inflamed T cells in psoriasis. To examine such a natural trajectory, also referred to as a pseudotime, we optimized a manifold consisting of Trms/Th17/Tc17 cells from all 15 samples and overlaid expression of key transcript sets. Briefly, we calculated an integrated expression value for each gene program within each cell, representing the summed standardized values of each constituent gene within a set. The cells were then organized in twodimensional space based on a rooted trajectory.
After batch correction, this trajectory shows evenly distributed contributions from all our samples ( Figure 5A), and unsurprisingly, its rooted, unsupervised gradient is defined by a transition from normal Trms to psoriatic Th17/Tc17 cells ( Figure 5B). This trajectory is not strongly biased by either CD4/CD8 expression ( Figure 5C), nor is it a product of central memory markers such as CCR7 or SELL ( Figure 5D). As we speculated, the major transition from normal to pathologic T cells is strongly anti-correlated with ZIST program magnitude ( Figure 5E) and is correlated with cytokine and inflammatory mediator expression ( Figures 5F and  5G), a phenomenon that is robustly represented on the individual-sample level ( Figure S6). The rooted trajectory is notable for its linear nature, with few outlying developmental sinks (Figure 5H). This analysis thus reinforces a model in which psoriatic skin-resident T cells harbor an inflammatory identity inversely related to expression of ZIST program members.
Therapeutic IL-23 blockade drastically reduces Th17/ Tc17 cell frequency, but residual Trms in healed psoriasis lesions show persistent ZIST program suppression We next sought to understand how treatment of psoriasis using IL-23 blockade would affect our inflammation-suppressive gradient. Although IL17A and IL17F cytokine expression should be largely quenched by IL-23 blockade, the pro-inflammatory identity represented by ZIST suppression might persist in residual T cells. To assess this hypothesis, we isolated and performed scRNA-seq on lesional T cells from three patients with psoriasis prior to initiation on the IL-23-blocker tildrakizumab. After at least 8 weeks of biologic treatment, which produced a 75% or greater reduction in psoriasis area and severity index score in all three patients, we re-biopsied and profiled clinically resolved skin adjacent to the original sample sites for comparison. Tildrakizumab dramatically reduces Th17/Tc17 cells in psoriatic lesions, an expected feature of blocking IL-23 signaling to skin-resident memory cells, both in absolute number ( Figure 6A) and, importantly, as a relative percentage of T cell class (Table S9).
However, in residual resident memory T cells of the Trm1 class isolated from clinically resolved psoriasis lesions, our scRNAseq analysis also shows persistent abnormal suppression of ZIST elements, with mid-treatment versus normal comparison Kruskal-Wallis p values <2 3 10 À16 for all transcripts ( Figure 6B; patient-level data in Table S9). Plotting ZIST program member expression against ZFP36L2 levels in Trms/Th17/Tc17 cells from these pre-and mid-treatment psoriasis lesions shows that despite resolution of visible inflammatory plaques, the overall inflammatory-suppressive gradient in T cells remains sharply suppressed ( Figure 6C). Some, but not all, inflammatory transcripts lacking AREs also appear persistently elevated in the context of clinical improvement, including CD2 and CD3G (Table S9). These data thus reveal that skin-resident memory cells under IL-23 blockade continue to harbor a molecularly primed state closely resembling that in inflamed skin. Such

DISCUSSION
Th17/Tc17 cytokines are well-established to play a central, therapeutically actionable role in the development of psoriasis, but these pathways only partly explain the clinical course of skin inflammatory disease. Dynamic, secondary mechanisms likely activate immune cell populations during cutaneous flares, which are rapid and unpredictable. Our analysis of single-cell sequencing data from psoriasis lesions reveals a reactive layer of heterogeneity embedded in cutaneous T cells. Specifically, a transcriptional single-cell identity involving multiple inflammation-suppressive regulators (described here as the ZIST program) appears to play a dominant role in defining the most pathogenically activated psoriatic T cells. ZIST program members like ZFP36L2 and ZFP36 post-transcriptionally suppress pre-formed Th17/Tc17 cytokine transcripts and are known to be rapidly de-activated, 42,43 nominating a mechanism by which chronic rashes might rapidly flare.
Some ZIST members outside of ZFP36L2 and ZFP36 are wellcharacterized global regulators of T cell activation, such as BTG1, or SOCS proteins, which appear to substantially inhibit JAK/STAT signaling. 44 Others, such as CCL22, may lower T cell inflammatory tone both by direct action on Tregs 24 and also dendritic cell populations. 45 Like the numerous genes that are upregulated in and augment the inflammatory capacity of inflamed immune cells, the exact delineation of such a program is partially arbitrary and likely to be substantially tissue specific.
Our findings reveal that ZIST inflammation-suppressive transcripts exist in a single-cell continuum. This gradient is then coordinately shifted downward in psoriatic lesions, relieving suppression of inflammation and driving cytokine expression in the most ZIST-depressed subpopulation. In fact, we find that IL17F and IL26 expression emerges only in those psoriatic T cells whose ZIST expression is depressed below the range observed in normal control samples ( Figure 2B), underscoring the tight link between this gradient and pathogenic inflammation. We investigate this regulatory pattern primarily in our Trm1 subpopulation, in which Th17/Tc17 cytokines are expressed in psoriasis samples and can be studied directly, but other inflammatory markers negatively correlate with ZFP36L2-deficiency in other T cell classes ( Figure 2D). A recent study also observed ZFP36L2 suppression in stimulated Tregs, as seen in our data, apparently elevating Treg tone and perhaps representing a feedback mechanism that is overwhelmed in the context of pathogenic inflammation. 33 We propose that in psoriatic lesion development, general, exogenous inflammatory insults suppress elements of the ZIST program in healthy Trms, leading to enhanced pathogenic psoriatic cytokine and inflammatory mediator transcript levels. An increasing inflammatory milieu then further depresses ZIST levels, as previously shown with in vitro suppression of ZFP36L2 and ZFP36 in peripheral T cells by T cell receptor (TCR) stimulation. 7,33 Our pre-and mid-tildrakizumab treatment data support this model, in which IL-23 signaling and ZIST regulation represent dual, distinguishable inputs, as IL-23 blockade sharply reduces Th17/Tc17 cell density but does not normalize ZIST suppression in Trms from healed psoriatic lesions. Although our patients represent a limited initial cohort that must be expanded, such a molecularly primed, inflammatory cell identity persisting despite IL-23 blockade might accelerate recurrence of psoriatic lesions after withdrawal of treatment.
ZFP36L2 and ZFP36 attenuation help explain cytokine induction in ZIST-deficient T cells but are unlikely to directly elevate the numerous other pro-inflammatory transcripts lacking 3 0 ARE sequences, their binding target. Such inflammatory markers include CD2, CD3G, and CD82, which are known to establish the T cell immunological synapse and maintain signaling downstream of the TCR, 46 and their upregulation is likely to further activate more greatly ZIST-deficient single T cells. Potential sources for regulation of these other inflammatory markers include ZIST members such as BTG1, whose suppression tracks specifically with CD82 elevation ( Figure 4D) and which is not affected by ZFP36L2/ZFP36 knockdown (Table S7), a relationship that has been mechanistically validated in murine T cells. 39 Therefore, the combinatoric model we present here may be used to predict and further experimentally validate other suppressor-inflammatory mediator relationships.
The persistence of ZIST gradient abnormalities in psoriasis despite IL-23 blockade nominates this regulatory mechanism as an intriguing, complementary therapeutic target in overactivity of the Th17/Tc17 axis, as does recent work indicating that fibroblasts may also suppress psoriatic inflammation through ZFP36L2. 47 It also remains to be understood whether this transcriptional gradient acts similarly on Th2 cytokines in atopic dermatitis or in other chronic skin inflammatory disorders. One key technical limitation of our study is the inability to enhance activity of ZFP36L2 or other ZIST elements in human skin. The tightly linked expression of ZIST program elements strongly suggests that a yet undiscovered global regulator produces this gradient cell identity, analogous to factors such as FOXO1 and KLF2 that guide the transition from naive to memory cells. A means to selectively activate ZFP36L2 and other ZIST genes in human T cells remains an important milestone for future work.  Table S6 for CD4 + /CD8 + ratios).

Limitations of the study
This study is limited by the relatively small number of samples analyzed, given the current costs and technical challenges of scRNA-seq. It is also not possible to formally exclude influences of tissue dissociation on the skin-resident T cell transcriptional profiles we report here. Finally, we lack the technical capability to directly genetically modify T cells isolated from human skin, limiting direct, in vitro validation of the inflammation suppression function of our candidates. These technical bounds in genetic manipulation also restrict our capacity to measure the aggregate impact of downregulation of multiple inflammation-suppressive transcripts in our T cell gradients.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

INCLUSION AND DIVERSITY
We worked to ensure gender balance in the recruitment of human subjects. We worked to ensure ethnic or other types of diversity in the recruitment of human subjects. One or more of the authors of this paper self-identifies as an underrepresented ethnic minority in science. One or more of the authors of this paper self-identifies as a member of the LGBTQ+ community. One or more of the authors of this paper received support from a program designed to increase minority representation in science. While citing references scientifically relevant for this work, we also actively worked to promote gender balance in our reference list. The author list of this paper includes contributors from the location where the research was conducted who participated in the data collection, design, analysis, and/or interpretation of the work. Data and code availability d Sequence data is submitted at the European Genome-Phenome Archive (EGA), which is hosted by the EBI and the CRG, under accession number EGA: S00001005271. Flow cytometry data reported in this paper will be shared by the lead contact upon request. d All analysis scripts are available at the online repository https://github.com/cpcook1/TTP. All statistical analysis and plotting of scRNA-seq and cell surface protein data were performed using Rstudio software (v1.2.5033). d Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Patient characteristics
Normal, healthy control skin obtained from patients undergoing mastectomies or reduction abdominoplasties and psoriasis lesional skin were obtained using protocols approved by the UCSF Human Research Protection Program Institutional Review Board. Informed consent was obtained from all patients donating samples based on these active protocols. Adult patients donating psoriasis samples from the trunk or proximal extremity were assessed at a Psoriasis Area Severity Index (PASI score) of 8 or more (moderate to severe disease). Histopathology was verified by a board-certified dermatopathologist. Patients undergoing IL23 blockade (tildrakizumab) treatment were enrolled under an IRB protocol as part of a study funded by Sun Pharma, had a pre-treatment PASI score of 7 or greater, and received at least two 100 mg treatments before mid-treatment biopsy of clinically resolved lesions at > 8 weeks. Comprehensive biographical metrics of all subjects from this study are available in Table S1.

METHOD DETAILS
Skin sample collection, single cell RNA-seq, and CITE-seq processing Skin samples were obtained by 6 mm punch biopsy and minced into fine pieces using surgical scissors and transferred into 3 ml of RPMI-1640 medium supplemented with 10% fetal bovine serum, 100 IU/mL Penicillin+100 mg/mL Streptomycin, 10mM HEPES, collagenase type IV (200U/ml; Worthington) and 20 mg/ml DNAse I (MilliporeSigma), then incubated for 16-18h at 37 C with 5% CO 2 . After incubation, the suspension was filtered through a 100 mm cell strainer, pelleted by centrifugation at 1500 rpm for 5 min, and resuspended in 100 mL of FACS buffer (PBS supplemented with 2% FBS) containing APC conjugated human CD45 Antibody at a 1:33 dilution. After staining at 4 C for 30 min, cells were washed twice in FACS buffer, passed through a 45 mm filter and sorted for live CD45 + cells following the addition of DAPI at 1 mg/ml. For samples where CITE-seq was performed (Table S2), murine splenocytes were added to comprise 5% of total cells (to allow for computational adjustment for non-specific antibody binding). Blocking was performed with the addition of 5% Human TruStain FcX (Biolegend) in a 100 mL cell suspension volume for 10 minutes at 4 C, followed by CITE-seq antibody addition for 30 min at 4 C (TotalSeq-A antibodies, 0.5 ug for each antibody per 1 million cells; Biolegend: Table S2). Cells were washed three times in Cell Stain Buffer (BioLegend) and concentrated by pelleting and resuspended in a smaller volume prior to Chromium chip loading. Single-cell RNA-sequencing libraries were then prepared from isolated cells using the Chromium Single cell 3' Solution V2 or V3 kit (10x Genomics) following the manufacturer's protocol by the Genomics Core Facility UCSF Institute for Human Genetics (Table S2). For CITE-seq samples, 0.2 pmol of ADT (antibody derived tag) additive primer was added at the RNA library cDNA amplification step. CITE-seq libraries were prepared according to the manufacturer's TotalSeq-A antibody protocol (BioLegend). In brief, 70 mL of ADT-containing cDNA amplification supernatant was purified with two rounds of 2X SPRI beads (Beckman-Coulter) and then amplified for 14-20 cycles using HiFi HotStart ReadyMix (KAPA, Roche Sequencing & Life Science) and 0.25 mM of oligos corresponding to SI PCR primer and Truseq Small RNA RPI1-6 primers. The resulting amplification products were purified by 1.2X SPRI bead cleanup, then quantified with Qubit dsDNA HS Assay Kit. Quality for scRNA and ADT libraries was assessed by a TapeStation D1000 ScreenTape (Agilent Technologies Inc.) and quantitated by the Kapa library quantitation kit prior to sequencing. mRNA and ADT libraries were sequenced on a HiSeq 4000 (Illumina) with a read length configuration of 150 PE.
Single cell RNA sequencing data processing FASTQ files derived from sequencing were aligned and quantified using Cell Ranger Software (v3.0.2, 10x Genomics) against the human GRCh38 reference transcriptome. A filtered data matrix was used which excludes empty droplets and doublets. Quality of cells were then assessed based on the total number of detected genes per cell and the percentage of mitochondrial gene counts. Cells were filtered if there were unique feature counts greater than 5000 or mitochondrial gene counts higher than 20%. Normalized counts were then computed based on feature counts for each cell divided by the total counts for that cell and multiplied by the scale factor (1e4) with natural-log transformation. Finally, about 500 to 6,000 cells per sample were retained for downstream analysis (Table S2). Based on these criteria, high-quality transcriptomic data was obtained from 68,835 single cells after removal of doublets and cells with high mitochondrial content or low number of genes. The median gene count per cell was 1351 postfiltering and the average total number of detected genes per donor was 18,048 (detailed sequencing metrics for each sample are in Table S2).

Dimensional reduction and unsupervised clustering
The harmony algorithm 48 was used for batch effect correction and integration of cells into a coordinated space for unsupervised clustering. To detect high variable features (HVFs) used for harmony, the variance stabilizing transformation parameter was used as a function input. After scaling the data, a PCA matrix was calculated using detected HVFs with 20 components based on the ElbowPlot function in Seurat v4.0.2. We then provided this PCA matrix to the RunHarmony() function in Seurat using samples as technical covariates for correction. The batch-corrected coordinated space then was used to compute the nearest neighbor graph by the FindNeighbors() function.
Indeed, across cell types prior to batch correction, samples clearly fell out as coherent clusters in the UMAP ( Figure S1D), indicating that a large proportion of transcriptome-wide variation correlated to sample. However, upon batch correction, sample-specific cells no longer cohered but disease-level clusters continued to be clear ( Figure S1D). We therefore concluded that sample-level batch correction successfully eliminated potential sample-level confounding while preserving target disease-level variation.
This nearest neighbor graph was used by FindClusters() using the Louvain algorithm for clustering of cell populations and dimensional reduction such as UMAP and tSNE. Based on clustree optimization, we utilized resolution = 0.4 to obtain 14 initial immune clusters (6 CD3 + or KLRB1 + lymphocyte and 8 HLA-DRA + antigen presenting cell clusters). We grouped the lymphocyte cell clusters independently and further subclustered them into 21 CD3 + or KLRB1 + clusters to achieve finer resolution. To filter out cells containing markers previously identified as associated with collagenase digestion, and thus representing potential artifacts 41 , we removed clusters 15 and 21 (as numbered in the Seurat object), which were highly enriched for JUN, FOS, HSPA1B, and HSP1B. We retained for analysis the six major T cell populations that each represented more than 5% of all T cells from our sample dataset. Cluster-specific differentially expressed genes were detected using the FindMarkers() function in the Seurat package using the MAST method and manually annotated based on the expression of known marker genes and then were checked using the singleR package v3.13. 49 The Human Primary Cell Atlas was used as a reference dataset for annotation.
To identify genes differentially expressed between psoriasis and matched normal samples, we developed pseudo-bulk counts for each cluster in each patient using the Bioconductor package muscat 12 , which then applies edgeR to detect statistically significant differences between the patient groups. 13 To deconvolute CD4 + and CD8 + T cells in the Trm1 population, we took a conservative transcriptional gating approach to ensure that we called high-confidence helper and cytotoxic T-cells. We used empiric unimputed data, filtering out non-expressing or double positive cells (double positive were defined as cells expressing CD4/CD8A and/or CD4/CD8B). We then defined CD4 + cells as those expressing detectable CD4 transcripts and CD8 + cells as those expressing either CD8A or CD8B transcripts. For analysis in Figure 6, all cells containing at least one IL17A + or IL17F + read were counted and then corrected for sampling of the overall. The correction factor was (total CD45 + cells isolated in the biopsy/total CD45 + cells profiled by scRNA-seq).

scRNA-seq correlation analysis
To identify genes associated with ZFP36L2 expression, we used Spearman rank correlations to search through cluster-specific sets of high variability genes. First, we identified the top 5000 variable genes across cells within sample-and cluster-specific subpopulations using the Seurat function FindVariableFeatures. Within these cell subpopulations, we then tested for the association between ZFP36L2 and each of these highly variable genes with Spearman correlations implemented in the cor.test function from the Stats package v3.6.2. We accounted for multiple testing by adjusting p-values with the Benjamini-Hochberg procedure.
For Figures 2A, 2B, and 4A, the standard imputed expression value shown on the x-axis represents count data that has been normalized and log-transformed by R/Seurat's default standardization and normalization algorithms, and imputed by R/iCellR/ run.impute(my.obj, dims = 1:10, nn = 10, data.type = pca).

Trm1 manifold and developmental trajectory inference
We analyzed Trm1 subpopulation cells with the goal of understanding how unbiased developmental trajectories could be inferred from them, using these spatially separate subjects as a space-for-time substitutions (i.e. how cells could develop from uninflamed to inflamed states using samples separated in space among different subject to infer temporal transitions in pseudotime). Since human subjects differ in a number of potentially confounding ways, we controlled for sample-level variation by implementing a mutual nearest neighbor batch correction routine to 'regress out' the effect of sample-specific variation, presuming that the residual variation would be produced by disease-level rather than sample-level variation. 59 Indeed, this proved to be the case since disease correlated strongly to axes 1 and 2 of the resulting uMAP, indicating that the plurality of variation in the data after accounting for sample-specific variation correlated to disease state. We then undertook a manifold optimization routine to balance local and global structure in the data based on: 1) global connectedness to support the inference of a global trajectory and 2) local clustering to isolate specific disease and cell-state subpopulations. First, we re-composed Trm1 cells in a uMAP, exploring 3,000 combinations in Nn (number of nearest neighbors [5,50] and Md (minimum nearest neighbor distance [0.1,0.9], finding 10 and 0.8 respectively to strike the best balance between global and local information on manual inspection. We then used Leiden-based clustering at default resolution (k=10, partition q-value=0.05) to call unsupervised clusters for downstream lineage inference. 50 Since the true branching structure and hierarchy in Trm1 cells are unknown, we used Monocle3 to infer a rooted trajectory since it has been benchmarked in controlled trials as the most informative and robust method to interrogate complex developmental topologies in single cells. 51 We then used the Monocle3/learn_graph function without partition-consideration to infer a global trajectory. 52, 53 We then e4 Cell Reports Medicine 3, 100715, August 16, 2022 Article ll OPEN ACCESS overlaid clinical information on the uninflamed or inflamed provenance of each cell onto the manifold, knowing that the true trajectory begins with uninflamed cells. Node 1 occurred in a highly dense cluster of uninflamed cells oriented along axes 1 and 2 away from the majority of inflamed cells, and thus, we rooted our global trajectory in this node and calculated pseudotemporal ordering from it using the order_cells function in Monocle3. 54 We then observed that pseudotemporal order and ZPF36L2 expression were broadly correlated to Axis 1 and testing for the Spearman rank correlation between them revealed a significant anticorrelation. Since this result concords with inflammation-associated downregulation in ZPF36L2 observed in clinical specimens, we considered this topology to reflect true ZFP36L2-associated dynamics in these cells. We then repeated these associations by integrating over specific sets of genes using the aggregate_gene function in Monocle3 without maximum or minimum cutoffs. These scaled, integrated gene expression values enabled the direct comparison of the total expression of sets of genes among cells with a single value, which were overlaid on the manifold and correlated with pseudotime as in Figure 5. The association between the integrated expression of these programs and pseudotime were not linear and thus we fit LOESS curves to visualize overall trends.
Combinatorial single-cell inflammatory mediator prediction analysis Since the aggregation of gene sets smoothed over gene-level variation within these sets, it was not clear which constituent genes within a set may be driving correlational patterns observed in Figure 6. In order to test this directly, we implemented a combinatoric approach to test the correlation of all possible, non-repeating combinations of genes of all set sizes between predictive sets and single response genes. The predictive set consisted of ZFP36L2, ZFP36, CCL22, SOCS1, SOCS3, BTG1, TXNIP, CD69; and the response genes were IL17A, IL17F, IL26, IFNG, CXCL13, CD2, CD3E, CD3G, CD82, TFRC, LAYN, GNLY, ENTPD1, CPM, NPDC1, ZEB2. Gene values were batch-corrected at the sample level using the CPCA method in the R package iCellR; missing gene values were independently imputed within inflamed and unflamed states of sample-aligned matrices using the PCA method in iCellR/run.impute. Resulting matrices were then passed to Monocle3, and all possible combinations of predictor genes were then summed using the aggregate_gene function, and their correlation with response genes calculated with Spearman rank correlations using the cor.test function.
To develop pairwise statistical tests of the predictive power of these ZIST gene combinations, we developed a stepwise model comparison algorithm using R/cocor. We used this package to implement a two-sided test between the correlations of independent variables (in this case the gene sets). This method relies on Hittner, May, and Silver's modification of Dunn and Clark's Z statistic 55 , which in turn uses a backtransformed Fisher's Z procedure. To control type 1 error by shrinking the test space, we ordered pairwise steps according to correlation coefficients of gene sets and tested for correlation differences only between two adjacent gene sets.
Donor T cell nucleofection T cells were isolated from donor leukopaks (STEMCELL Technologies) by negative selection with CD3 magnetic beads (EasySep Human T Cell Isolation Kit; STEMCELL Technologies) according to the manufacturer's protocol and frozen in LN 2 at 20x10 6 cells/ mL for storage. Donor T cells were then thawed, washed twice in X-VIVO 15 media (Lonza) then resuspended at 1x10 6 cells/mL in T cell medium (X-VIVO 15 supplemented with 5% FBS, 55 mM 2-mercaptoethanol and 200 IU/mL recombinant human IL2) prior to the addition of CD3/CD28 conjugated magnetic activation beads (Dynabeads; Invitrogen) at a 1:1 bead-to-cell ratio. Forty-eight hours later the cells were washed 1x in PBS, resuspended in 20 mL of supplemented P3 nucleofection solution (Lonza) at 50x10 6 cells/ mL, and immediately added to the RNPs and mixed by gentle pipetting. The cell/RNP mixture was placed in one well of a 16-well cuvette strip (P3 Primary Cell 4D-Nucleofector X Kit S; Lonza) and electroporated using program EH115 on a Lonza X-Unit nucleofector. Immediately after, 100 mL of warm X-VIVO+5% FBS was added to the well, and the cells were placed in a 5% CO 2 incubator at 37 C for 15 min, then plated in 4 mL of T cell medium. Cells were monitored daily and split 1:1 when the density reached 2x10 6 cells/ mL or media showed signs of acidification.