Mice<\/p>\n
The mice used in this study were housed in pathogen-free facilities at the University of Chicago and Stanford University. All mice were housed in positively pressurized, individually ventilated cage racks and changed in biological safety cabinets. Cage supplies were sanitized using hot water (82\u2009\u00b0C). Bedding and shredded-paper enrichment were autoclaved and cages were provided with irradiated food. Reverse Osmosis water was provided by an automated watering system directly to each cage. Rodent housing rooms were maintained at a 12\u2009h:12\u2009h light:dark cycle. Temperature and humidity were within the Guide for the Care and Use of Laboratory Animals recommended ranges: 20\u201326\u2009\u00b0C and 30\u201370% humidity. All experiments and animal-use procedures were conducted in compliance with the Guide for the Care and Use of Laboratory Animals and were approved by the Institutional Animal Care and Use Committee (IACUC) at the University of Chicago. B6.129-Trp53LSL-L25Q,W26S,F53Q,F54S heterozygous mice27<\/a>,61<\/a> were provided by Laura Attardi (Stanford University) and were bred with B6-Foxn1cre homozygous mice62<\/a> purchased from Jackson Laboratories to generate Trp53LSL-L25Q,W26S,F53Q,F54S\/wt;Foxn1cre\/wt and Trp53wt\/wt;Foxn1cre\/wt littermates. Trp53fl\/fl mice were purchased from Jackson Laboratories and bred with B6-Foxn1cre mice to generate Trp53fl\/fl;Foxn1cre\/wt mice. C57BL\/6J mice were purchased from Jackson Laboratories. mTECs and thymocytes were collected from mice 4\u20135 weeks old. Sex-matched littermates were used for all comparisons of genetic perturbations.<\/p>\n Isolation, sorting and analysis of mouse mTECs<\/p>\n Thymic epithelial cells were isolated as previously described63<\/a> with minor modifications. In brief, thymi from 4\u20136-week-old mice were removed and connective tissue was removed. Stromal tissue was perforated using scissors and incubated with rotation in DMEM-F12 (Gibco) at room temperature for 10\u2009min to liberate the thymocytes. The remaining stromal tissue was enzymatically digested (0.5\u2009mg\u2009ml\u22121 Collagenase D (MilliporeSigma), 0.2\u2009mg\u2009ml\u22121 DNaseI (MilliporeSigma), 0.5\u2009mg\u2009ml\u22121 Papain (Worthington Biochemical)). Cells were stained with anti-EpCAM antibodies conjugated to APC-Cy7 (clone G8.8, BioLegend, 3\u2009\u00b5l per 100 million cells) and EpCAM+ cells were enriched by positive selection using magnetic anti-Cy7 beads (Miltenyi, 10\u2009\u00b5l per 100 million cells). The enriched fraction was stained with the appropriate panel of fluorochrome-conjugated antibodies to CD45 (clone 30-F11, Invitrogen, 1:100), Ly-51 (clone 6C3, BioLegend, 1:100), MHC-II I-A\/I-E (clone M5\/114.15.2, Invitrogen, 1:100), CD104 (clone 346-11A, BD Biosciences, 1:200), GP2 (clone 2F11-C3, MBL, 1:10), CD177 (clone 1171\u2009A, R&D, 1:25), Ly-6D (clone 49-H4, Invitrogen, 1:200), Sca-1 (clone D7, BioLegend, 1:200), AIRE (clone 5H12, Invitrogen, 1:500), Ki-67 (clone SolA15, Invitrogen, 1:100), SynCAM (clone 3E1, MBL, 1:100), CD171\/L1CAM (clone 555, Miltenyi, 1:25) along with fluorescein-labelled UEA-I (Vector Labs, 1:100), Zombie Aqua (BioLegend, 1:500) and DAPI (Invitrogen, 1:20). Intracellular staining for AIRE and Ki-67 was subsequently done using the eBioscience FoxP3 transcription factor staining kit (Invitrogen) according to the manufacturer\u2019s instructions. Intracellular staining for MDM2 (clone EPR22256-98, Abcam, 1:25) was also done using the eBioscience FoxP3 transcription factor staining kit (Invitrogen) according to the manufacturer\u2019s instructions with the addition of a 1-h incubation in blocking buffer (eBioscience permeabilization buffer with 5% normal donkey serum) before a secondary stain (BV412 donkey anti-rabbit, Jackson Immuno, 1:50). Cells were sorted using FACS Symphony S6, FACSAria Fusion or FACSAria II equipped with a 100-\u03bcm nozzle (BD Biosciences). Flow-cytometry data for thymic mimetic cells were acquired using a Cytek Aurora. All other flow-cytometry data were acquired using a BD LSRII or Fortessa. All flow-cytometry data were analysed using FlowJo (v.10).<\/p>\n Human thymic tissue acquisition and processing<\/p>\n Thymus fragments were obtained from a 12-week-old human patient with no known genetic abnormalities undergoing standard-of-care cardiac surgery. The patient was de-identified on receipt with written informed consent for the release of genomic sequence data in accordance with IRB protocol 20\u20131392 approved by the Biological Sciences Division and University of Chicago Medical Center Institutional Review Boards at the University of Chicago and protocol 2020-203 approved by the Advocate Aurora Health Research Subject Protection Program and Advocate Aurora Health Care Institutional Review Board. Connective tissue was removed and the remaining tissue was minced, then incubated with rotation in DMEM-F12 (Gibco) at 4\u2009\u00b0C for 20\u2009min to liberate the thymocytes. Stromal tissue was enzymatically digested using 0.5\u2009mg\u2009ml\u22121 Collagenase D (MilliporeSigma) and 0.2\u2009mg\u2009ml\u22121 DNase I (MilliporeSigma) at 37\u2009\u00b0C for 20\u2009min. The remaining fragments were incubated with rotation in 0.5\u2009mg\u2009ml\u22121 Papain (Worthington), 0.25\u2009mg\u2009ml\u22121 Collagenase D and 0.1\u2009mg\u2009ml\u22121 DNase I at 37\u2009\u00b0C for 20\u2009min. Cells were stained with anti-EpCAM antibodies conjugated to APC-Cy7 (clone 9C4, BioLegend, 1:100) and EpCAM+ cells were enriched by positive selection with magnetic anti-Cy7 beads (Miltenyi). The enriched fraction was stained with DAPI (Invitrogen, 1:20), CD45 (clone 2D1, BioLegend, 1:100), LY51\/CD249 (clone 2D3\/APA, BD Biosciences, 1:00) and HLA-DRA (clone L243, BioLegend, 1:100) and sorted on a Symphony S6 (BD Biosciences).<\/p>\n Flow cytometry of thymocytes and splenocytes<\/p>\n Thymi from 4\u20136-week-old mice were removed and small cortical incisions were made before mechanical agitation with wide-bore glass pipettes in DMEM\/F-12 (Gibco) to liberate the thymocytes. Spleens from mice aged 4\u2009weeks to 12\u2009months old were isolated in RPMI (Gibco) supplemented with 10% FCS. Cells were liberated by mincing with a syringe plunger and filtered through a 40-\u03bcm strainer. Following red blood cell lysis (BD PharmLyse), cells were stained with fluorochrome-conjugated antibodies specific for mouse CD4 (GK1.5, 1:100), CD8\u03b1 (53-6.7, 1:100), CD25 (PC61, 1:100), CD44 (IM7, 1:100), CD69 (H1.2F3, 1:100), CD62L (MEL-14, 1:100), TCR\u03b2 (H57-597, 1:100) and DAPI (Invitrogen, 1:20). Intracellular staining for FoxP3 (clone FJK-16s, eBioscience, 1:100) was done using an eBioscience FoxP3 transcription factor staining kit (Invitrogen) according to the manufacturer\u2019s instructions. Flow-cytometry data were acquired using a BD LSRII or Fortessa and analysed using FlowJo (v.10).<\/p>\n Bulk RNA-seq sample preparation<\/p>\n We FACS-sorted 75,000 primary mTECs directly into RULT lysis buffer (Qiagen RNEasy UCP Micro Kit) and total RNA was extracted following the manufacturer\u2019s instructions. The mRNA was enriched and RNA-seq libraries were constructed using an Illumina TruSeq Stranded mRNA kit. Paired-end, dual-index sequencing was performed on an Illumina NovaSeq 6000 platform.<\/p>\n Bulk RNA-seq data processing<\/p>\n RNA-seq reads were mapped to the mm10 mouse genome assembly using TopHat (v.2.1.1) with the setting \u2013microexon-search. Unmapped, unpaired and low-quality reads (MAPQ\u2009\u2264\u20095) were removed using samtools (v.1.9) view with settings -q 5 -f 2. Paired reads were counted for each gene using featureCounts from Subread (v.2.0.1). TPM values were calculated for each gene to quantify the relative abundance of transcripts for clustering analysis. The trimmed mean of M values was calculated for each gene for differential comparisons across samples using edgeR (v.4.0.2) (calcNormFactors()). Common dispersions were estimated using estimateCommonDisp() and Benjamini\u2013Hochberg FDRs were calculated for pairwise comparisons using the exactTest(). Genes with FDR\u2009\u2264\u20090.05 were regarded as significant.<\/p>\n Definition of tissue-specific and AIRE-dependent genes<\/p>\n Previously published transcriptional data64<\/a> from Aire wild-type and Aire-knockout mTEChi were analysed according to the bulk RNA-seq pipeline outlined above. Genes that exhibited at least 1.5-fold induction in Aire wild type relative to Aire knockout and had Benjamini\u2013Hochberg FDR\u2009\u2264\u20090.05 were regarded as Aire-induced. TSGs were classified as previously64<\/a>, and \u03b1TSGs were taken to be the intersection of these two gene sets. For human TSGs, GTEx65<\/a> expression counts (median TPM), Shannon entropy \\(\\left(S=-\\sum p{\\log }_{2}p\\right)\\) across tissues was calculated for each gene. Genes with an entropy S\u2009\u2264\u20093 were included for downstream analyses.<\/p>\n Multiome sample preparation and sequencing<\/p>\n For all Multiome experiments, we used an ATAC\u2009+\u2009GEX single-cell kit and protocol (10X Genomics 1000236 with protocol CG000338<\/a> RevE) with minor modifications to sample preparation. In brief, 40,000 mTECs were FACS-sorted into 1\u00d7 PBS supplemented with 2% BSA and centrifuged at 300g for 5\u2009min. Cells were gently washed in 50\u2009\u03bcl lysis buffer (10\u2009mM Tris, 10\u2009mM NaCl, 3\u2009mM MgCl2 in nuclease-free water) and centrifuged at 300g for 5\u2009min. Cells were resuspended in 50\u2009\u03bcl permeabilization buffer (10\u2009mM Tris, 10\u2009mM NaCl, 3\u2009mM MgCl2, 0.1% Tween20, 0.01% digitonin and RNase inhibitor (Invitrogen) in nuclease-free water) and incubated for 5\u2009min on ice. Nuclei were gently washed with wash buffer (10\u2009mM Tris, 10\u2009mM NaCl, 3\u2009mM MgCl2, 0.1% Tween20 and RNase inhibitor in nuclease-free water) and centrifuged at 500g for 5\u2009min. Finally, nuclei were resuspended in 5\u2009\u03bcl chilled diluted nuclei buffer (10X Genomics) and added to the transposition mix. Paired-end, dual-index sequencing was performed on an Illumina NovaSeq 6000 platform.<\/p>\n Multiome data quality control<\/p>\n After sequencing, bcl files were converted to fastq using cellranger-arc (v.2.0.2) mkfastq. FASTQ files were aligned to the mm10 or hg38 genome assembly using cellranger-arc count. ATAC-seq fragment files were used as inputs to the ArchR66<\/a> (v.1.0.2) analysis pipeline in R (v.4.3.2). Transcript count matrices were used as inputs to the Seurat (v.5.1.0) gene expression analysis pipeline. For gene expression quality control, cells with nFeature_RNA\u2009\u2265\u2009250 and \u2264\u20096,000, nCount_RNA\u2009\u2264\u200925,000 and percent_mitochondrial \u2264\u200925 were included for downstream analyses. Transcript counts were log-normalized. For scATAC-seq quality control, cells with n_ATAC_Frags\u2009\u2265\u20093,000 and TSS_Score\u2009\u2265\u200910 were included for downstream analyses. Doublet inference was conducted using ArchR addDoubletScores(), and presumed doublets were excluded. Cells that passed each filter were admitted for downstream analyses. Finally, based on gene expression markers, contaminating cells (thymocytes) and putative mTEC mimetic cells were excluded from analysis (except for targeted analyses of mimetic compartments). In the wild-type multiome (Fig. 1<\/a>), a further cluster of cells that exhibited uncharacteristically low TSS enrichment scores was excluded.<\/p>\n Multiome data processing<\/p>\n Dimensionality reduction, scATAC-seq clustering, projections, pseudotime, transcription factor motif enrichment (except for scATAC-seq fragments or genomic tiles, which was computed using HOMER2 (v.5.1) findMotifsGenome.pl with settings -size given), and transcription factor footprinting were performed using the ArchR pipeline with default parameters. For UMAP plots overlaid with continuous colour scales, MAGIC67<\/a> (v.2.0.3) imputation was used for data smoothing to facilitate better visualization. MAGIC-imputed values were used for UMAP display purposes only; imputed values were not used anywhere else in the analysis of scATAC-seq or scRNA-seq datasets (such as violin plots or heatmaps). For scATAC-seq peak calling, the standard ArchR workflow was used using MACS2 (v.2.2.9.1). To maximize the detection of open chromatin regions specific to each sample and stage in the mTEC developmental trajectory, fixed-width 501-bp scATAC-seq peaks were called (extendSummits\u2009=\u2009250) on the Tn5-corrected single base insertions (shift\u2009=\u2009\u221275, extsize\u2009=\u2009150, \u2013nomodel) for each scATAC-seq cluster identified per sample (groupBy\u2009=\u2009Clusters, reproducibility\u2009=\u20091) using the ArchR wrapper function addReproduciblePeakSet(). The significance of each called peak was calculated as a false discovery rate (q-value) comparing the observed number of Tn5 insertions in the sliding window (300\u2009bp) and the expected number of insertions (total number of insertions\/genome size (\u2013nolambda)). A q-value cutoff (cutOff\u2009=\u20090.1) and an upper limit for the number of peaks called per cell (peaksPerCell\u2009=\u20091,000, minCells\u2009=\u2009100) were applied to prevent consideration of low-quality peaks. We also excluded peaks that mapped to the mitochondrial or Y chromosomes (excludeChr\u2009=\u2009c(chrM, chrY)). Peak sets called from each scATAC-seq cluster from respective samples were combined and trimmed for overlap using an iterative procedure that discarded any peak that directly overlapped with the most significant peak66<\/a>. The resultant \u2018union peak set\u2019 was applied to all cells for WIP and OOP count-based and motif-based analyses. The fraction of fragments within peaks was computed automatically as a product of the addReproduciblePeakSet() function. Subnucleosomal and mononucleosomal fractions for each cell or sample were computed as the fraction of the cell\u2019s scATAC-seq fragments whose length L\u2009\u2264\u2009100\u2009bp (subnucleosomal) or 100\u2009<\u2009L\u2009\u2264\u2009200\u2009bp (mononucleosomal). To ensure reproducibility of bioinformatic analysis results, for each dataset, a single script was used for all the quality control and pre-processing, including purging of low-quality cells, doublet removal, peak calling, motif enrichment, dimensionality reduction and clustering. A file representing the full processed data was saved using saveArchRProject() and loaded for all subsequent analyses (this file was not edited after pre-processing). More individual scripts were used to load processed data and perform specific analyses or generate specific figures.<\/p>\n Peak-centric differential accessibility analysis<\/p>\n Differential chromatin accessibility analysis across peaks was done using ArchR getMarkerFeatures() with the following arguments: useMatrix\u2009=\u2009PeakMatrix, bias\u2009=\u2009c(TSSEnrichment, log10(number of scATAC-seq fragments)), testMethod\u2009=\u2009wilcoxon.<\/p>\n Processing of OOP scATAC-seq fragments<\/p>\n For each Multiome dataset, WIP and OOP fragments near genes of interest (such as \u03b1TSGs, housekeeping genes and maturation-induced genes) were retrieved using the ArchR and GenomicRanges R packages. For each gene: first, a search window, search_window, was established around the \\({\\rm{TSS}}({\\rm{search}}\\_{\\rm{window}}={\\rm{TSS}}\\pm {\\ell })\\); and second, scATAC-seq fragments intersecting the search_window were retrieved from cells of interest, cell_subset, using the ArchR getFragmentsFromProject() function with arguments subsetBy\u2009=\u2009search_window and cellNames\u2009=\u2009cell_subset. Fragments were then partitioned based on whether they overlapped the data\u2019s union peak set using subsetByOverlaps() with arguments invert\u2009=\u2009FALSE to retrieve WIP fragments, or invert\u2009=\u2009TRUE to retrieve OOP fragments. Finally, fragments were binned and\/or tallied for the specific application (see below).<\/p>\n Analyses comparing \u03b1TSGpos and \u03b1TSGneg mTECs<\/p>\n Cells from early mature, mid mature and late mature clusters expressing any \u03b1TSGi\u2009>\u20090 were selected as the \u03b1TSGpos cohort and a size-matched cohort of \u03b1TSGneg cells was sampled randomly from the remaining cells from the same three clusters. These cohorts were then used as inputs to getMarkerFeatures()in ArchR for differential accessibility of peaks between \u03b1TSGpos and \u03b1TSGneg mTECs. For local OOP and WIP analysis, ATAC-seq fragments within peaks and outside of peaks from \u03b1TSGpos and \u03b1TSGneg cohorts were intersected with a \u00b15\u2009kb sliding window with 1\u2009kb increments, normalized to the total number of ATAC-seq fragments per cell, and tallied in each window within a region flanking \u03b1TSGi\u2009. For \u03b1TSG coexpression analysis, the probability of detecting each \u03b1TSGi\u2009 neighbouring \u03b1TSG0\u2009 within the specified length scale (or a randomly selected alternative \u03b1TSG as a control) was computed for each of the \u03b1TSGpos and \u03b1TSGneg cohorts.<\/p>\n Regression analysis<\/p>\n For each \u03b1TSGi, the total number of OOP and WIP scATAC-seq fragments within the characteristic window of instability \\(({\\ell }=\\pm 50\\,{\\rm{kb}})\\) was computed for each mTEC in the early mature, mid mature and late mature clusters. A logistic regression framework was used (glm() with family\u2009=\u2009binomial) to estimate the probability of expressing a given \u03b1TSG based on the number of log10(OOP\u2009+\u20091) or log10(WIP\u2009+\u20091) fragments using log10(n_ATAC_Frags) per cell as a covariate. P-values for regression coefficients were generated using the Wald-\u03c72 test (anova(test\u2009=\u2009\u2018LR\u2019)).<\/p>\n CUT&RUN sample preparation<\/p>\n CUT&RUN was performed as previously described28<\/a> with minor modifications. In brief, 350,000\u2013500,000 cells were washed 3 times in wash buffer (20\u2009mM HEPES pH 7.5, 150\u2009mM NaCl, 0.5\u2009mM spermidine, 1\u00d7 EDTA-free protease inhibitor cocktail (Roche)) then bound to Concanavalin-A beads (Bangs Laboratories) according to the manufacturer\u2019s instructions. Cells were incubated with 1:100 dilution of anti-p53 antibody (Leica NCL-L-p53-CM5p) for 2\u2009h or overnight at 4\u2009\u00b0C in permeabilization buffer (1\u00d7 permeabilization buffer (eBioscience), 0.5\u2009mM spermidine, 1\u00d7 EDTA-free protease inhibitor cocktail, 2\u2009mM EDTA). The sample was then incubated with 700\u2009ng\u2009ml\u22121 pA-MNase (S. Henikoff) in permeabilization buffer at 4\u2009\u00b0C for 1\u2009h. Digestion was done in 0.5\u00d7 permeabilization buffer supplemented with 2\u2009mM CaCl2 at 4\u2009\u00b0C for 1\u2009h. The reaction was stopped by the addition of 2\u00d7 stop buffer (final concentration 100\u2009mM NaCl, 10\u2009mM EDTA, 2\u2009mM EGTA, 20\u2009\u03bcg\u2009ml\u22121 glycogen, 25\u2009\u03bcg\u2009ml\u22121 RNase A (Thermo Fisher)) and the sample was incubated at 37\u2009\u00b0C for 20\u2009min. Protein in the sample was then digested in 0.1% SDS and 250\u2009\u03bcg\u2009ml\u22121 Proteinase K (New England Biolabs) for 2\u2009h at 56\u2009\u00b0C, shaking gently. CUT&RUN fragments were purified by phenol chloroform extraction. CUT&RUN libraries were generated using NEBNext UltraII DNA Library Prep Kit for Illumina coupled with NEBNext Multiplex Oligos for Illumina (New England Biolabs) with modifications optimized for small fragments, as detailed in https:\/\/doi.org\/10.17504\/protocols.io.wvgfe3w<\/a>. Paired-end, dual-index sequencing was performed on the Illumina NextSeq500 platform.<\/p>\n CUT&RUN data processing<\/p>\n CUT&RUN reads were mapped to mm10 mouse genome assembly using Bowtie2 (v.2.2.9) with settings –local –very-sensitive-local \u2013no-unal \u2013no-mixed \u2013no-discordant \u2013phred33 -I 10 -X 700. PCR duplicates were removed using Picard (v.2.21.8) MarkDuplicates REMOVE_DUPLICATES=true VALIDATION_STRINGENCY\u2009=\u2009LENIENT. Reads with MAPQ scores below 30 were purged and excluded from downstream analysis using samtools (v.1.9) view -b -q 30 -f 2 -F 1804. Peaks were called for each sample using MACS2 (v.2.2.7.1) with settings –shift 0 –extsize 200 –nomodel –call-summits –keep-dup all -p 0.01. For each sample, a 301-bp fixed-width peak set was generated by extending the MACS2 summits by 150\u2009bp in both directions. Peaks were ranked by significance (MACS2 peak score) and overlapping peaks with lower peak scores were removed iteratively to create non-overlapping sample peak sets. Peaks mapping to chrY, as well as any that spanned genomic regions containing \u201cN\u201d nucleotides, were removed. Robust peaks were defined by a score per million (SPM) (each peak score divided by the sum of all peak scores in the sample, divided by 1 million), and we retained only those peaks with SPM\u2009\u2265\u20095. We defined p53 CUT&RUN peaks by further filtering for peaks that overlapped with known p53-binding motifs (HOMER2, v5.1) from samples with characterized p53 activity (mTEClo samples). CUT&RUN fragment counts across regions of interest were normalized by the number of unique fragments in the sample library.<\/p>\n ChIP\u2013seq data processing<\/p>\n ChIP\u2013seq reads were mapped to mm10 mouse genome assembly using Bowtie2 (v.2.2.9) with settings –very-sensitive -X 2000. PCR duplicates were removed using Picard (v.2.21.8) MarkDuplicates REMOVE_DUPLICATES=true VALIDATION_STRINGENCY\u2009=\u2009LENIENT. Reads with MAPQ scores below 30 were purged and excluded from downstream analysis using samtools (v.1.9) view -b -q 30 -F 1796. ChIP\u2013seq read counts were normalized by the number of unique reads in the sample library.<\/p>\n Histopathology<\/p>\n