Mice

Foxp3Thy1.1, Foxp3GFP-DTR, Foxp3CreER-GFP and Il10fl mice have been previously described and were maintained in-house8,33,70,71. Gt(ROSA)26SorLSL-YFP and Rorcfl have been previously described and were purchased from Jackson Laboratories34,72. Tcrafl and Maffl have been previously described55,73. Maffl mice were provided by D. R. Littman and Tcrafl mice were provided by M. Schmidt-Supprian. Il10FM mice were generated by intercrossing Foxp3Thy1.1, Gt(ROSA)26SorLSL-YFP and Il10tdTomato-CreER mice (see below) to homozygosity for each allele. Littermates were used in all experiments and were distributed among experimental groups evenly whenever possible, with different experimental groups co-housed. In experiments with different genotypes, all genotypes were represented in each litter analyzed. All mice were maintained at the Research Animal Resource Center for Memorial Sloan Kettering Cancer Center (MSKCC) and Weill Cornell Medicine under specific-pathogen-free conditions, with controlled humidity and temperature, a 12 h/12 h light/dark cycle and ad libitum access to diet (LabDiet 5053) and reverse-osmosis-filtered water. For studies in which treatments were given in drinking water, the same reverse-osmosis-filtered water was used as the vehicle. All studies were under protocol 08-10-023 and approved by the MSKCC Institute Institutional Animal Care and Use Committee. All animals used in this study had no previous history of experimentation and were naive at the time of analysis. Both sexes were used in all experiments unless otherwise noted, as no sex differences in IL-10 expression were detected.

Generation of Il10 tdTomato-CreER and Il10 tdTomato-Cre mice

Il10tdTomato-CreER mice were generated by insertion of a targeting construct into the Il10 locus by homologous recombination in embryonic stem cells on the C57BL/6 background. The targeting construct was generated by inserting a sequence containing exons 2–5 of the Il10 gene into a plasmid backbone containing a PGK promoter driving expression of diphtheria toxin A subunit followed by BGHpA sequence (modified PL452 plasmid). A SalI restriction enzyme site was simultaneously engineered into the Il10 3′ UTR between the stop codon and the polyadenylation site. The Clontech Infusion HD Cloning system was used to generate in the pUC19 plasmid backbone sequence encoding (in order from 5′ to 3′) encephalomyocarditis virus IRES; tdTomato; T2A self-cleaving peptide from Thosea asigna virus; Cre recombinase fused to the estrogen receptor ligand binding domain (CreER); followed by a FRT site-flanked PGK-Neomycin resistance gene (Neo)-BGHpA cassette. The IRES-tdTomato-T2A-CreERT2-FRT-Neo-BGHpA-FRT sequence was PCR-amplified and inserted into the SalI site in the Il10 3′ UTR in the modified PL452 backbone. The resulting plasmid was linearized with the restriction enzyme NotI before electroporation into embryonic stem cells. Il10tdTomato-CreER mice were bred to Gt(ROSA)26SorFLP1 mice (MSKCC Mouse Genetics Core) to excise the Neo cassette and backcrossed to C57BL/6 mice to remove the Gt(ROSA)26SorFLP1 allele. Il10tdTomato-Cre mice were generated in an identical manner except that the targeting vector contained a codon-optimized NLS-Cre encoding sequence after the T2A.

Generation of Foxp3 LSL-DTR mice

Foxp3LSL-DTR mice were generated by Biocytogen. First, a guide RNA targeting the 3′ UTR of the Foxp3 gene was designed and validated (GGAAAGTTCACGAATGTACCA). Then, a targeting vector was constructed including 1,400 bp homology upstream and downstream of an SspI site in the Foxp3 3′ UTR. The following sequence was inserted into the SspI site: loxP-IRES-thy1.1-pA-loxP-IRES-DTR-eGFP. Cas9 protein, in-vitro-transcribed sgRNA and the targeting vector were then micro-injected into C57BL/6N zygotes. Founder pups were then bred and confirmed to have the proper integration by PCR and Southern blot analysis.

Mouse treatments

