Human regulatory T cells locally differentiate and are functionally heterogeneous within the inflamed arthritic joint

Abstract Objective Tregs are crucial for immune regulation, and environment‐driven adaptation of effector (e)Tregs is essential for local functioning. However, the extent of human Treg heterogeneity in inflammatory settings is unclear. Methods We combined single‐cell RNA‐ and TCR‐sequencing on Tregs derived from three to six patients with juvenile idiopathic arthritis (JIA) to investigate the functional heterogeneity of human synovial fluid (SF)‐derived Tregs from inflamed joints. Confirmation and suppressive function of the identified Treg clusters was assessed by flow cytometry. Results Four Treg clusters were identified; incoming, activated eTregs with either a dominant suppressive or cytotoxic profile, and GPR56+CD161+CXCL13+ Tregs. Pseudotime analysis showed differentiation towards either classical eTreg profiles or GPR56+CD161+CXCL13+ Tregs supported by TCR data. Despite its most differentiated phenotype, GPR56+CD161+CXCL13+ Tregs were shown to be suppressive. Furthermore, BATF was identified as an overarching eTreg regulator, with the novel Treg‐associated regulon BHLHE40 driving differentiation towards GPR56+CD161+CXCL13+ Tregs, and JAZF1 towards classical eTregs. Conclusion Our study reveals a heterogeneous population of Tregs at the site of inflammation in JIA. SF Treg differentiate to a classical eTreg profile with a more dominant suppressive or cytotoxic profile that share a similar TCR repertoire, or towards GPR56+CD161+CXCL13+ Tregs with a more distinct TCR repertoire. Genes characterising GPR56+CD161+CXCL13+ Tregs were also mirrored in other T‐cell subsets in both the tumor and the autoimmune setting. Finally, the identified key regulators driving SF Treg adaptation may be interesting targets for autoimmunity or tumor interventions.


INTRODUCTION
Regulatory T cells (Tregs) comprise a subset of CD4 + T cells crucial in preserving immune homeostasis by antagonising immune responses. The transcription factor (TF) FOXP3 characterises Tregs, and mutations in the FOXP3 gene lead to severe inflammation in both mice and humans. 1,2 In recent years, potential therapeutic strategies targeting Tregs in both the autoimmune and the tumor setting have been explored. In autoimmunity, the number and/or functionality of Tregs should be enforced, whereas in the tumor milieu, the suppressive capacity of Tregs should be dampened. 3,4 This can, amongst others, be achieved by expanding Tregs, or by blocking Treg functioning via immune checkpoint blockade. PD-1 and CTLA-4 blockades are employed in the cancer setting, and several other targets are currently tested in clinical trials. 4 Understanding the heterogeneity and plasticity of Tregs is essential in understanding immunodynamics at play in health and disease. This knowledge can then be exploited to develop and improve (potential) Treg-based therapeutic strategies.
Tregs are not identical in every tissue of residence but are tailored to the environment in which they have to function, and can alter their phenotype according to micro-environmental changes over time. 5,6 This plasticity enables continuous adaptation to changing immunodynamics. Upon activation, Tregs gain an effector profile (eTreg) and can initiate transcriptional programmes associated with the dominant T-helper (Th) response at the site of inflammation to enable Treg survival and suppression of the respective Th cells. 5,7,8 These cotranscriptional programmes include Th1 (T-bet), Th2 (GATA-3), Th17 (RORc) and T follicular regulatory cell (Tfr, Bcl6) programmes. 9 We have recently demonstrated that in synovial fluid (SF) of patients with juvenile idiopathic arthritis (JIA), a predominantly Th1/Th17-associated disease, 10 Tregs retain a functional Treg core profile and obtain a Th1-skewed cotranscriptional profile on the proteomic, transcriptomic and epigenetic level. 7 Furthermore, it has been shown that some Tregs can acquire additional functions including stimulation of tissue repair in the intestine 11 and hair follicle growth in the skin. 12 These additional functions further elucidate the local function and relevance of Tregs. An outstanding knowledge gap regarding human Tregs is the degree of heterogeneity of Tregs present within inflammation, if heterogeneity has consequences for Treg function, and how this relates to clonality and thus the T-cell receptor (TCR)repertoire of Tregs.
Tregs in tissues and peripheral blood (PB) have been shown to be heterogeneous, 5,7,8 but data on local inflammatory environments in humans are lacking. Single-cell (sc)RNA-sequencing enables us to map cellular states within an environment to facilitate our understanding of the phenotypic plasticity and functional diversity within the Treg population. SF-derived Tregs from JIA patients provide us with a relevant model of local autoimmune inflammation to determine the functional differentiation and clonality of inflammation-derived Tregs and its regulators. In this study, we aimed at dissecting the heterogeneity of SF-derived Tregs by employing scRNA-sequencing to further elucidate the immunodynamics at play in an inflammatory environment.

