Skin-resident CD8+ T cells comprise distinct IFN-γ- (TRM1) and IL-17-producing (TRM17) subsets that differentially contribute to immune responses. However, whether these populations employ common mechanisms to establish tissue residence is unknown. Here, we show that TRM1 and TRM17 cells navigate divergent trajectories to acquire tissue residency in skin. While TRM1 cells depend on a T-bet-Hobit-IL-15 axis, TRM17 cells develop independently of these factors. Instead, c-Maf commands a tissue-resident program in TRM17 cells parallel to that induced by Hobit in TRM1 cells, with an ICOS-c-Maf-IL-7 axis pivotal to TRM17 cell commitment. Accordingly, targeting this pathway enables ablation of skin TRM17 cells without compromising their TRM1 counterparts. Thus, skin-resident T cells rely on distinct molecular circuitries, which can be exploited to strategically modulate local immunity.

Mice

C57BL/6, B6.SJL‐PtprcaPep3b/BoyJ (CD45.1), B6.SJL‐PtprcaPep3b/BoyJ×C57BL/6 (CD45.1×CD45.2), OT-I×B6.SJL-PtprcaPep3b/BoyJ (OT-I.CD45.1), T-bet-ZsGreen×Rorgt-E2Crimson, Tgfbr2^flox/flox×dLck-cre (Tgfbr2^–/–), Tbx21^–/–, Tbx21^–/–×Eomes×Cd4^cre,  Nfil3^–/–, Hobit^Tomato, Blimp1^Tomato, Il15^–/–, lox-stop-lox-Tgfbr1^CA (Tgfbr1^CA)×Cd8(E8i)^cre, Trac^flox/flox×Rosa26^creERT2/+, Zfp683^–/– (Hobit-KO), Maf^flox/flox×Cd4^cre (c-Maf-cKO) and Icos^–/– mice were bred and maintained in the Department of Microbiology and Immunology at the The University of Melbourne, Australia. Hobit^Tomato and Blimp1^Tomato were generated by the Kallies laboratory by inserting tandem tomato expression cassettes into the endogenous Zfp683 and Blimp1 loci, respectively. For parabiosis experiments, C57BL/6 and B6.SJL‐PtprcaPep3b/BoyJ (CD45.1) mice were obtained from either the NIAID Taconic exchange program or purchased from The Jackson Laboratory and maintained at the National Institutes of Health, Bethesda, USA. Bowie×Rag1^−/− transgenic mice (6) were bred and maintained at the National Institutes of Health, Bethesda, USA. Trac^flox/flox mice were kindly provided by D. Mucida (The Rockefeller University, New York, USA). Tcf7^flox/flox×dLck^cre and Tcf7^flox/flox×Cd4^cre (Tcf7 cKO) mice were kindly provided by L. Mielke (Olivia Newton-John Cancer Research Institute) and were bred and maintained at La Trobe University, Bundoora, Victoria, Australia. For wound-healing assays, male Maf^flox/flox×Cd4^cre mice were bred and maintained at the National Institutes of Health, Bethesda, USA. All other mice used were female and all mice were aged between 6–18 weeks at the beginning of experiments. All animal experiments were approved by the University of Melbourne Animal Ethics Committee or performed in an American Association for the Accreditation of Laboratory Animal Care (AAALAC)-accredited animal facility at the NIAID and housed in accordance with the procedures outlined in the Guide for the Care and Use of Laboratory Animals.

Mouse and human tissue processing

Lymphocytes were isolated from mouse spleens by dissociating tissues through a metal mesh to create single-cell suspensions and then treating them with 1X Red Blood Cell Lysis Buffer (eBioscience) for 2–3 min. Except where indicated, T cells were isolated from mouse flank skin by incubating skin (3 cm²) in Liberase TL Research Grade solution (Sigma, 0.25 mg/ml) in Hank’s balanced saline solution (HBSS) at 37°C for 25 min. The epidermis was then separated from the dermis before finely chopping and incubating for a further 60 min at 37°C. In certain experiments, T cells were isolated from both ear pinnae of mice by separating the dorsal and ventral sides of the ear, incubating in Liberase TL Research Grade solution (Sigma, 0.25 mg/ml in HBSS) for 90 min at 37°C and 6.5% CO₂ before chopping and mixing in RP-10 solution. Skin T_RM cell enumeration is combined from two ear pinnae and 3 cm² of associated flank skin. T cells were isolated from human PBMC by separation on a Ficoll (Sigma) gradient and from human skin by incubating skin for 90 min in Dispase solution (Roche, 2.5 mg/ml) at 37°C and, peeling the epidermis from the dermis and then incubating the skin for a further 90 mins in Collagenase Type III (Worthington, 3 mg/ml) at 37°C. Suspensions were filtered through 70-μm meshes prior to stimulation.