For tamoxifen treatment, mice were gavaged with 8 mg tamoxifen dissolved in 200 μl corn oil (Sigma-Aldrich). Tamoxifen was dissolved by gentle agitation at 37 °C overnight. Aliquots were frozen (−80 °C) and thawed as needed throughout the experiments. We found that freezing and storing at −80 °C rather than −20 °C greatly reduced the tendency of the tamoxifen to precipitate when thawed. For Il10iΔTCR experiments, mice were treated on days 0 or 11 and analyzed on day 21. For Foxp3iΔIl10 experiments, mice were treated on days 0 and 2 and analyzed on day 16. For Il10iΔRorc, Il10iΔMaf and Foxp3iΔIl10 (long-term) experiments, mice were treated on days 0, 4, 11, 18, 25 and 32 and analyzed on day 35. DT was reconstituted in sterile PBS at 1 mg ml−1 and frozen at −80 °C in single-use aliquots. Aliquots were thawed and diluted in 995 μl PBS. For inactivated control (bDT), this 1 ml dilution was heated at 95–100 °C for 30 min. Both active and control DT were filtered through 0.22 μm syringe-driven filters. Mice were injected intraperitoneally with 200 μl of this dilution for the first two doses (1,000 ng DT) or with 200 μl of a 1:1 dilution with PBS for subsequent doses (500 ng DT), except for experiments depicted in Extended Data Fig. 9h–k, in which 1,000 ng DT was administered for each dose. Bleomycin was dissolved in PBS at a concentration of 5.7 U ml−1, sterile-filtered and frozen (−80 °C) in single-use aliquots. Aliquots were diluted with sterile PBS immediately before use. Mice were anesthetized with isofluorane (3% in O2, 3 l min−1; Covetrus), and 0.1 U bleomycin in 35 μl PBS was administered intranasally using a micropipette. Mice were exposed to isoflurane for at least 5 min before delivery of bleomycin, and the mouth was gently pressed shut during delivery to prevent swallowing. Bleomycin was given drop-wise, with pauses between drops to ensure inhalation. For antibiotic treatment, a solution of 1 g l−1 ampicillin sodium salt, 1 g l−1 kanamycin sulfate, 0.8 g l−1 vancomycin hydrochloride, 0.5 g l−1 metronidazole and 2.5 g l−1 sucralose (Splenda) was prepared in acidified, reverse-osmosed water and sterile-filtered. The control solution contained only sucralose but was otherwise treated similarly. Solutions were replaced every 7 days for the duration of the experiment. For chemically induced colitis (Fig. 1), 15 g of DSS salt (molecular weight, ~40,000) was dissolved in 50 ml distilled deionized water and then sterile-filtered. This solution was then diluted in acidified, reverse-osmosed water, resulting in a final concentration of 3% (w/v) DSS. Control groups received the same amount of sterile-filtered distilled deionized water diluted into acidified, reverse-osmosed water. For these experiments, female mice were used, as male mice proved to be highly sensitive to even lower concentrations of DSS74. For chemically induced colitis (Fig. 7 and Extended Data Fig. 6), 7.5 g of DSS (molecular weight, ~40,000) was dissolved in 50 ml distilled deionized water and then sterile-filtered. For Foxp3CreERIl10fl/KO male mice (Extended Data Fig. 6c–e), 5 g of DSS was dissolved in 50 ml distilled water and then sterile-filtered. The stock solution was then diluted in 450 ml of acidified, reverse-osmosed water, resulting in a final concentration of 1.5% or 1% (w/v) DSS. The solution was replaced after 7 days. Foxp3CreERIl10fl/KO mice were orally administered with two doses of 8 mg tamoxifen dissolved in 200 μl corn oil (Sigma-Aldrich) 48 h apart and, 7 days after the last dose of tamoxifen, they were administered 1.5% w/v (females) or 1% w/v (males) DSS in drinking water, a relatively low dose which causes minimal weight loss in control Il10CreFoxp3LSL-DTR mice treated with bDT.

Cell isolation for flow cytometry

