This project used unbiased systems immunology approaches including RNA sequencing, mass cytometry (CyTOF), and TCR repertoire analytical methods to test the hypothesis that the Tregs that respond to injury represent a subset of CD4+FoxP3+ T cells in mice. All experiments were performed using a well-defined mouse burn trauma model following ARRIVE guidelines. Age- and sex-matched mice were randomized into control and experimental groups, with sample sizes chosen based on statistical power analysis and previous experience. Experimental replication for each experiment is indicated in the figure legends. The investigators were not blinded, and no data were excluded from analysis.

C57BL/6 wild-type (stock #000664), Foxp3^DTR (stock #016958) and BALB/c (stock #000651) mice were purchased from Jackson Laboratory (Bar Harbor, ME). Eight to16-weeks old mice were acclimated for at least 1 week before being used for experiments. All procedures performed in this study were reviewed and approved by the Brigham and Women’s Hospital IACUC (2020N000458) and were found to be in accordance with guidelines set by the U.S. Department of Agriculture (Washington, DC) and the National Institutes of Health (Bethesda, MD).

A mouse burn injury model was used to model traumatic injury in mice as previously described (24). Briefly, mice were anesthetized by intraperitoneal (IP) injection with 125 mg/kg ketamine (Fort Dodge Animal Health, Fort Dodge, IA) and 10 mg/kg xylazine (Lloyd Laboratories, Shenandoah, IA). Buprenorphine at 0.6mg/kg was injected subcutaneously at the time of injury. The mice had their dorsal fur shaved and were placed in a plastic mold exposing 20% of their total body surface area. Injury was induced by immersing the exposed part of the dorsum for 9 seconds in a 90˚ C water bath. This approach causes a full-thickness and well-demarcated anesthetic injury due to complete loss innervation. Uninjured mice underwent the same procedure but were exposed to room temperature water for 9 seconds. All animals were resuscitated with an IP injection of 1 mL of 0.9% pyrogen-free normal saline containing ANTISEDAN (atipamezole, Zoetis, US) at 1mg/kg. The mortality from this burn trauma model is less than 2%.

Mice were euthanized by CO₂ asphyxiation. Injury-site draining lymph nodes (axillary, brachial, and inguinal) were harvested and immediately placed in ice-cold culture medium (RPMI 1640 supplemented with 5% heat-inactivated fetal bovine serum (FBS), 1 mM glutamine, 10 mM HEPES, 100 μM nonessential amino acids, penicillin/streptomycin/fungizone, and 50 μM 2-mercaptoethanol, all purchased from Gibco-Invitrogen, Grand Island, NY). To prepare single cell suspensions, the lymph nodes were mechanically minced through sterile 70- μm cell strainers. Cell preparations were washed twice in culture medium and then strained to remove debris.

Fluorescent conjugated antibodies were purchased from BioLegend (San Diego, CA) or eBioscience (Waltham, MA): APC/Cy7-labeled anti-CD3 (145-2C11), APC-labeled anti-CD44 (IM7), eFluor450-labeled anti-FoxP3 (FJK-16s), and PE/Cy7-labeled anti-CD4 (GK1.5), FITC-labeled anti-CD25(PC61). FITC-labeled Anti-Mouse TCR Vβ Screening Panel was purchased from BD bioscience. Zombie NIR™ Fixable Viability Kit (BioLegend) was used for cellular viability staining. Flow cytometry stains were performed in PBS with 1% BSA and 0.1% sodium azide at room temperature. Cells were plated and Fc-block reagent (TruStain FcX Fc receptor blocking reagent, BioLegend, San Diego, CA) was added for 10 min to minimize nonspecific antibody staining before adding cell-surface fluorochrome-labeled antibodies. For intracellular stains, cells were permeabilized using the FoxP3/transcription factor staining buffer set from eBioscience (eBioscience/ThermoFisher Scientific, Waltham, MA). Stained samples were fixed in 4% PFA in PBS, washed by centrifugation, then reconstituted in PBS for flow cytometry analysis on a MACS-Quant Analyzer (Miltenyi Biotech, San Diego, CA). Data analysis was performed using the FlowJo software program (Tree Star, Ashland, OR). The flow cytometry gating schemes for T cells, CD4⁺ and CD4^- T cells, CD44^high and CD44^low Tregs as well as TCRVβ are presented in Supplemental Figures 6 , 7 . For sorting Tregs, the stains were performed in sorting buffer (RPMI 1640 without phenol red supplemented with 0.5% BSA, 1 mM glutamine, 10 mM HEPES, 100 μM nonessential amino acids, penicillin/streptomycin/fungizone, and 50 μM 2-mercaptoethanol). CD4+ T cells were purified by negative selection using Miltenyi MACS CD4⁺GFP⁺ Tregs or CD4⁺GFP⁺CD44^high and CD4⁺GFP⁺CD44^low cells were sorted respectively according to experimental requirements.

CD4⁺CD25^- T cells (Tconv) for Treg functional assays were negatively selected using LS Columns (Miltenyi Biotec, 130-042-401) with MidiMACS™ Separator (Miltenyi Biotec,130-042-302) by staining the cells with Biotin-labeled anti-mouse CD8a (BioLegend 100704), CD25 (BioLegend 102004), TER119 (BioLegend 116204), CD45R/B220 (BioLegend 103204), CD49b (BioLegend 108904), and CD11b (BioLegend 101204). The negatively selected cells were collected and stained with Zombie NIR™ Fixable Viability kit and then CellTrace™ Violet Cell Proliferation Kit (ThermoFisher, C34571) for the Treg functional assay. CD44^high and CD44^low Tregs were sorted respectively from sham or day 7 after injury Foxp3^DTR mice, co-cultured with CellTrace Violet labeled Tconv, and stimulated with mouse T cell activator anti-CD3/CD28 beads (Gibco, 11452D). Purified Tconv (8 x 10⁴) were co-cultured with CD44^high or CD44^low Tregs at 1:1, 1:2, 1:4, 1:8 Treg : Tconv ratios or control Tconv alone. Proliferation of CD3/CD28 bead stimulated Tconv without added Tregs was set as control proliferation with 0% suppression. Proliferation was measured using the CellTrace Violet dilution assay by flow cytometry after 4 days co-culture that detects proliferated cells by reduced fluorescence intensity in cells that undergo one or more cell divisions. Suppression% was calculated as: 100 − ( Percent proliferated cells in Treg : Tconv Percent proliferated cells Tconv only × 100 )

Foxp3^DTR mice were treated with 40ng/kg diphtheria toxin (DT) at 2h and 24h after burn injury to deplete Tregs. Male C57BL/6 mice underwent burn trauma injury. At 7 days after injury, CD4 T cells were purified from injury-site draining lymph nodes (axillary, brachial, and inguinal) into CD44^high or CD44^low Tregs using anti-CD25 and anti-CD44 antibodies. Tregs from injured mice were transfused into the anesthetized injured Foxp3^DTR mice depleted of Tregs at 50,000 cells/mouse by intracardiac injection in 0.2 ml of PBS. Mice were subsequently challenged with intranasal Pseudomonas aeruginosa 1 day after Treg-transfusion. Survival was monitored over a week period. For lung infections, Pseudomonas aeruginosa (ATCC 27853) were grown for 16 hours with gentle agitation in brain-heart infusion (BHI) broth medium at 37°C and were harvested by centrifugation at 450 x g for 10 minutes, then washed once with PBS by centrifugation. Based on absorption spectroscopy measurements (ABS₆₀₀), bacteria were diluted to contain 3-5 x 10⁷ CFUs/ml in PBS. Mice were anesthetized by IP injection using Ketamine/Xylazine at 125/10 mg/kg. Mice were then held with their nares upright, and 40 μl of bacteria suspension was administered by intranasal route (1.2-2 x 10⁶ CFUs). This inoculum dose was found to cause 50% mortality in male C57BL/6 mice over a 7 day period, with deaths first occurring at days 2-3 after infection. Bacteria CFUs were quantified by drop-plating of serial dilutions on Luria broth (LB) agar plates, and colonies were counted the following day after incubating the plates at 37°C.

Lymph node CD4⁺CD44^highFoxP3-GFP⁺ and CD4⁺CD44^lowFoxP3-GFP⁺ T cells or FoxP3-GFP⁺ T cells were purified by FACS sorting at 7 days after injury from injured or uninjured Foxp3^DTR mice (strategy shown in Supplementary Figure S7 ). RNA was purified from 50,000 cells from each sample using the RNeasy Micro Kit (Qiagen, CA) according to manufacturer’s protocol. RNA sequencing (RNAseq) was performed by the Molecular Biology Core Facilities (MBCF) at Dana-Farber Cancer Institute (DFCI). cDNA was synthesized using Takara SMART-Seq v4 PLUS kit and was fragmented to ~150bp by Covaris Adaptive Focused Acoustics^®-AFA^® technology. The cDNA library prepared with Swift 2S™ Turbo DNA library kit from 2ng of cDNA was submitted for next-generation sequencing (NGS).

Analysis of RNA sequencing data was performed using Visualization Pipeline for RNA-seq (VIPER) workflow (25). The VIPER workflow provides sample to sample correlation and sample to feature heatmap plots to judge the correlation and clustering patterns of all samples ( Figures S2A, B ). Gene count stabilized variance of RNAseq data was determined using the vst function of DESeq2 (26). Principal component analysis (PCA) plots were generated based on the variance stabilized counts. The top 2000 genes that were most differentially expressed in CD44^high versus CD44^low Tregs were identified using DESeq2 and input into the gene ontology (GO) molecular function term enrichment analysis webtool (27–29). Heat maps were generated using gene count Z-scores of highly variable cytokines, cell surface markers, and transcription factor genes. Hierarchical clustering of samples was done using the hclust function in R. Differentially expressed genes between Treg subsets and the effects of injury on gene expression were identified by DESeq2 and volcano plots showing differential expression of genes were generated using the Enhanced Volcano R package (30).

Viable CD4⁺CD44^highGFP⁺ and CD4⁺CD44^lowGFP⁺ cells were sorted from lymph nodes harvested at 7-days after injury from burned or uninjured Foxp3^DTR mice (cell numbers are detailed in Supplementary Table 1 ). All RNA extraction, cDNA synthesis, amplification, NGS library preparation and sequencing were performed by iRepertoire, Inc. (Huntsville, USA), RNA was extracted from FACS sorted cells using Qiagen’s RNeasy Micro kit with DNase Treatment. One-third of the total RNA from each sample were used for the construction of mouse TCRα and β chain libraries by reverse transcription polymerase chain reaction (RT PCR) using Qiagen’s OneStep RT PCR mix (Qiagen). First strand cDNA was selected, and unused primer was removed by SPRIselect beads (Beckman Coulter) followed by a second round of binding and extension with the V-gene primer mix. After binding and extension, SPRIselect beads were used to purify the first and second strand synthesis products. Library amplification was performed with a pair of primers that are specific for communal sites engineered onto the 5’ end of the C- and V- primers used in first and second strand synthesis. The final constructed library includes Illumina dual index sequencing adapters, a 10-nucleotide unique molecular identifier region, and an 8-nucleotide internal barcode associated with the C-gene primer. The amplified libraries were multiplexed and pooled for sequencing on the Illumina MiSeq platform using a 600-cycle kit and sequenced as 250 pair-end read. The portion of TCRA and TCRB receptor genes that were sequenced included the second framework region stretching to the beginning of the TCR constant region, including the CDR2 and CDR3 regions.

Sequencing of raw data were analyzed using the previously described iRmap program (11, 31). Sequence reads were de-multiplexed according to both Illumina dual indices incorporated during the amplification process and barcode sequences at the 5’ end of reads from the constant region. Reads were then trimmed according to their base qualities with a 2-basesliding window. If either quality value in this window was lower than 20, this sequence stretching from the window to the 3’ end, was trimmed out from the original read. Trimmed pair-end reads were joined together through overlapping alignment with a modified Needleman-Wunsch algorithm (32). If paired forward and reverse reads in the overlapping region were not perfectly matched, both forward and reverse reads were thrown out without further consideration. The merged reads were mapped using a Smith-Waterman algorithm to germline V, D, J and C reference sequences using an IMGT reference library. To define the CDR3 region, the position of CDR3 boundaries of reference sequences from the IMGT database was migrated onto reads through mapping results, and the resulting CDR3 regions were extracted and translated into amino acids.

Treemap plots were generated by iRweb tools (iRepertoire, Inc. Huntsville, USA). V, D and J gene usage and CDR3 sequences were identified and assigned. In treemap plots, each unique CDR3 is shown as a colored rectangle. The size of each rectangle corresponds to the abundance of each CDR3 within the repertoire and the positioning is determined by the V region usage. To compare the relative diversity of TCR libraries, we used iRweb tools to generate diversity-50 (D50) values, diversity plots, and diversity index (Di) values. D50 is the percent of unique TCR CDR3 sequences that account for the cumulative 50% of the total CDR3s counted in the sample. D 50 =   n u m b e r   o f   u C D R 3   s e q u e n c e s   i n   t o p   50 %   ×   100   n u m b e r   o f   u C D R 3   i n   t o p   10 , 000   C D R 3   s e q u e n c e s . Low D50 values indicate higher clonality. The diversity index (Di) was calculated as: D i = [ 1 −   ∑ ​ n ( n − 1 ) N ( N − 1 )   ] × 100 . TCR diversity index values are proportional to the values where high values indicate high diversity, while low values indicate low diversity. Diversity index plots were generated by plotting the percentage of total TCR reads versus the percentage of unique CDR3 sequences. The curve illustrates the overall diversity of the population with full diversity represented by the dashed line in the middle. The higher the TCR diversity in the cell population, the closer the curve will be to the dashed line in the middle. The D50 values are indicated on these diversity index plots to illustrate the differences in TCRA and TCRB clonality in Treg populations from uninjured and injured mice.

For single cell RNAseq (scRNAseq) and single cell TCRseq (scTCRseq), viable lymph node FoxP3 GFP⁺ cells were sorted from injured or uninjured Foxp3^DTR mice at 7 days after injury. The sorted FoxP3 GFP⁺ cells were washed and resuspended in sorting buffer at a cell concentration of 1000 cells/µL. About 17,000 mouse cells were loaded onto a 10x Genomics Chromium™ instrument (10x Genomics) according to the manufacturer’s recommendations. The scRNAseq libraries were processed using Chromium™ single cell 5’ library & gel bead kit (10x Genomics). Matched scTCRseq libraries were prepared using 2 µL of post cDNA amplification material and Chromium Single Cell V(D)J Enrichment Kit, Mouse T Cell. Quality controls for amplified cDNA libraries and final sequencing libraries were performed using Bioanalyzer High Sensitivity DNA Kit (Agilent). The sequencing libraries for scRNAseq and scTCRseq were normalized to 4nM concentration and pooled using a volume ratio of 4:1. The pooled sequencing libraries were sequenced on Illumina NovaSeq S4 300 cycle platform. The sequencing parameters were: Read 1 of 150bp, Read 2 of 150bp and Index 1 of 8bp. The sequencing data were demultiplexed and aligned to GRCm38 using cell ranger version 3.1.0 pipeline (10x Genomics).

scRNAseq reads were aligned to the mouse reference genome (GRCm38) with the Ensembl GRCM28.91 GTF file using Cell Ranger (10x Genomics; sample statistics in Supplementary Table 3 ). The uninjured and 7D after injury samples were aggregated using the aggr function in Cell Ranger. Cell Ranger gene expression matrices were further analyzed with Seurat (v3.2.2). Matrices were filtered to exclude low-quality cells; cells with < 500 or > 4000 unique features. Cells with > 7% reads mapping to the mitochondrial genome were filtered out (thresholds chosen using analysis shown in Supplementary Figures S4 C, D ). Feature counts were normalized using the NormalizeData function of Seurat. Per standard methods, we identified the 2,000 most variable genes and did PCA on these genes. Gene counts were scaled using the ScaleData function of Seurat. Clustering was done using the FindNeighbors (dims = 1:20) and FindClusters (resolution = 0.5) functions in Seurat. Cluster 6 was then removed as these cells expressed B cells markers. TCR clonotype metadata from the output of Cell Ranger VDJ was added to each cell using custom scripts. Cells with a CDR3 clonotype found in 2 or more cells were labeled as “Expanded”. Clustering was redone with parameters to generate many clusters, FindNeighbors (dims = 1:25) and FindClusters (resolution = 2.0) functions in Seurat. Clusters composed of 15% or more expanded cells were labeled as “Expanded” phenotype. Normalized expression was visualized with t-SNE plot projections and violin plots using Seurat.

For scRNAseq, the raw sequencing data generated by Cell Ranger 3.1.0 were loaded into the third-party tools for secondary analysis: Loupe Cell Browser 4.1.0 and Loupe VDJ Browser 3.0.0 (https://www.10xgenomics.com/). TCR sequencing data were aligned to the mm10 reference genome and RefSeq gene models using Cell Ranger VDJ TCA and TCB sequences from individual cell were used to infer clonotypes. The clonotype comparison feature in Loupe Cell Browser was then applied to pool TCR clonotypes across groups by matching CDR3 amino acids of both TCRA and TCRB. TCR clonotypes were designated as “expanded” if the paired CDR3 sequences were found in 2 or more cells.

CyTOF staining was performed at room temperature as previously described (33). Cells harvested from lymph nodes were incubated with cisplatin live/dead stain for 2 min (Fluidigm) and blocked with TruStain FcX Fc receptor blocking reagent for 10 min before antibody staining. For surface staining, cells were stained with the CyTOF antibody cocktail for 30 min. After washing twice with staining buffer (calcium/magnesium-free PBS, 0.2% BSA, 0.05% sodium azide), cells were fixed and permeabilized using eBioscience™ Fixation/Permeabilization kit (Invitrogen). Subsequently, cells were barcoded using a palladium-based barcode reagent (33). After washing out excess barcode reagent, samples were pooled together and stained with intracellular CyTOF antibodies. Details of all antibodies used for CyTOF staining are listed in Supplementary Table 4 . Single-cell CyTOF data was collected using a Helios mass cytometer (Fluidigm). Data normalization and deconvolution of barcoded staining data was conducted using the normalizer and the single-cell-debarcoder software developed in the Nolan Lab (Stanford) (33, 34). Data analysis was conducted using Cytobank (35). A nonlinear dimensionality reduction algorithm, viSNE, was run on all markers for dimensionality reduction (36). A hierarchical clustering algorithm, SPADE was run on viSNE parameters to cluster cells for phenotypic identification (37).

Statistical analysis was performed using GraphPad Prism 8.4.3 (GraphPad Software Inc. San Diego, CA, USA) or custom R scripts for RNAseq data. Data comparing multiple groups were analyzed by one-way ANOVA with Dunnett statistical hypothesis testing to correct for multiple comparisons and assuming Gaussian distribution of residuals. Data consisting of one variable and two factors were analyzed by two-way ANOVA and the Sidak or Tukey tests were applied for multiple comparison correction as appropriate. For protein expression validation studies by CyTOF, two-tailed unpaired t-tests were applied, and Welch’s correction was used if two populations had different standard deviations. Analyzed data with P <0.05 were considered as significantly different.

We first assessed the function of CD44^high Tregs in a mouse model of burn injury. Injury-site draining lymph nodes were harvested from mice to measure CD44^high and CD44^low CD4⁺FoxP3⁺ T cell (Treg) abundances at 1 and 7 days after burn trauma injury using the scheme shown in Figure 1A . By day 7, there was a significant increase in CD44^high Tregs in injury-site draining lymph nodes compared to uninjured controls ( Figures 1B, C ). In contrast, there was no difference in the number of CD44^low Tregs between injured and control mice ( Figure 1D ). Next, we utilized a fluorescence-activated cell sorting (FACS) approach to purify CD44^high and CD44^low Treg subsets to assess Treg function and for subsequent RNA sequencing studies ( Figure 1E ). In Treg-mediated suppression assays, both CD44^high and CD44^low Tregs significantly suppressed the proliferation of CD4⁺CD25^- T cells (Tconv). The CD44^high Tregs were identified as being much more potent than CD44^low Tregs at suppressing T cell proliferation and injury did not further enhance the suppressive function of CD44^high Tregs ( Figures 1F, G ). However, we did observe that CD44^low Tregs from injured mice were significantly more suppressive at a 1:1 Treg to Tconv cell ratio. Taken together, both Treg subsets are affected by injury showing either expansion (CD44^high) or enhanced Treg suppressive function (CD44^low).

Because infectious complications are a leading cause of death in trauma patients (38, 39), we next tested whether CD44^high or CD44^low Tregs were beneficial or harmful in a secondary infection after burn injury ( Figure S1A ). In brief, Foxp3^DTR mice were Treg depleted by administration of diphtheriae toxin after burn trauma as previously described (40). Treg depletion was validated by flow cytometry ( Figure S1B ). CD44^high and CD44^low Tregs sorted from injury-site draining lymph nodes of burn trauma mice were then adoptively transferred into Treg-depleted mice 2 days after injury. Mice were challenged intranasally with Pseudomonas aeruginosa (1.2-2 x 10⁶ CFU per mouse) the next day and were observed over a week for survival (9, 10, 12). Injured mice given CD44^high Tregs demonstrated a trend towards lower survival (43%) compared to mice given CD44^low Tregs (75%, probability of survival, P=0.0617) ( Figure S1C ). Thus, CD44^high Tregs suppress anti-microbial immune function more so than CD44^low Tregs and their expansion may increase susceptibility to secondary infections following traumatic injury.

Next-generation RNA-sequencing (RNAseq) was used as an unbiased approach to compare gene expression variation and profiles between CD44^high and CD44^low Tregs from injured and uninjured mice. RNA from flow-sorted CD44^high and CD44^low Tregs from injured and uninjured FoxP3-GFP mice were prepared for RNA sequencing ( Figure 1E ). The resulting RNAseq data underwent quality control analysis to generate intra- and inter-group sample variability plots, sample-sample correlation plots, and sample-feature hierarchical clustering plots ( Figures S2A, B ). Principal component analysis (PCA) of sequencing data demonstrates that CD44^high versus CD44^low Tregs separated along principal component 1 (PC1) and accounted for the majority (89%) of gene expression variation ( Figure 2A ). Injury accounted for 5% of the observed variance in PC2. Gene Ontology (GO) analysis by enrichment of functional terms of DE genes indicated that genes associated with cell division and chemokine receptors were enriched in CD44^high as compared to CD44^low Tregs ( Figure 2B ). Given these findings, we specifically compared gene counts of the top differentially expressed chemokine receptors between CD44^high and CD44^low Tregs, which showed distinct differences in gene expression profiles ( Figure 2C ). Heatmaps showing expression levels of the most variably expressed cytokine, cell surface markers, and transcription factor were generated from RNAseq data to illustrate some of the key markers that distinguished these Treg subsets ( Figure 2D ). CD44^high Tregs expressed higher levels of known Treg-related genes such as Icos, Ctla4, Il10ra, Il10, Itgae, Tigit, Fgl2, Havcr2, Cd83, Nt5e, and Klrg1 than CD44^low Tregs. Moreover, some notable cytokine and chemokine receptor genes such as Il22ra2, Tgfbr1, Ildr1, Il18r1, Il10ra, Il3ra, Il17rb, Il23r, Il1r2, Il12rb, Ccl5, Cxcl10 Cxcl2, Cxcr3 were highly expressed in CD44^high Tregs, while others like Il2ra, Il6ra, Tgfbr3, Il4ra and Trnfsf8 were highly expressed in CD44^low Tregs. CD44^high Tregs from injured mice had higher expression levels of Cd22, Cxcr5 and Il7r, as well as genes encoding components of TNF receptor signaling pathways, such as Tnfrsf18, Tnfrsf9, and Tnfrsf4, when compared to uninjured mice. Tregs from injured mice showed higher counts of Cd69 and Ifnr1 in both CD44^high and CD44^low groups. Notable transcription factor genes highly expressed in CD44^high subsets were Arnt2 (involved in responses to environmental stimuli), Prdm1 (tissue-resident memory T cell signaling), Atf6 (endoplasmic reticulum [ER]stress sensor), Mybl2 (cell survival, proliferation, and differentiation); Crem (cAMP responsive element), T cell related genes (Tbx21, Gata3, Rora, Irf4, Maf, Id2, Nfil3), as well as cell cycle-related genes and some zinc finger proteins ( Figure 2C ). There were 24 differentially expressed (DE) genes increased and 25 DE genes decreased due to injury in the CD44^high subset, and 35 increased and 15 decreased in the CD44^low subset ( Figure 2D ). DE genes in CD44^high versus CD44^low Tregs in both the uninjured and injured groups are shown in Figures S2C, D . Among genes with a P value < 10^-6, there were 1006 genes with an increase of 2-fold or more and 486 genes with a decrease of 2-fold or more in the uninjured group CD44^high Tregs. In the injured group, gene counts were 743 and 313, respectively. Upon burn injury, 24 genes were upregulated and 25 downregulated in the CD44^high Tregs, for the CD44^low Tregs, 35 upregulated and 15 downregulated ( Figure 2E ). Taken together, the bulk RNA sequencing results indicate that CD44^low and CD44^high Tregs are transcriptionally distinct populations.

The experimental scheme for TCRα and TCRβ repertoire analysis is shown in Figure 3A . Treg TCR clonality and diversity were identified by RNA sequencing using the iRepertoire analytical platform. A summary of the TCR sequencing dataset is provided in Supplementary Table 2 . Without injury, CD44^high Tregs were found to have higher TCRα/TCRβ clonality than CD44^low Tregs as determined lower D50 and diversity index (Di) values ( Figure 3B ). The D50 value indicates the percentage of TCR CDR3 sequences that account for 50% of the total unique CDR3 sequence reads, so low D50 values would indicate high clonality (41, 42). Interestingly, the D50 value decreased for TCRα in CD44^high Tregs (uninjured = 9.7, injured = 5.5), but D50 for TCRβ did not change following injury suggesting that TCRα clonality is influence by injury. The Di values and curves support that CD44^high Tregs have lower TCR diversity than CD44^low Tregs ( Figure 3B ). The Di values were much lower for CD44^high Tregs (TCRα/13.3, TCRβ/10.8) as compared to CD44^low Tregs (TCRα/31.3, TCRβ/35.2). The shape of the diversity curves shown in Figure 3B illustrates that injury increased TCR diversity in CD44^high Tregs by altering the curve to be closer to the diagonal dashed line, which identifies full TCR diversity. Accordingly, CD44^low Tregs showed higher diversity than CD44^high Tregs since the diversity curves are much closer to the diagonal line. Treemap plots illustrate that injury increased the diversity of CD44^high Tregs, which is likely due to emergence of new CD44^high Tregs or movement of CD44^high Tregs from other tissues into lymph nodes following injury ( Figures 3C, D ). Plots comparing the top 10 TCRA and TCRB sequence percentages between CD44^high and CD44^low Tregs further support differences in TCRα and TCRβ clonality between these Treg populations ( Figures 3E, F ).

Next, specific TCR-Vβ chain expression on Tregs prepared from 2 different inbred lines of injured mice were compared by flow cytometry. Among the 14 different TCR-Vβ chains tested, Tregs expressing TCR-Vβ 3, 4, and 6 were expanded in the CD44^high Treg population from injured C57BL/6 mice and TCR-Vβ 4, 6, 8.1/8.2, 8.3, and 14 were expanded CD44^high Treg from injured BALB/c mice. No TCR-Vβ expansion was observed on CD44^low Tregs or conventional CD4⁺ T cells in either inbred mouse lines ( Figures 3G, H ; Figures S3A, B ). Thus, the oligoclonal expansion of CD44^high Tregs is not restricted to a single inbred mouse strain with different TCR repertoires and only CD44^high Tregs expanded in response to injury.

Paired single-cell TCR and RNA sequencing on FACS sorted FoxP3-GFP⁺ cells was performed to identify αβ TCR clonotypes and matched transcriptome expression profiles in Tregs from injured and uninjured mice ( Figure 4A ). Quality control was performed on single-cell RNA sequencing data to exclude doublets and dying cells ( Figure S4 ). Consistent with above TCR repertoire analysis done by the iRepertoire sequencing approach, we observed 182 expanded clonotypes (2 or more paired TCRαβ clones) in Tregs from burn injured mice and 83 clonotypes in uninjured mice ( Figure 4B ). There were no overlapping paired TCR clonotypes among the top 10 clonotypes identified in CD44^high or CD44^low Tregs from injured or uninjured mice. However, the frequency of Tregs with expanded TCR clonotypes was markedly higher in injured mice as compared to uninjured mice ( Table 1 ). We next classified Treg scRNAseq clusters that have 15% or more of their cells with 2 or more paired TCR clonotypes as having an “expanded phenotype” ( Figure S4E, F ). Treg-related genes such as Klrg1, S100a4, Tigit, S100a6, Icos and Itgae were found to be significantly differentially expressed in expanded as compared to non-expanded Tregs ( Figure 4C ). Moreover, t-SNE plots showed that Klrg1, S100a4, Tigit, S100a6, Icos and Itgae expression overlapped in the areas of the t-SNE cell atlas where cells with expanded TCR clonotypes were found ( Figure 4D ). Expression Cd44 was also found to be correlative with the expanded populations. Gene ontology (GO) analysis of the top differentially expressed genes between the expanded and non-expanded populations found the most significant GO term hits to be SRP-dependent translational protein targeting to membrane, T cell activation, negative regulation of immune system process, and T cell selection ( Figure 4E ). These results suggest that injury expanded Tregs acquire a highly activated phenotype.

Density contour plots of CyTOF staining data analyzed by tSNE for dimensionality reduction shows injury-induced increases in cell density in regions identified as FoxP3⁺ and Helios⁺ CD4⁺ T cells ( Figures 5A, C ). CyTOF data was clustered using the SPADE algorithm and cluster 5 was identified as CD44^high Treg and cluster 13 as CD44^low Treg clusters based on CD44, ICOS, CD25, Foxp3, CLTA-4 and OX40 expression patterns ( Figure 5B and Figures S5A, B ). The abundance of cells within the CD44^high Treg cluster increased at 7 days after injury but there was no significant increase in the CD44^low Treg cluster ( Figure 5D ). Expression of KLRG1, ICOS, CTLA-4, ITGAE (CD103), Helios, Galectin-3, PYCARD, ICAM-1, CXCR6, and CRLF2 proteins in the CD44^high Treg cluster was significantly higher than that in the CD44^low Treg cluster ( Figure 5E ). Deeper analysis of marker expression levels on Tregs showed that injury further upregulated Helios, Galectin-3 and PYCARD expression in the CD44^high Treg cluster ( Figure 5F ).