Flow cytometric analysis

For single-cell sequencing experiments, cells were sorted from the skin or spleen using a BD FACS Aria III (BD Biosciences) with the following gating strategies: DAPI^–→TCRβ⁺TCRγδ^–→CD44^hi→CD8β⁺CD4^–→CD45.1⁺CD45.2⁺ or CD45.2⁺ (spleen), DAPI^–→TCRβ⁺TCRγδ^–→CD44^hi→CD8β⁺CD4^–→CD69⁺CD103⁺→CD45.1⁺CD45.2⁺ or CD45.2⁺ (skin), or DAPI^–→TCRβ⁺TCRγδ^–→CD44^lo→Vα2⁺CD8β⁺ (naïve gBT-I cells).

Bone marrow (BM) chimeras

Mixed BM chimeras were generated by irradiating recipient mice twice with 550 rads before injecting 0.25-1×10⁷ congenically marked bone marrow cells from WT or KO mice at a 1:1 ratio to CD45.1⁺ recipients i.v. For single-cell-sequencing experiments, irradiated mice were singly transferred Maf^flox/flox×Cd4^cre (c-Maf-cKO) or Hobit^–/– (Hobit-KO) or WT cells. Mice were rested for at least 8 weeks prior to S. epidermidis association. For tamoxifen treatment, mice were administered 2 mg of tamoxifen diluted in sunflower seed oil (both from Sigma) i.p. daily for a total of five injections.

Bulk RNA-sequencing analysis

RNA sequencing datasets were obtained as either summarized raw counts or raw sequence read data files from NCBI Gene Expression Omnibus (GEO) or Bioproject datasets with accession numbers PRJNA41936 (T_RM1 versus T_RM17 cells) (Ref. 1), GSE70813 (Hobit or Blimp1 knockout (KO) versus WT and T_RM cell core signature) (Ref. 2) and GSE120427 (Maf-KO versus WT) (Ref. 3). To generate T_RM1 and T_RM17 gene expression profiles, S. epidermidis Tc17 and Tc1 samples were selected for analysis. To generate Hobit or Blimp1 KO versus WT (i) or T_RM core signatures (ii) (Ref. 2), either WT, Hobit-KO, and Blimp1-KO bone marrow NKT cell samples (i) or HSV skin T_RM cell, LCMV gut/liver T_RM cell, and HSV/LCMV T_EM and T_CM cell (combined and designated T_CIRC) samples (ii) were selected for analysis. For T_RM1 versus T_RM17 comparisons, raw sequence read data was processed by mapping reads to the Mus musculus (GRCm38/mm10) genome using Rsubread (Ref 4.) (v2.2.6) and reads mapping to NCBI RefSeq annotated genes were counted and summarized with featureCounts (Ref 5). For all data sets, genes were annotated according to the NCBI database and genes with retired Entrez IDs or that failed to achieve a count above 10 in all samples in at least one biological group were removed. To generate Hobit or Blimp1 KO versus WT and T_RM cell core signatures, genes designated rRNA were also excluded. T_RM cell core signatures (Ref 2.) were processed further using the imputation strategy previously described (Ref. 6, 7). Each data set was upper-quartile normalized using EDAseq (Ref. 8) and log₂-transformed with an offset of 1, before applying RUVIII (Ref. 9) with biological replicates nominated as technical replicates, mouse housekeeping genes (Ref. 4) nominated as ‘negative control’ genes, and k factors of unwanted variation where: k = 1 (for c-Maf versus WT), k = 5 (for Hobit/Blimp1versus WT) or k = 3 (for T_RM core signature and T_RM1 versus T_RM17). edgeR (v 3.30.3) (Ref. 10) was used to calculate upper quartile normalization factors and fit gene-wise negative binomial generalized linear models with the output from RUVIII as additional covariates where a prior count of 1 was used for Hobit or Blimp1 KO versus WT or T_RM cell core signature generation or default values applied elsewhere. Differentially expressed (DE) genes were called when likelihood ratio tests returned a Benjamini and Hochberg (Ref. 11) adjusted P-value (false discovery rate) <0.05. Standardized log₂-fold-change (logFC) values (empirical Bayes moderated t-statistics) were calculated by computing log₂ counts-per-million values with edgeR and using limma (Ref. 12) (v3.44.3) to fit gene-wise linear models with the output from RUVIII as additional covariates, where the ‘trend’ and ‘robust’ (Ref. 12) options were specified. To generate Hobit or Blimp1KO versus WT signatures, the corresponding main effects contrast for a 2×2 factorial design was used. To generate the T_RM cell core signature, skin T_RM versus T_CIRC cell (HSV) and liver or gut T_RM versus T_CIRC cell (LCMV) samples were compared and genes that were DE up or DE down in all three comparisons defined the core up and down T_RM cell signatures, respectively. Tests of association of expression changes, as displayed in scatter plots, were performed by constructing a two-way contingency table of the number of genes with concordant/discordant logFCs, where genes with small logFCs (between –0.25 and 0.25) in either comparison were ignored and Fisher’s exact test was applied to the resulting table. Gene sets were selected from lists ranked by adjusted P-value. PCA plots and heatmaps were produced from RUVIII normalized data, where heatmaps used gene-wise standardization to produce a Z-score and genes clustered by Pearson correlation. Barcode enrichment plots were generated using limma by producing separate plots for top up genes and top down genes and then stacked into a single plot and P-values were calculated using the camera method (Ref. 14) applied to the top up and down genes separately. Venn diagram visualization was performed in R using the venneuler package (v1.1.3). Where specified, RNA sequencing data sets were processed to impute missing values. For each gene, when a biological group had (post-filtered) non-zero counts for some replicates but zero counts for other replicates, the zeros were regarded as inconsistent, and thus inferred to be missing data.