Mice were injected retro-orbitally with 1.5 μg anti-mouse CD45.2 (Brilliant Violet 510 conjugated; BioLegend, 109838) in 200 μl sterile PBS 3 min before the mice were killed to label and exclude blood-exposed cells. All centrifugations were performed at 800g for 3 min at 4 °C. SLOs were dissected and placed in 1 ml wash medium (RPMI 1640, 2% FBS, 10 mM HEPES buffer, 1% penicillin–streptomycin, 2 mM l-glutamine). Tissues were then mechanically disrupted with the back end of a syringe plunger and then passed through a 100 μm, 44% open area nylon mesh. For skin and lung, both ears and all lung lobes were collected. Ears were peeled apart to expose the dermis and cut into six total pieces. Tissues were then placed in 5 ml snap-cap tubes (Eppendorf, 0030119401) in 3 ml wash medium supplemented with 0.2 U ml−1 collagenase A, 5 mM calcium chloride and 1 U ml−1 DNase I along with three ¾ inch ceramic beads (MP Biomedicals, 116540424-CF). The tubes were shaken horizontally at 250 RPM for 45 min at 37 °C for the lung and for two rounds of 25 min for skin, replacing collagenase solution in between. Digested samples were then passed through a 70 μm strainer (Milltenyi Biotec, 130-095-823) and centrifuged to remove the collagenase solution. Lungs were then treated with ACK buffer (155 mM ammonium chloride, 10 mM potassium bicarbonate, 100 nM EDTA pH 7.2) to lyse red blood cells and then washed by centrifugation in 40% Percoll (ThermoFisher, 45-001-747) in wash medium to remove debris and enrich for lymphocytes. For colon, the cecum and large intestine were dissected and, after the removal of fat and the cecal patch, opened longitudinally and vigorously shaken in 1× PBS to remove luminal contents. Colon tissue was then cut into 1–2 cm pieces, placed in a 50 ml screw-cap tube with 25 ml wash medium supplemented with 5 mM EDTA and 1 mM dithiothreitol and shaken horizontally at 250 RPM for 15–20 min at 37 °C. After a 5 s vortex, epithelial and immune cells from the epithelial layer were removed by filtering the suspension through a tea strainer. The remaining tissue was placed back in 50 ml tubes, washed with 25 ml wash medium, strained again and replaced in 50 ml tubes. Then, 25 ml wash medium supplemented with 0.2 U ml−1 collagenase A, 4.8 mM calcium chloride and 1 U ml−1 DNase I was added along with four ¾ inch ceramic beads, and tissues were shaken horizontally at 250 RPM for 35 min at 37 °C. The suspension was then passed through a 100 μm strainer, centrifuged to remove debris and collagenase solution and then washed by centrifugation in 40% Percoll in wash medium. The small intestine was processed with the same steps as the colon, except that the pieces were cleaned by shaking in corn starch and then rinsed with PBS before EDTA treatment. All enzymatically digested samples were washed by centrifugation in 5 ml wash medium.

Flow cytometry

To assess cytokine production after ex vivo restimulation, single-cell suspensions were incubated for 4 h at 37 °C with 5% CO2 in the presence of 50 ng ml−1 PMA and 500 ng ml−1 ionomycin with 1 μg ml−1 brefeldin A and 2 μM monensin to inhibit endoplasmic reticulum and Golgi transport. For flow cytometric analysis, cells were stained in 96-well V-bottom plates with antibodies and reagents used at concentrations indicated in Supplementary Table 1. All centrifugations were performed at 900g for 2 min at 4 °C. Staining with primary antibodies was carried out in 100 μl for 25 min at 4 °C in staining buffer (PBS, 0.2 % (w/v) BSA, 2 mM EDTA, 10 mM HEPES, 0.1% (w/v) NaN3). Cells were then washed with 200 μl PBS and then concurrently stained with Zombie NIR Fixable Viability dye and treated with 20 U ml−1 DNase I in DNase buffer (2.5 mM MgSO4, 0.5 mM CaCl2, 136.9 mM NaCl, 0.18 mM Na2HPO4, 5.36 mM KCl, 0.44 mM KH2PO4, 25 mM HEPES) for 10 min at room temperature (18–23 ºC). Cells were washed with 100 μl staining buffer, resuspended in 200 μl staining buffer and passed through a 100 μm nylon mesh. For cytokine staining, cells were fixed and permeabilized with BD Cytofix/Cytoperm per the manufacturer’s instructions. Intracellular antigens were stained overnight at 4 °C in 1× Perm/Wash buffer. Samples were then washed twice in 200 μl 1× Perm/Wash buffer, resuspending each time, resuspended in 200 μl staining buffer and passed through a 100 μm nylon mesh. All samples were acquired on an Aurora cytometer (Cytek Biosciences) and analyzed using FlowJo (v.10) (BD Biosciences).

Histopathological analysis

Sections of colon (~1 cm) were fixed in 4% PFA for >48 h. Tissue embedding, sectioning and staining was carried out by Histowiz Inc. A blinded pathologist scored sections based on the Simplified Geboes Score rubric described in Supplementary Table 2 (ref. 63).

Flow cytometric identification of immune cell populations

Generally, the following populations were identified with the associated markers (all immune cells were first gated as CD45+ and ZombieNIR– and excluded doublets):

Treg cells: CD90.2+CD5+SSCloFSCloTCRβ+TCRγδ–CD4+CD8α–Thy1.1+

TH cells (CD44hiCD4+ cells): CD90.2+CD5+SSCloFSCloTCRβ+TCRγδ–CD4+CD8α–Thy1.1–CD44+

Macrophages: CD11b+CD64+CD90.2–CD19–NK1.1–Gr-1–SiglecF–Ly6C–/lo

Monocytes: CD11b+CD64–/loCD90.2–CD19–NK1.1–Gr-1–SiglecF–Ly6Chi

Plasma cells: CD19loCD44hiCD64–CD11b–/loCD11c–/loCD90.2–NK1.1–SiglecF–Gr-1–