Heterogeneity within inflammatory synovial fluid Tregs
To assess the heterogeneity at the site of inflammation in JIA patients, SF Tregs (live CD3 + CD4 + CD127 low CD25 high ) from three patients with oligoarticular JIA were sorted for single-cell transcriptome analysis. Dimensionality reduction of 980 Tregs after quality control revealed the presence of four clusters within the Treg population ( Figure 1a), with each cluster present in all three patients (Figure 1b). All clusters were of Treg origin with the vast majority of cells (97%) expressing at least one FOXP3 mRNA molecule and in 99.96% of the Tregs the human core Treg signature as defined by Ferraro et al. 13 was enriched (Supplementary figure 2b). There was no cluster-specific association with the cell cycle phase the Tregs resided in (Supplementary figure 2c; Pearson's Chi-squared test, P-value = 0.8512).
The largest cluster (37.24%), cluster 1, was characterised by genes that are downregulated upon activation and maturation of T cells. These genes included CCR7, LEF1, KLF2, KLF3 and TCF7, suggesting a relatively quiescent/resting phenotype and probably representing Tregs that only recently migrated into the inflamed joint. The other three Treg clusters all displayed an activated gene signature including expression of many MHC class II genes (e.g. HLA-DR, -DQ, -DP and -DM), but also DUSP4, CTSC, CTSB, ITM2A and LMNA amongst others. Additional markers separated these clusters from each other. Both cluster 2 (31.22%; TIGIT, CTLA4, IKZF2, LAYN) and 3 (23.57%; LGALS1, CXCR6, CCR5, TNFRSF8, GZMA) showed high expression of genes associated with highly suppressive Tregs. The smallest activated Treg cluster, cluster 4 (7.96%), expressed a set of genes not commonly or previously associated with Tregs including CXCL13, GPR56, MYO7A, BHLHE40, PTPN13 and KLRB1 (Figure 1c and d;  Supplementary figures 2d and 3, Supplementary  table 1). Overall, the three activated Treg clusters showed the expression of a wide array of (a) cluster 1 cluster 2   LEF1  CCR7  KLF2  SMCHD1  KLF3  TXNIP  RASA3  TCF7  DGKA  ZFP36L2  SOD2  FAIM3   TIGIT  TBC1D4  SAT1  IKFZ2  CXCR6  GBP2  LAYN  ANTXR2  HLA-DRB1  MAF  PHACTR2  PHTF2   LGALS1  LMNA  CLIC1  APOBEC3C  CCR5  HLA-DRB6  APOBEC3G  HCST  COTL1  CORO1A  LDHA  PHLDA1   CXCL13  MYO7A  GPR56  BHLHE40  RBPJ  PTPN13  KLRB1  FKBP5  LMO4  ID2  PDCD1  IFNG   cluster 4   (c) Top 12 upregulated genes, based on the P-adjusted value, per cluster based on MAST differential gene expression analysis. (d) UMAPs of the combined expression of 2-4 selected differentially expressed genes per cluster shown in nebulosa density (kernel density estimation to handle sparsity of scRNA-sequencing data). The scale ranges from blue to yellow, with in yellow the highest kernel density, thus the highest (estimated) expression of all combined selected genes. (e) Gene set analysis of a gene module downregulated in naive versus memory CD4 T cells (GSE61077). Enrichment of a gene set is calculated per cell; grey signifies no enrichment of the gene set and yellow to red represents increasing enrichment. (f) Similar to e, but for a human shared tissue Treg signature 56 and a human tumor-infiltrating Treg signature. 57 costimulatory and co-inhibitory markers to help suppress immune activation in a mostly contactdependent manner but with a per cluster different potential dominant mode of suppression: for example CTLA4 in cluster 2, GZMA in cluster 3 and LAG3 in cluster 4 (Supplementary table 1).
Microenvironmental cues can shape the transcriptomic signature of Tregs with a resulting cotranscriptional Th cell programme. 5,7 These can be distinguished based on the expression pattern of the chemokine receptors CXCR3, CCR4, CCR6, CCR10 and CXCR5. 15 CXCR5, linked to Tfr, was not expressed in any of the SF Tregs. Cluster 1 harboured Tregs expressing mixed chemokine receptor profiles associated with Th2, Th22, Th1 and Th17 cells. Overall expression of CCR4, CCR6 and CCR10 was low. By contrast, cluster 2-4 Tregs were for 51.7% positive for the Th1-associated CXCR3 (Supplementary figure 4). These data indicate that Tregs in the inflammatory SF environment are heterogeneous, and that local cues preferentially induce a Th1 cotranscriptional programme in these Tregs.
Clonotype sharing amongst synovial fluidderived Treg clusters Next, we assessed whether the TCR of individual Tregs skews differentiation to a certain phenotype upon triggering. Therefore, we employed a 109 Genomics data set, including both 5 0 gene expression and TCR sequences, published by Maschmeyer et al. 16 containing SF-derived Tregs from JIA patients. Dimensionality reduction revealed similar Treg clusters as in our data set indicating that the clustering is robust (Supplementary figure 5a).
On average, 55% of the detected clonotypes (full length combined TRA and TRB chain) were unique within a patient (range 34.8-72.7%), showing the presence of clonal expansion. Exploratory analysis revealed that cluster 1 contained mostly single frequency clonotypes, whereas for cluster 2-4, this was less (average of 65.7% single clonotypes, range of 42.0-89.6% for cluster 1 and average of 49.2%, 44.3% and 61.4% for cluster 2-4, respectively) ( Figure 2a). Clonal expansion was analysed by dividing the detected clonotypes in six groups based on frequency. The proportional space filled by the most expanded clones revealed a gradient from cluster 1 to 4. In cluster 4, the most prevalent clonotypes composed of the greatest proportion of the total clones present ( Figure 2b). A wide spread was observed for cluster 1 for the number of clonotypes comprising 10% of the total repertoire with an average of 16 clonotypes (Figure 2c). For cluster 2 and 3, only five clonotypes, and for cluster 4, merely two different clonotypes were observed on average ( Figure 2c). Several diversity indices indeed indicated that cluster 1 had the most diverse clonal repertoire than the three activated Treg clusters (Supplementary figure 5b).
We also determined the presence of clonal sharing between the clusters, indicative of a shared origin. The overlap coefficient was calculated based on overlap of the complete TRA/ TRB nucleotide sequences. This showed a clonotype overlap of 25% and 32.4% of cluster 1 with clusters 2 and 3, respectively. The latter two clusters shared 32.3% of their clonotypes. However, interestingly, there was < 15% clonotype sharing on nucleotide level between cluster 4 and the other clusters, although on TRA/ TRB chain level, this was 21.1-28.5%. Overall, based on the Morisita similarity index, a measure that takes the total number of cells into account, it was indeed clear there was little similarity between clusters 1-3 and cluster 4 Tregs ( Figure 2d). For the TRA chain alone, the similarity between all four clusters was very high, suggesting that the differences are primarily formed by the TRB (Figure 2e). We also assessed whether any of the SF Treg clusters showed a similar TCR repertoire to SF CD4 non-Treg cells.
For none of the four SF Treg clusters, a similar TCR repertoire to SF CD4 non-Treg cells was observed (Supplementary figure 5c). This suggests that there is minimal in situ differentiation of CD4 T cells to Tregs for any of the SF Treg clusters. Altogether, these data indicate that cluster 2 and 3 eTregs are relatively similar, as also suggested by their transcriptomic profile, whereas cluster 4 Tregs contain more distinct clonotypes primarily skewed by the TRB.

Non-linear differentiation of Tregs within synovial fluid
To further explore how the SF Treg clusters are related, we employed pseudotime analyses to estimate Treg differentiation within SF based on transcriptional similarities. We applied Monocle v3 17,18 to perform trajectory inference and pseudotime plotting (Figure 3a and b). Mathematical assessment of the potential starting node for the differentiation trajectory pointed towards cluster 1 Tregs. Those Tregs were indeed highest in (combined) the expression of the genes CCR7, LEF1, TCF7 and KLF2 associated with a naive state (Figure 3c). The identified trajectory suggests that upon arrival within the SF environment, Tregs are skewed towards a classical eTreg phenotype (clusters 2 and 3) or towards cluster 4 Tregs. Cluster 4 Tregs may, however, also pass a cluster 3 phenotype along the differentiation trajectory (Figure 3a and b).
Transcription factors interact with cis-regulatory elements and regulate cell differentiation, and together with their target genes, it forms a regulon. SCENIC can be employed for gene regulatory network analysis to deduce active regulons per cell. Here, we compared cluster 1 recently migrated Tregs with cluster 2-4 activated (e)Tregs. There was no clear binarisation of the regulons per cluster indicating that there is a gradual differentiation of SF Tregs. The top regulons were all upregulated in cluster 2-4 Tregs and associated with differentiation or maintenance of the Treg phenotype indicating that upon migration to the SF-environment cluster 1 Tregs differentiate towards cluster 2-4 Tregs and not in the opposite way (Figure 3d). BATF is a known key driver in eTreg differentiation, 7,19,20 and indeed, our analysis revealed BATF as the primary local regulon for differentiation of Tregs that recently migrated to SF. RUNX1 and NFYC have also been previously associated with Treg development, maintenance and differentiation, and TBX21 is a possible driver of the Th1-associated cotranscriptional programme, 7,21,22 whereas for example, SP3 is associated with generic processes including proliferation, apoptosis and metabolism. 23 BHLHE40 was identified as a novel regulon for eTreg differentiation, and target genes mostly composed of cluster 4-associated genes including CXCL13, GPR56 and KLRB1. However, also amongst its target genes are genes required for Treg differentiation and survival such as ID2. 24 Other BHLHE40 associated genes, such as RBPJ and RFX1, have been linked to inhibiting Th2 and Th17 differentiation, respectively. 25,26 Another novel identified regulator, that is JAZF1, primarily drives SF Treg differentiation towards cluster 2-3 Tregs. Its target genes include the (e)Treg genes TIGIT and CXCR6, 6,7 but also POU2F2 which regulates a shift towards aerobic glycolysis 27 ( Figure 3e). Overall, these data indicate that Tregs arriving in the inflamed SF-environment can differentiate into activated (e)Tregs via different routes in a non-linear fashion, facilitated by key regulators including BATF, TBX21, RUNX1 and the newly (e)Treg-associated regulons BHLHE40 and JAZF1.
GPR56 + /CD161 + /CXCL13 + synovial fluid Tregs are highly differentiated and suppressive Cluster 4 Tregs were defined by genes not commonly associated with Tregs which warranted further investigations. On protein level, we could confirm the presence of Tregs expressing GPR56 and/or CD161 (average of 18.1%, range 10.3-32.3%). This subset was specifically enriched for CXCL13 expression. In PB Tregs, the expression of GPR56 and/or CD161 averaged 4.6%, with no CXLC13 + Tregs present (Figure 4a). In line with the transcriptomic data, GPR56 + CD161 + CXCL13 + SF Tregs expressed high levels of PD-1 and LAG3 (Figure 4b). We also confirmed FOXP3 protein expression within these Tregs. Although the intensity of FOXP3 expression in GPR56 and/or CD161-positive CXCL13 + SF Tregs was lower than other SF Tregs, it was higher than in PB Tregs and PB and SF non-Tregs (Figure 4c). In addition, a gradual decline of Helios-expressing Tregs was observed from PB Tregs to CXCL13 + SF Tregs. However, compared with non-Tregs, about 5.59 more cells expressed Helios: 53.9% for CXCL13 + SF Treg versus 9.9% for CXCL13 + SF non-Treg (Figure 4d).
Even though Helios cannot distinguish thymic-derived and peripherally derived Tregs, its expression in FOXP3 + Tregs indicates a stable and activated phenotype, 28,29 suggesting that CXCL13 + Tregs maintain fundamental Treg characteristics.
Genes associated with TCR stimulation showed the highest enrichment in clusters 3 and 4 ( Figure 4d). Upregulation of TCR signalling was also reflected in the gene ontology analysis and their previously shown clonality (Supplementary  TIGIT  CD2  CARD16  SRGN  TNFRSF18  CXCR6  CTSC  PKM  IL2RB  OAS1   RBPJ  CXCL13  GPR56  MYO7A  ID2  LMO4  LAG3  TRPS1  TSPAN2  NOTCH2   ID2  PPP3CA  LAT2  ADORA2A  CAPN1  LAG3  LAIR2  APOBEC3D  IRF2BP2  TNIP3   TIGIT  CXCR6  HLA-DPB1  SSR4  S100A4  COMMD3  LMNA  GPX1 ANXA2 POU2F2  CD4 + non-Treg cells can transiently express FOXP3 upon activation, while maintaining effector functions. 30 To verify that our Tregs are bona fide suppressor cells, we assessed Tregspecific characteristics. On transcriptome level, we observed higher expression of BHLHE40, IFNG and ID2 in cluster 4 than in the other SF Tregs, which are TFs associated with cytokine expression. However, on protein level, we could not detect significant levels of IFNc (1.63%), IL-2 (3.88%) or IL-17 (2.2%) produced by GPR56/CD161 + CXCL13 + SF Tregs (Supplementary figure 6a) indicating that these cells are true Tregs. In addition, we performed a suppression assay as developed by Long et al. 31 to assess the suppressive capacity of cluster 4 Tregs. We sorted cluster 4 Tregs based on CCR7 À GPR56 and/or CD161 expression (for the gating strategy, see Supplementary figure 6b), which includes CXCL13 expressing Tregs. By contrast, within the CD161 À GPR56 À Treg population, CXCL13 + Tregs composed of only 0.8% of the cells. In addition, we included CCR7 + Tregs for resting/quiescent Tregs and CCR7 À Tregs negative for both GPR56 and CD161 from both SF and control PB. As expected, we observed that a higher ratio of Tregs to CD4 effector T cells resulted in a higher percentage of CD4 effector T-cell suppression. Furthermore, similar levels of CD4 effector T-cell suppression were observed by all sorted SF Treg subsets (Figure 4f). These data show that Tregs from cluster 4 are proficient suppressor cells without any clear differences amongst Treg clusters present in SF.
setting, but this remains to be fully elucidated. These data support the notion that Tregs within SF follow different and unique routes of adaptation.
In both mice and humans, studies have explored PB to tissue Treg differentiation. However, data on functional adaptation pathways of Tregs in inflamed tissues are lacking. Here, we show that infiltrating Tregs are heterogeneous in chemokine receptor expression upon migration to SF indicating partial unspecific homing in response to inflammatory signals, with local cues and antigen(s) inducing preferential differentiation to, and expansion of, Th1-skewed Tregs. The acquisition of a Th-associated cotranscriptional programme is crucial for Treg survival and function in those inflammatory environments. 37 The demonstrated heterogeneity also suggests that upon local cues Tregs can induce subtle shifts in phenotype, possibly to suppress immune responses under changing inflammatory conditions. Especially for cluster 2 and 3 Tregs, this seems likely since they share many clonotypes but show phenotypic differences in functional Treg markers. For GPR56 + CD161 + CXCL13 + Tregs (cluster 4), however, clonotypes are less overlapping suggesting a more predetermined route of these Tregs, perhaps based on TCRspecificity or affinity, when migrating to the site of inflammation. Thereby, the question rises where the imprinting of Treg differentiation occurs. Chemokine receptor expression is determined in secondary lymphoid structures and functional marker heterogeneity seems to be partially predetermined in the periphery. [38][39][40] Our data suggest that there is further extensive local plasticity and differentiation, which is supported by the observation that different T-cell subsets show a similar phenotypical adaptation in local inflammatory environments.
Several scRNA-sequencing studies in autoimmune and tumor settings have shown the presence of CD4 and CD8 T-cell subsets characterised by CXCL13 and inhibitory immune checkpoints 16,[41][42][43][44] indicating a partial mirrored phenotype amongst inflammation-derived T cells. However, for FOXP3 + Tregs, this phenotype has not been described before. These CXCL13 + subsets express high levels of PD-1 and other inhibitory immune checkpoints, proliferate, are clonally expanded and have lost classical CD4 and CD8 T-cell effector functions including cytokine production. 16,[41][42][43][44] These cells have a T peripheral helper-like phenotype with respect to CXCL13 and PD-1 expression as well as the absence of CXCR5 and Bcl6, 43 but do not produce IL-21. Genes such as PDCD1, CCL5, GNLY, HAVCR2 and CCL4 expressed by these cells suggest a terminal stage of differentiation. However, Li et al. 41 showed that these dysfunctional cells are tumor-reactive, and gaining CXCL13 expression might indicate acquisition of novel functions. 44 TCR stimulation, presence of inflammatory cytokines including type I interferons and especially TGFb are known to induce CXCL13 expression, 45 which are all abundant in inflammatory environments. That these cells are not found in the periphery indicates that CXCL13 expression is acquired locally. CXCL13 expression suggests a role in ectopic lymphoid structure (ELS) formation as CXCL13 is a known B-cell attractant. Although, Tregs can regulate B-cell maturation within lymphoid structures, it is unknown whether they can also help forming ELS via CXCL13 expression. Supportive hereof could be expression of extracellular matrix organisation-related genes, such as ADAM19, COL6A3, COL9A2, CTSH, CCL4 and CCL5 by a fraction of these CXCL13 + Tregs. Expression of these genes could also indicate an involvement in tissue repair via extracellular matrix organisation after inflammation-related damage in JIA patients. 46 Furthermore, in a murine model, tissue repair recently has been shown as a function of Tregs at the site of myocardial infarctions. 47 Additional markers expressed by CXCL13 + cells, such as GPR56 for SF Tregs, seem more cell-specific since they are not co-expressed in CXCL13 + (SF) non-Tregs. Furthermore, for SF GPR56 + CD161 + Tregs (containing CXCL13 + Tregs), there is no loss of classical functioning as these Tregs remain proficient suppressors. This indicates that CXCL13 expression and expression of inhibitory immune checkpoints is likely dependent on the T-cell subset concerned. Whether these CXCL13 + T-cell subsets, and specifically CXCL13 + Tregs, in the autoimmune setting are beneficial or pathogenic needs to be elucidated.
Limitations of this study include the risk of contamination of the sorted Tregs with non-Tregs. However, CD127 high CD25 low Tregs were sorted with a purity > 98% with a FOXP3 intracellular staining confirming > 90% FOXP3 expression and 96.7% of the Tregs contained at least one FOXP3 transcript. Another concern could be the sample size of three JIA patients for the scRNA-sequencing. Therefore, we confirmed our findings in a public 109 Genomics sequenced SF T-cell data set including Tregs. We could also confirm the presence of CXCL13 + Tregs at protein level further supporting our transcriptomic findings. Furthermore, there is still a limited understanding of gene signatures representing functions or states in human Tregs. For example, a definition of exhausted Tregs is lacking with many of the genes defined for CD8 and CD4 T cells actually signifying enhanced Treg activity such as PDCD1, ENTPD1 (CD39) and LAG3. It is of importance to further explore (human) Tregs in inflammatory settings to improve our understanding of their (dys)function.
In recent years, several Treg-targeted therapeutic strategies have been implemented in the clinic or are in clinical trial including inhibitory immune checkpoints such as CTLA-4 and LAG3. 4 Our data suggest that not all Tregs will be targeted equally since the expression of inhibitory immune checkpoints differs amongst Tregs in one environment. It would be valuable to study patient inter-and intravariability (if heterogeneity is dynamic) regarding Treg heterogeneity since this could aid in improving personalised Treg-based therapeutic strategies; is for example PD-1, LAG3 or a combination, better suited for a patient? 48 Additionally, when aiming to prevent adverse side effects, it might be more specific to target chemokine receptors such as CCR4 in clinical trial 4 because the expression of this chemokine receptor depends on the environment.
In conclusion, our study reveals a heterogeneous population of Tregs at the site of inflammation in JIA. SF Treg differentiate to a classical eTreg profile with a more dominant suppressive or cytotoxic profile that share a similar TCR repertoire, or towards GPR56 + CD161 + CXCL13 + Tregs with a more distinct TCR repertoire. The latter cluster of Tregs is also mirrored in other T-cell subsets at the site of inflammation. Finally, the novel Treg regulon BHLHE40 seems to drive differentiation towards primarily GPR56 + CD161 + CXCL13 + Tregs and JAZF1 towards the classical eTreg phenotype.

Patient samples
For scRNA-sequencing, patients with oligo JIA (n = 3, 3/3 female) were enrolled in the paediatric rheumatology department at the University Medical Center of Utrecht (the Netherlands). The average age was 9.3 years (range 7-12 years) with a disease duration at the time of inclusion of 7 years (range 4-10 years). Two patients were without medication, and one received methotrexate maintenance therapy. HLA-B27 status has not been assessed. For flow cytometry, patients with oligo JIA (n = 17; 40% female, average age 14.2 [2-19 years]) of which 30% had extended and 55% persistent oligoarticular disease were included. Ten patients were without medication, five on methotrexate, two on NSAIDs and three on anti-TNF maintenance therapy. In addition, PB of controls (n = 21; 58% female, average age 39.3 [25-62 years]) were included.
Active disease was defined by physician global assessment of ≥ 1 active joint (swelling, limitation of movement), and inactive disease was defined as the absence hereof. During an outpatient clinic visit, SF was obtained by therapeutic joint aspiration of the affected joints, and blood was withdrawn via vein puncture or an intravenous drip catheter. The study was conducted in accordance with the Institutional Review Board of the University Medical Center Utrecht (approval no. 11-499/C). PB from healthy adult volunteers was obtained from the Mini Donor Service at University Medical Center Utrecht. The research was carried out in compliance with the Declaration of Helsinki. Informed consent was obtained from all the participants and/or from their parents/guardians/legally authorised representatives.
Synovial fluid of JIA patients was incubated with hyaluronidase (Sigma-Aldrich, St Louis, USA) for 30 min at 37°C to break down hyaluronic acid. Synovial fluid mononuclear cells (SFMCs) and peripheral blood mononuclear cells (PBMCs) were isolated using Ficoll Isopaque density gradient centrifugation (GE Healthcare Bio-Sciences, Chicago, USA).

Single-cell mRNA-sequencing
Live CD3 + CD4 + CD25 + CD127 low cells were sorted from fresh SF (Supplementary figure 1a) into 384-well hard shell plates (Bio-Rad, Hercules, USA) with 5 lL of vapour-lock (QIAGEN, Hilden, Germany) containing 100-200 nL of RT primers, dNTPs and synthetic mRNA Spike-Ins and immediately spun down and frozen to À80°C. Cells were prepared for SORT-seq as previously described. 49 Illumina sequencing libraries were then prepared with the TruSeq small RNA primers (Illumina, San Diego, USA) and sequenced single-end at 75 basepair read length with 60 000 reads per cell on a NextSeq500 platform (Illumina). Sequencing reads were mapped against the reference human genome (GRCh38) with BWA.

Single-cell mRNA-sequencing analysis
Quality control was performed in R with the scater package v1.12.2 50 and cells were dropped when the number of genes, number of UMI's and/or the percentage of mitochondrial genes was over 3 median absolute deviations under/above the median. Afterwards, principal component analysis (PCA) outliers were removed with the package mvoutlier v1. The raw data expression matrices were subsequently analysed using Seurat v2-4 51-53 following the outline provided by the distributor (https://satijalab.org/ seurat/). Each data set was normalised, and the cell cycle was regressed out using SCTransform. 54 Thereupon, the SF data sets were integrated with PrepSCTIntegration followed by FindIntegrationAnchors and IntegrateData for SCTransform-processed data.
For dimensionality reduction, first the principal components (PCs) were calculated (RunPCA) and clustering was performed with UMAP (RunUMAP: 30 dimensions; FindNeighbors: clustering resolution of 1). One cluster was removed from further analyses (Supplementary figure 2a) since it was of ambiguous origin (e.g. hybrid, transferred extracellular vesicles or doublets). Subsequent differential gene expression was performed using the MAST test (standard settings) with a P-adjusted value < 0.05 considered statistically significant. For visualisation, the functions DoHeatmap, Dimplot, Featureplot and plot_density were employed. UMAPs were plotted with raw mRNA counts or using the new nebulosa algorithm 55 based on kernel density estimation to handle sparsity of scRNAsequencing data.
Gene set enrichment analysis was performed with gene sets derived from the literature (Ferraro et al. 13 for a human Treg signature, Niedzielska et al. 56 for a shared tissue Treg signature, De Simone et al. 57 for a tumorinfiltrating Treg signature, and Li et al. 41 for a dysfunctional signature of CD4 + T cells) or the MSigDB C7 database (GSE61077 for CD44 hi CD62L lo versus CD44 lo CD62L hi murine Tregs, GSE11057 for naive T cells versus effector memory human CD4 T cells, GSE16835 for CD3/CD28 stimulated versus ex vivo Treg) with AUCell. 14 The cut-off for enrichment of a gene set in a cell was defined using the AUC. Per cluster the proportion of cells enriched for the gene set was calculated and compared with the Chi-square test. Gene ontology pathway analyses implementing the probability density function were performed using ToppFun (https://toppgene.cchmc.org/enrichment.jsp) with as input the differentially expressed genes belonging to each Treg cluster, with a false discovery rate (FDR)-corrected P-value < 0.05 defining significance.
For pseudotime trajectory analysis, Monocle v3 17 was used with a trajectory predicted using standard settings based on the clustering previously performed with Seurat. The principal root node was estimated mathematically with the function get_earliest_principal_node. Network inference analysis employing SCENIC 14 was performed in python v3.6 with Jupyter Notebook v6.1.5 using the UMAP clustering as starting point to define regulons. In short, co-expression based on the raw count data and DNA motif analysis is used to obtain transcription factors and their target genes using standard settings. Activity of these potential TFs and targets (regulons) are analysed per cell, and finally, the scores per cell were combined to compare clusters with each other to define regulators that might drive differentiation within the SF environment. The scRNAsequencing count data generated for this study have been submitted to a public repository on GitHub (https://github. com/lutterl/JIA-synovial-fluid-Tregs-scRNAseq). Raw data files will be made publicly available before publication.

Single-cell TCR-sequencing analysis
RNA-and TCR-sequencing profiling data of single-cell SF Tregs employing 109 genomics were downloaded from GSE160097 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi? acc=GSE160097). TCR-sequencing data were analysed using scRepertoire following the outlined guidelines. 58 In short, T-cell receptor alpha locus (TRA) and TCR beta locus (TRB) data were combined based on the cell barcode. If there was more than one TRA and/or TRB chain detected, the most prevalent chain was selected for integration. RNAsequencing and TCR-sequencing data were combined, and subsequent T-cell clustering was performed with Seurat as described above. TCR-sequencing analysis was performed on the nucleotide level.

Statistical analysis
Statistical analyses were performed with Pearson's Chisquared test, a Friedman's test with Dunn's post hoc, or a one-way ANOVA with Tukey's post hoc test if applicable. Paired data comparisons with missing values were analysed with a mixed-effect model (restricted maximum likelihood). Analyses were performed in Graphpad Prism v7.04 and v8.3 (San Diego, USA), Excel Office v2017 (Washington, USA) and R v3.5.2-4.1.0 (Cincinnati,).