Network analysis

Protein–protein network interaction analysis was performed by selecting for differentially expressed genes upregulated in either T_RM1 or T_RM17 cells and annotated with either of the Gene Ontology (GO) terms GO:0003700 (DNA-binding TF activity) or GO:0043565 (sequence-specific DNA binding), as well TFs as per Fig. 1D. Interaction networks were generated in string-db (Ref. 15) with a minimum required interaction score of 0.400 then visualized in Cytoscape v3.2.0. TF with less than 2 interactions in the network or log Fold-Change <0.5 and >–0.5 were excluded from the network. PageRank scores for TF in the network were calculated in igraph with edges weighted by STRING network score. Expression score rankings were generated by assigning each TF a ranking based on adjusted P-value (FDR) and a ranking based on logFC then combining both to produce an average ranking.

Single-cell RNA sequencing library preparation

CD69⁺CD103⁺ CD8⁺ T_RM cells and CD8⁺CD44^hi splenic T cells were sorted from bone marrow chimeras (CD45.1⁺) reconstituted with either (i) Hobit^–/– (Hobit-KO, CD45.2⁺) or B6 (WT, CD45.2⁺) cells or (ii) c-Maf^flox/flox×CD4^cre+/− (c-Maf-cKO, CD45.2⁺) or c-Maf^+/+×CD4^cre+/– (WT, CD45.2⁺) cells 14 days after S. epidermidis association. Sorted cells were stained with TotalSeq-C antibodies against CD49a, CCR6 and CD62L and barcoded with individual Hashtags then pooled. A total of 2–3.6×10⁴ cells per experiment were filtered through a Flowmi cell strainer and loaded onto a 10x chromium controller then prepared for sequencing using a Chromium Next GEM Single Cell 5′ kit with feature barcoding and immune receptor mapping (v2, Dual Index) from 10x and VDJ enrichment kit for mouse T cells (v1.1) from 10x. Libraries were generated according to the manufacturer’s instructions. Libraries were profiled on an Agilent Tapestation and quantified using a KAPA Library Quantification Kit (KK4824) before sequencing on an Illumina NextSeq 6000.

Single-cell RNA sequencing analysis