Germinal center B cells: CD19hiCD44loCD73+IgD–CD64–CD11b–/loCD11c–/lo CD90.2–NK1.1–SiglecF–Gr-1–

γδT cells: CD90.2+SSCloFSCloTCRβ–TCRγδ+

CD8eff cells (CD44hiCD8+ T cells): CD90.2+SSCloFSCloTCRβ+TCRγδ–CD4–CD8α+CD44hiCD62L–

Natural killer cells: NK1.1+CD90.2+/–SSCloFSCloTCRβ–TCRγδ–CD19–CD64–CD11b–CD127–

Neutrophils: Gr-1+CD11b+CD64–/loCD90.2–CD19–NK1.1–SiglecF–

Eosinophils: SiglecF+CD11b+CD64–/loCD90.2–CD19–NK1.1–

Cell sorting for sequencing

Cell isolation was performed as described above, except that samples were not washed with 40% Percoll. Staining was performed as described above, except the buffer contained 2 mM l-glutamine and did not contain NaN3, with the staining volume adjusted to 500 μl, washes adjusted to 5 ml and staining performed in 15 ml screw-cap tubes. ‘Hash-tag’ antibodies (1 μg; BioLegend, 155801, 155803, 155805, 155807) were added to the extracellular antigen stain for scRNA-seq sorting, and scRNA-seq sort samples were not treated with DNase I. Samples were resuspended in wash buffer supplemented with 5 mM EDTA for sorting. Samples were double sorted, with the first sort enriching for all Thy1.1+ Treg cells and the second sort separating Il10neg, Il10recent and Il10stable or tdTomato+ versus tdTomato– cells. For bulk RNA-seq, samples were sorted directly into Trizol-LS per the manufacturer’s instructions in 1.5 ml microcentrifuge tubes. For scRNA-seq, samples were sorted into PBS with 0.04% BSA (w/v) in 1.5 ml Protein LoBind tubes (Eppendorf, 0030108442). For ATAC–seq, samples were sorted into wash medium in 1.5 ml Protein LoBind tubes. All sorting was performed on an Aria II (BD Biosciences).

Preparation of reference genome

The mm39 mouse genome assembly and NCBI RefSeq annotation information (GTF file) were downloaded from the UCSC Genome browser75,76,77,78. To account for the presence of the Il10tdTomato-CreER, Foxp3Thy1.1 and Gt(ROSA)26SorLSL-YFP targeted mutations, the corresponding sequences were inserted into the appropriate locations of the mm39 genome using the ‘reform’ script, creating the ‘reformed mm39’ genome79. The GTF file was modified to appropriately extend the Il10, Foxp3 and Gt(ROSA)26Sor transcript and gene annotations and to shift all other affected annotations, resulting in a ‘reformed GTF’ using a custom R script, relying on the ‘GenomicRanges’ and ‘rtracklayer’ packages80,81,82. The reformed mm39 and reformed GTF were used for bulk RNA-seq and ATAC–seq alignment and analyses after generating a STAR genome index using STAR (v.2.7.3a)83.

Bulk RNA-seq

A total of 5,000 cells were sorted per population per replicate for bulk RNA-seq, with each replicate pooled from two mice. RNA was extracted and libraries were prepared using SMARTer Stranded RNA-Seq Kits according to the manufacturer’s protocols (Takara) by the Integrated Genomics Operation (IGO) Core at MSKCC. Paired-end 50 bp reads (20–30 million per sample) were sequenced on an Illumina HiSeq 3000 by IGO.

Bulk RNA-seq data processing

Samples were processed and aligned using Trimmomatic (v.0.39), STAR (v.2.7.3a) and Samtools (v.1.12), with the following steps, where Sample represents each Il10neg, Il10recent or Il10stable replicate83,84,85.

TrimmomaticPE Sample_R1.fastq.gz Sample_R2.fastq.gz -baseout Sample.fastq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

STAR–runThreadN 6–runMode alignReads–genomeLoad NoSharedMemory–readFilesCommand zcat–genomeDir mm39_100_RNA–readFilesIn Sample_1P.fastq.gz Sample_2P.fastq.gz–outFileNamePrefix Sample–outSAMtype BAM Unsorted–outBAMcompression 6–outFilterMultimapNmax 1–outFilterMismatchNoverLmax 0.06–outFilterMatchNminOverLread 0.35–outFilterMatchNmin 30–alignEndsType EndToEnd

samtools sort -@ 4 -n -o Sample.bam SampleAligned.out.bam

samtools fixmate -@ 4 -rm Sample.bam Sample.fixmate.bam

samtools sort -@ 4 -o Sample.resort.bam Sample.fixmate.bam

samtools markdup -@ 4 -l 1500 -r -d 100 -s Sample.resort.bam Sample.duprm.bam

samtools index -@ 4 -b Sample.duprm.bam

This procedure resulted in the retention of all uniquely aligning reads, with PCR and optical duplicates removed, to be used for downstream analysis. Reads aligning to genes derived from the reformed GTF were then counted using a custom R script relying on the ‘GenomicAlignments’, ‘GenomicRanges’ and ‘GenomicFeatures’ packages with default counting parameters81. Differential expression analysis was carried out using the ‘DESeq2’ package, with the formula ‘~ Celltype + Replicate’, in which Celltype was either Il10neg, Il10recent or Il10stable and replicates were the separate samples from which each of the three populations were sorted86. Fragments per kilobase mapped (FPKM) normalized counts were extracted using the fpkm function of DESeq2. Differential expression analysis and statistical testing were performed for all pairwise comparisons of ‘Celltype’: Il10neg, Il10recent and Il10stable. Differential expression analysis was performed on all genes, but genes with FPKM counts below the mean FPKM count of Cd8a (a gene functionally not expressed in Treg cells), genes with zero counts in the majority of samples or genes corresponding to immunoglobulin or TCR variable, diversity or junction segments were eliminated from subsequent analyses. This process did not remove any significantly differentially expressed genes except immunoglobulin or TCR variable, diversity or junction segments, whose differential expression was not interpretable. K-means clustering was performed with R using per-gene Z-score-normalized counts of genes differentially expressed (adjusted P < 0.05) in any pairwise comparison between the three cell populations, with seven clusters chosen based on preliminary hierarchical clustering. TCR-activated and repressed genes were defined as genes that lost and gained expression in Treg cells ablated of the Tcra gene52.

scRNA-seq

Uniquely ‘hash-tagged’ samples from different tissues were pooled during sorting as separate tdTomato+ and tdTomato– samples or total Treg cells from large intestine lamina propria (LILP) and small intestine lamina propria (SILP) (Thy1.1+CD4+TCRβ+). The tdTomato+ sample had 48,000 cells (35,000 from LILP; 1,100 from lung; 11,000 from mesLN; and 900 from mediastinal lymph node (medLN)) and the tdTomato– sample had 60,000 cells (30,000 from LILP; 10,000 each from lung, mesLN and medLN) or total 16,537 (Thy1.1+CD4+TCRβ+) Treg cells from LILP and SILP. Samples were centrifuged and resuspended in 30 μl PBS with 0.04% BSA (w/v). Libraries were then prepared following the 10× Single Cell 3′ Reagent Kit (v.3) or 5′ kit with V(D)J enrichment for immune profiling (10× Genomics) following the manufacturer’s instructions, incorporating the BioLegend TotalSeq-A hash-tag oligonucleotide (HTO) protocol. Samples were sequenced on an Illumina NovaSeq platform by IGO.

scRNA-seq processing