Sequencing data were generated using the cellranger multi and cellranger vdj pipelines in 10x CellRanger (v7.0.0) with reads aligned to the mm10 genome. Data were processed and analysed in R (version 4.2.1) using the Seurat package (v4.1) with the default parameters used unless otherwise noted. Cells were demultiplexed using HTODemux with margin set to 1 and custom positive.quantile values assigned to each experiment (the minimum value for which RidgePlots were bimodal) and negative droplets and doublets removed. Cells with <1000 or >20,000 counts, <250 features and >5% mitochondrial reads were excluded from analysis. Remaining cells were processed using SCTransform v2 with 4000 variable features and regression of mitochondrial read percentage. Integration was performed using reciprocal PCA and excluded TCR genes, ribosomal protein genes, Hsp and Dnaj family genes and a set of genes upregulated during heat treatment (Ref. 16) from the integration anchor features to minimize effects of digestion-induced stress on clustering. Standard processing and normalization steps were then run on the integrated assay (RunPCA, RunUMAP, FindNeighbors) and FindClusters run with a final resolution of 0.3 (for spleen and skin samples) or 0.4 (for skin samples alone). Clusters with low quality (low RNA reads and high mitochondrial reads) that resolved at high resolution and contaminant clusters comprising γδT cells were removed from the object. Gene expression was visualized using SCTransform normalized counts generated prior to integration. Marker genes for each cluster and differentially expressed genes were determined using PrepSCTFindMarkers, FindAllMarkers or FindMarkers on the unintegrated SCT assay with an adjusted p-value cutoff of 0.05, avg_log2FC cutoff of 0.125, min pct=0.02 and no umi recorrection. Volcano plots were generated in ggplot2 or EnhancedVolcano (v1.14.0). For heatmap visualization, average expression was calculated using the function AverageExpression on normalized counts and heatmaps generated using pheatmap (v1.0.12). Venn diagram visualization was performed using the venneuler package (v1.1.3). Gene module score enrichment was performed using AddModuleScore in Seurat. TCR Clonotype analysis was performed using scRepertoire (v1.6.0) using combineTCR to generate clonotypes and remove cells with multiple TCRα or multiple TCRβ chains (removeMulti = TRUE). Individual clones were called using CTstrict (identical CDR3 nucleotide sequences and the same TCRα and TCRβ chains). Frequency barplots were generated using dittoSeq (v1.8.1). Gene ontology analysis was performed in Metascape (Ref. 17) for GO and Reactome pathways using differentially expressed genes with a P-value<0.01 (with ribosomal protein genes excluded) and extracting the top Summary pathways from the Metascape output ranked by P-value.

Statistics and reproducibility

All experiments were performed at least twice with either similar results obtained, and representative plots shown, or data pooled from multiple experiments. Mice of the same group were randomly allocated to groups, without the need for blinding and sample sizes were selected based on prior research to provide informative results and perform statistical testing.

References

1.              J. L. Linehan et al., Non-classical Immunity Controls Microbiota Impact on Skin Immunity and Tissue Repair. Cell 172, 784-796.e718 (2018).

2.              L. K. Mackayet al., Hobit and Blimp1 instruct a universal transcriptional program of tissue residency in lymphocytes.Science352, 459-463 (2016).

3.              M. K. Zuberbuehler et al., The transcription factor c-Maf is essential for the commitment of IL-17-producing γδ T cells. Nat Immunol 20, 73-85 (2019).

4.              Y. Liao, G. K. Smyth, W. Shi, The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. Nucleic Acids Research 41, e108-e108 (2013).

5.              Y. Liao, G. K. Smyth, W. Shi, featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923-930 (2014).

6.              S. N. Christo et al., Discrete tissue microenvironments instruct diversity in resident memory T cell function and plasticity. Nat Immunol 22, 1140-1151 (2021).

7.              R. Fonseca et al., Runx3 drives a CD8(+) T cell tissue residency program that is absent in CD4(+) T cells. Nat Immunol 23, 1236-1245 (2022).

8.              D. Risso, K. Schwartz, G. Sherlock, S. Dudoit, GC-content normalization for RNA-Seq data. BMC Bioinformatics 12, 480 (2011).

9.              R. Molania, J. A. Gagnon-Bartsch, A. Dobrovic, T. P. Speed, A new normalization for Nanostring nCounter gene expression data. Nucleic Acids Res 47, 6073-6083 (2019).

10.           M. D. Robinson, D. J. McCarthy, G. K. Smyth, edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139-140 (2010).

11.           Y. Benjamini, Y. Hochberg, Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B Stat Methodol 57, 459-463 (1995).

12.           M. E. Ritchie et al., limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 43, e47 (2015).

13.           B. Phipson, S. Lee, I. J. Majewski, W. S. Alexander, G. K. Smyth, Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression. Ann Appl Stat 10, 946-963 (2016).

14.           D. Wu, G. K. Smyth, Camera: a competitive gene set test accounting for inter-gene correlation. Nucleic Acids Res 40, e133 (2012).

15.           D. Szklarczyk et al., STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res 47, D607-d613 (2019).

16.           C. H. O'Flanagan et al., Dissociation of solid tumor tissues with cold active protease for single-cell RNA-seq minimizes conserved collagenase-associated stress responses. Genome Biol 20, 210 (2019).

17.           Y. Zhou et al., Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun 10, 1523 (2019).