Reads from the tdTomato+ and tdTomato– samples were processed, aligned to the mm39 genome and demultiplexed using Cell Ranger software (10× Genomics, v.7.0) with default parameters. Reads for the HTOs of the tdTomato+ and tdTomato– samples were processed and demultiplexed using Cell Ranger software with default parameters. Samples were further processed and analyzed with a custom R script relying on the ‘Seurat’ (v.4) package87. First, genes detectable in fewer than 0.1% of cells were removed. Second, HTO identities (that is, LILP, lung, mesLN, medLN) were assigned using the HTODemux function, and cells without an unambiguous HTO identity or those determined to be a doublet were excluded88. Then, cells with mitochondrial genes accounting for >10% of gene counts, presumed to be dead or dying, as well as cells in the top or bottom 2% of total counts were eliminated. This latter cutoff was chosen based on a percentile rather than an arbitrary absolute value to account for different cell numbers and different median unique molecular identifier counts across the two samples. Afterwards, tdTomato+ and tdTomato– samples were merged and analyzed together. First, the top 2,000 variable genes were identified and scaled. A principal component analysis (PCA) was performed on these genes and the top 30 principal components were used to assign k-nearest neighbors, generate a shared nearest-neighbor graph and then optimize the modularity function to determine clusters, at resolution = 0.5 (ref. 89). Based on these original clusters, a small population of cells dominated by high type I interferon signaling was excluded, and subsequent analyses were performed only on cells with LILP or mesLN HTO identities. The remaining cells had gene counts scaled again. A PCA was performed on the 3,000 most variable genes and the top 30 principal components were used to assign k-nearest neighbors, generate a shared nearest-neighbor graph and then optimize the modularity function to determine clusters, at resolution = 0.5. The shared nearest-neighbor graph was used as input for the Python-based algorithm ‘Harmony’ (600 iterations) to generate a two-dimensional force-directed layout for visualization90. The 30 principal components were used as input for the Python-based algorithm ‘Palantir’ to determine ‘pseudotime’ and ‘entropy’ values42. To reconcile scRNA-seq clusters and bulk RNA-seq populations, the ‘Seurat’ function AddModuleScore was used to assign scores for each bulk gene cluster to each cell. Mean scores for every cell cluster were calculated. For enrichment testing, the phyper function of R was used, in which q represents genes in a given bulk RNA-seq k-means cluster and also significantly over-expressed or under-expressed in a given scRNA-seq subset (combination of cell cluster and tdTomato+ or tdTomato- identity); m represents genes in a given bulk RNA-seq k-means cluster; n represents all other genes with detectable expression in a given scRNA-seq subset; and k comprises all genes significantly over-expressed or under-expressed in a given scRNA-seq subset. For determining the relationship between IL-10-expressing Treg cell subsets, we performed paired scRNA-seq and V(D)J-seq (Extended Data Figs. 9 and 10). Reads from the LILP and SILP samples were processed, aligned to the custom mm10 genome that contained Il10tdTomato-CreER, Foxp3Thy1.1 and Gt(ROSA)26SorLSL-YFP targeted mutations and demultiplexed using Cell Ranger software (10× Genomics) with default parameters. Reads for the HTOs of the LILP and SILP samples were processed and demultiplexed using Cell Ranger software with default parameters. Samples were further processed and analyzed with a custom Python script using the ‘scanpy’ package91. First, cells without an unambiguous HTO identity (LILP or SILP) or those determined to be a doublet were excluded. Next, cells with mitochondrial genes accounting for >5% of total genes were considered dead or dying and were removed. All genes encoding ribosomal proteins and genes expressed in less than 0.1% of cells were also removed. Then, using log-normalized data, the top 3,000 variable genes were identified to perform PCA. The top 100 principal components were used to generate a shared neighbors graph and a uniform manifold approximation and projection visualization. Initial clustering was performed using the ‘leiden’ function of scanpy with resolution = 1. Cells enriched in Malat1 expression and with low library size were deemed as low-quality and were removed92. Gene counts of the remaining cells were scaled again, and a uniform manifold approximation and projection embedding was generated after calculating a PCA with the top 3,000 variable genes and creating a neighbors graph with 100 principal components. Unsupervised clustering was performed using resolution = 0.75. Bulk RNA-seq gene cluster (bI–bIV) (Fig. 2c) enrichment score for each annotated scRNA-seq Treg cell cluster was calculated using the ‘score_genes’ function within scanpy. Differential gene expression analysis was performed using the Python-based ‘rpy2’ and R-based ‘MAST’ packages. Up to the top 50 differentially expressed genes (adjusted P < 0.05 and log10FC > 0.5) by either Gata3+Il10+ or Rorc+Il10+ Treg cells from both tissues were plotted. To determine the extent of transcriptional similarity between and predict the developmental trajectory of annotated LILP or SILP Treg cell clusters, the Python-based ‘PAGA’93 analysis was performed using distances computed on a diffusion map. To this end, after removing TCR-related genes, a diffusion map was created using the ‘diffmap’ function of scanpy and a PAGA map was generated using the ‘paga’ function of scanpy. ‘Pseudotime’ and ‘entropy’ values were calculated using the ‘Palantir’ algorithm with 30 principal components. Clonotypes with identical nucleotide-level CDR3 regions (ranging from 1–3 matching chains) were called using the built-in Cell Ranger enclone software, with 4,693 out of 7,868 cells (~60%) that had passed upstream RNA-level quality control being assigned to a clone. To measure clonal relatedness between phenotypic clusters, we computed the Jaccard overlap, defined as:

$$J(_,_)=\,\frac_\cap _|}_|+|_|-|_\cap _|}$$

where \(_\) is the set of clonotypes belonging to phenotype 1 and \(_\) is the set of clonotypes belonging to phenotype 2. We computed both the pooled overlap and the mouse-level overlaps. To measure the reproducibility of this clonal structure, we measured the correlation between the overlap matrices of pairs of mice. This correlation was quantitatively probed using the Mantel matrix permutation test94, which takes as a null hypothesis that any two overlapping matrices are uncorrelated. Visualizations for paired scRNA-seq and V(D)J-seq were generated using the Python-based ‘matplotlib’ package.

ATAC–seq

A total of 40,000 cells were sorted per population per replicate for ATAC–seq, with replicates one and two originating from a single mouse each and replicate three representing two pooled mice. ATAC–seq libraries were prepared as previously described, with some modifications95. Cells were pelleted in a fixed rotor benchtop centrifuge at 500g for 5 min at 4 °C. Cells were then washed in 1 ml cold PBS and pelleted again. The supernatant was aspirated and the cells were resuspended in 50 µl ice-cold cell lysis buffer (10 mM Tris-Cl pH 7.4, 10 mM NaCl, 3 mM MgCl2, 0.1% NP-40) to disrupt plasma membranes. Nuclei were pelleted at 1,000g for 10 min. The supernatant was aspirated and the nuclei were resuspended in 40 µl of transposition reaction mixture (Illumina Tagment Kit: 20 µl TD buffer; 2 µl TDE1; 18 µl ddH2O). Samples were incubated in a ThermoMixer at 1,100 RPM for 45 min at 42 °C. DNA was then purified using a MinElute Reaction Cleanup Kit, according to the manufacturer’s instructions. DNA was eluted in 10 µl buffer EB. Libraries were then barcoded and amplified with NEBNext High-Fidelity Master Mix and primers described in a previous publication96 (50 µl reaction with 10 µl DNA and 2.5 µl of 25 µM primers; one cycle of 5 min at 72 °C, 30 s at 98 °C; five cycles of 10 s at 98 °C, 20 s at 63 °C, 1 min at 72 °C). A qPCR analysis on the product determined that an additional seven cycles (10 s at 98 °C, 20 s at 63 °C, 1 min at 72 °C) were required. The library was purified and size-selected with AMPure XP beads: 45 µl of PCR product was incubated with 18 µl beads and the supernatant was collected (beads bound larger than ~2,000 bp fragments). The supernatant (63 µl) was then incubated with an additional 63 µl of beads for 5 min, the supernatant was removed, the beads were washed twice with 75% ethanol and DNA was eluted into 50 µl H2O by incubating for 2 min. Samples were quality-control-checked and quantified on an Agilent BioAnalyzer by IGO. Paired-end 50 bp reads, 20–30 million per sample, were sequenced on an Illumina HiSeq 3000 by IGO.

ATAC–seq data processing

Samples were processed and aligned using Trimmomatic (v.0.39), STAR (v.2.7.3a) and Samtools (v.1.12) with the following steps, where Sample represents each Il10neg, Il10recent or Il10stable replicate:

TrimmomaticPE Sample_R1.fastq.gz Sample_R2.fastq.gz -baseout Sample.fastq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

STAR–runThreadN 6–runMode alignReads–genomeLoad NoSharedMemory–readFilesCommand zcat–genomeDir mm39_100–readFilesIn Sample_1P.fastq.gz Sample_2P.fastq.gz–outFileNamePrefix Sample–outSAMtype BAM Unsorted–outBAMcompression 6–outFilterMultimapNmax 1–outFilterMismatchNoverLmax 0.06–outFilterMatchNminOverLread 0.35–outFilterMatchNmin 30–alignIntronMax 1–alignEndsType Local

samtools sort -@ 4 -n -o Sample.bam SampleAligned.out.bam

samtools fixmate -@ 4 -rm Sample.bam Sample.fixmate.bam

samtools sort -@ 4 -o Sample.resort.bam Sample.fixmate.bam

samtools markdup -@ 4 -l 1500 -r -d 100 -s Sample.resort.bam Sample.duprm.bam

samtools sort -@ 4 -n -o Sample.byname.bam Sample.duprm.bam

samtools index -@ 4 -b Sample.duprm.bam

This procedure resulted in the retention of all uniquely aligning reads, with PCR and optical duplicates removed, to be used for downstream analysis. Then, peaks were called across the three replicates of each cell population individually using Genrich (v.0.5), with the following command (in which Celltype stands in for each cell population, Sample (r1–r3) represents the three replicates and X is 0.002% of the mean number of uniquely aligned reads for each cell population)97:

Genrich -t Sample_r1.byname.bam,Sample_r2.byname.bam,Sample_r3.byname.bam -o./Gen_out/Peak/Celltype.narrowPeak -j -d 25 -g 5 -v -q 0.01 -a X.

Peak atlases for each population were then concatenated, sorted and clustered using bedtools (v.2.27.1) to identify overlapping peaks with the following command90:

bedtools cluster -d -1 -i combined_sort.narrowPeak > clustered.narrowPeak

A custom R script was then used to merge the atlases according to the following principles. If all the peaks in a cluster entirely overlapped, defined as all peak summits falling within the maximal start and minimal end positions of the cluster, the merged peak was defined as the mean start, summit and end of all peaks in that cluster. This was the case for ~82% of all peaks. In the other cases, clusters had multiple distinct summits. These clusters were divided into distinct peaks, with one for each distinct summit, and the boundaries were defined by the most proximal downstream and upstream start and end positions within the cluster. Peaks assigned to regions of the assembly not corresponding to any chromosome and peaks with a width >3,500 bases were eliminated. The getfasta function of bedtools and a custom R script were used to identify and eliminate peaks with >70% repetitive elements. The final combined atlas contained 70,323 peaks. The R packages ‘GenomicRanges’ and ‘ChIPpeakAnno’ were used to assign peaks to the closest gene according to the following principles81,98,99. Peaks 2,000 bases upstream or 500 bases downstream of a transcription start site were considered ‘promoter’ peaks. Non-promoter peaks within the body of a gene were considered ‘intragenic’ peaks. Peaks 100,000 bases up or downstream of a gene body were considered ‘intergenic’. All other peaks were not assigned to any specific gene. Differential accessibility analysis was carried out using the ‘DESeq2’ package, with the formula ‘~ Celltype + Replicate’, where Celltype is either Il10neg, Il10recent or Il10stable and replicates are the separate samples from which each of the three populations were sorted. Differential accessibility analysis and statistical testing were performed for all pairwise comparisons of ‘Celltype’: Il10neg, Il10recent and Il10stable.

Motif identification and model generation

Motif discovery in peaks and model creation were carried out as previously described44, with some modifications. Motifs for all mouse TFs were downloaded from CisBP (v.2.00)100,101. TFs with mean FPKM > 1 in any cell population were used for further analysis. This resulted in 319 motifs for 188 TF-encoding genes. The ‘AME’ software from the MEME suite (v.5.3.0) was used to identify which of these motifs were enriched in the sequences corresponding to the combined peak atlas102,103. At this stage, the best motif for each TF-encoding gene, defined as being detected in the highest fraction of peaks, was chosen for further analysis. Then, the ‘Tomtom’ software from the MEME suite was used to determine closely related motifs104. A custom R script was used to group motifs according to the following principles. Motifs were grouped if their ‘Tomtom’ E-value was <0.00001 and if they were in the same protein family (according to CisBP). Then, groups were merged based on overlapping members until each gene belonged to at most one group. This resulted in 142 motif groups. Finally, motifs not significantly enriched in the peak atlas (‘AME’-adjusted P < 0.01) were excluded. This resulted in 76 TF-encoding genes organized into 58 groups. The ‘FIMO’ software from the MEME suite was used to identify individual instances (P < 0.0001) of each of the 76 motifs across the entire peak atlas105. Motif families occurring in fewer than 2% of peaks were eliminated. This resulted in a final set of 57 motifs within 40 groups. Finally, a peak-by-motif matrix was generated, in which 1 indicated at least one instance of a motif belonging to that family and 0 indicated no motif.

The ‘ridge’ package in R was used to fit a linear ridge regression for the log2FC in accessibility at each peak as a function of the peak-by-motif matrix45,46. This package applies an algorithm for semi-automatically determining the optimal ridge parameter(s) to use to maximize model performance and also performs significance testing using the method of Cule45. This was done separately for the Il10stable vs Il10recent log2FC (the svr model) and Il10stable vs Il10neg log2FC (the svn model). Motif families with significant coefficients in the models (P < 0.001) were used in subsequent analyses. At the same time, linear ridge regressions were fit as above, except with each motif individually removed from the matrix, generating a series of ‘zeroed-out’ models. Then, for sets of peaks of interest (for example, those associated with a specific cluster of genes), the correlations between the actual log2FC and the log2FC predicted by the svr or svn model as well as the correlations between the actual log2FC and the log2FC predicted by the ‘zeroed-out’ svr or svn models were determined. Decreased correlation for the ‘zeroed-out’ model was assumed to be indicative of the ‘zeroed-out’ motif disproportionately contributing to the model’s predictiveness at those specific peaks, and therefore potentially regulating accessibility.

Plotting RNA-seq and ATAC–seq tracks

The UCSC utility ‘faCount’ was used to determine the effective genome size of the ‘reformed mm39’ genome106. The bamCoverage function of deeptools (v.3.5.1) was used to generate ‘bigwig’ files for ATAC–seq and RNA-seq tracks, with the following command107:

bamCoverage -b Sample.duprm.bam–effectiveGenomeSize x -bs 1–maxFragmentLength y–scaleFactor z -o Sample.normdt.bw, where Sample represents each Il10neg, Il10recent or Il10stable replicate; x is the effective genome size as defined above; y is the equivalent value determined by STAR during alignment ((2ˆwinBinNbits)*winAnchorDistNbins) for ATAC–seq or the default for RNA-seq; and z is the inverse of the size factors determined by DESeq2.

TCR deletion RNA-seq

A total of 5,000 cells were sorted per population per replicate for bulk RNA-seq. RNA was extracted and libraries were prepared using SMARTer Stranded RNA-Seq Kits according to the manufac