Single-cell multiomics identifies Tcf1 and Lef1 as key initiators of early thymic progenitor fate

Abstract

Bone marrow-derived multipotent hematopoietic progenitors seed the thymus and generate early thymic progenitors (ETPs). However, the factors governing ETP formation remain poorly defined. Using single-cell RNA sequencing (scRNA-seq) and single-cell assay for transposase-accessible chromatin with sequencing (scATAC-seq), we dissected the heterogeneity of transcriptomic and chromatin accessibility landscapes in murine ETPs. Whereas Tcf1^- ETPs exhibited higher proliferative capacity, Tcf1⁺ ETPs appeared to be immediate, more robust precursors to T lineage-specified early thymocytes. Prethymic ablation of Tcf1 and its homolog Lef1 severely impaired ETP formation in vivo. Whereas ablating Tcf1 alone had limited impact, loss of both Tcf1 and Lef1 impaired transcriptional activation of Notch1 and Notch pathway effector molecules, including Hes1 and Hhex, accompanied by aberrantly induced B cell and myeloid gene programs. Acute deletion of both factors compromised Notch pathway, glycolysis, and T cell gene programs in emergent ETPs ex vivo. Thus, Tcf1 and Lef1 act upstream of the Notch pathway, functioning as prethymic initiators of ETP fate and intrathymic gatekeepers of ETP identity and T lineage potential.

Figure 1.. Single-cell RNA-seq analysis of WT DN1 thymocytes resolves transcriptomic heterogeneity.

A. Uniform manifold approximation and projection (UMAP) plot of scRNA-seq data on sort-purified WT Lin⁻TCRβ⁻CD8⁻CD44^hiCD25⁻ DN1 thymocytes (pooled from two mice), with 18 clusters identified with Seurat and marked in distinct colors. B–E. UMAP plots showing single-cell transcript levels of characteristic genes in All ETPs (B), proliferative ETPs (C), cells with nonT potentials (D), and T-ILC like cells (E). The color scale denotes log-normalized expression values with arbitrary units. F. scRNA-seq UMAP showing functional subsets, as denoted with distinct colors. G. Heat map showing expression of selected genes in DN1 functional subsets as defined in F, with each column corresponding to a single cell and the color scale representing z-normalized transcript levels with arbitrary units. H. Pseudotime analysis of WT DN1 thymocytes using Monocle3, with the black line denoting the trajectory.

Figure 2.. Single-cell ATAC-seq analysis of WT DN1 thymocytes resolves heterogeneity in chromatin accessibility.

A. UMAP plot of scATAC-seq data on sort-purified WT DN1 thymocytes (pooled from three mice), with 19 clusters identified with Seurat and marked in distinct colors. B. UMAP plot showing scATAC-seq functional subsets based on gene activity scores and projection from scRNA-seq functional subsets. C. Heat map of marker ChrAcc sites in each scATAC-seq functional subset. D. Heat map of gene activity score based on marker ChrAcc sites in each scATAC-seq functional subset. E. Heat map of marker ChrAcc sites in regrouped ETP, nonT_potential and T-ILC like subsets using more stringent criteria (≥2-fold difference in one cluster over the other two clusters). The color scales denote z-normalized signal strength (C, E) or gene activity score (D) with arbitrary units. F. Bar graph showing the odds ratio of link between subset marker gene promoters (as determined with scRNA-seq) and subset marker ChrAcc sites (as determined with scATAC-seq). P values were determined with Fisher exact tests over intersection between marker ChrAcc sites and marker gene-linked sites in each functional subset. G. Pseudo-bulk ChrAcc tracks at select gene loci in ETP subset, with red arrows indicating ETP marker ChrAcc sites, green dotted lines denoting gene transcription start sites (TSS), and arches denoting coaccessibility of TSS and ChrAcc sites. Arches in red highlight coaccessibility of TSS and ETP marker ChrAcc sites.

Figure 3.. Thymic DN1 subsets exhibit distinct capacity and kinetics in generating CD25⁺ cells ex vivo.

A. Gating strategy for identifying and sorting of thymic DN1 subsets. B–C. Detection of DN1-DN3 cells in Lin⁻ cells generated from c-Kit⁺Tcf1⁻ and c-Kit⁺Tcf1⁺ thymic subsets seeded on OP9-DL1 monolayers. At the indicated time periods, CD11c⁻, CD11b⁻, CD19⁻ and NK1.1⁻ lymphocytes were analyzed for CD44 and CD25 expression. Representative contour plots (B) were from 2 independent experiments with 4 biological replicates examined. Cumulative data (C) are means ± s.d. in bar graphs. *, p<0.05; **, p<0.01; ***, p < 0.001 by two-tailed Student’s t-test. D. Detection of CD44 and CD25 expression in Lin⁻ cells generated from CD3ε⁺Tcf1⁺ and CD3ε⁺Tcf1⁻ thymic subsets seeded on OP9-DL1 monolayers for 4 days. E, F. Detection of CD11c⁺, CD11b⁺ (E), CD19⁺ or NK1.1⁺ cells (F) generated from the four thymic DN1 subsets seeded on OP9-DL1 monolayers. Representative contour plots were from 2 independent experiments with 4 biological replicates examined. G. Cumulative data on frequency of CD11c⁺, CD11b⁺, and NK1.1⁺ cells generated from the four thymic DN1 subsets after 4 and 6 days in culture. Data are means ± s.d. Statistical significance for multiple-group comparisons was first determined with one-way ANOVA, and Tukey’s test was used as post hoc correction for indicated pairwise comparison. ***, p < 0.001; ns, not statistically significant.

Figure 4.. Pre-thymic ablation of Tcf1 and Lef1 modestly affects hematopoietic progenitors but impairs ETP formation.

A. Detection of LSKs in Lin⁻ BM cells (top) and SLAM HSCs in LSKs (bottom) in Tcf7^−/−Lef1^−/− mice and WT littermates. B. Detection of LSKs in Lin⁻IL-7Rα⁻ BM cells (top) and Flt3⁺ MPPs in LSKs (bottom). C. Detection of Sca1^medc-Kit^med cells in Lin⁻IL-7Rα⁺ BM cells (top) and Flt3⁺ CLPs in Sca1^medc-Kit^med cells (bottom). Contour plots in A–C are representative from 4 independent experiments, and cumulative data in bar graphs are means ± s.d. *, p<0.05; **, p<0.01; ***, p < 0.001; ns, not statistically significant by two-tailed Student’s t-test. D. Total thymocytes numbers. E. Analysis of DN subsets in thymocytes of WT, Tcf7^−/−, and Tcf7^−/−Lef1^−/− mice following the Lin⁻CD8⁻TCRβ⁻ gate strategy to extract DN cells. F. Detection of c-Kit expression in DN1 and DN2 thymocytes of WT, Tcf7^−/−, and Tcf7^−/−Lef1^−/− mice. Representative contour plots are from ≥3 independent experiments, and cumulative data are means ± s.d. Statistical significance for multiple-group comparisons in D and F was first determined with one-way ANOVA, and Tukey’s test was used as post hoc correction for indicated pairwise comparison. *, p < 0.05; **, p < 0.01; ***, p < 0.001.

Figure 5.. Pre-thymic ablation of Tcf1 and Lef1 impairs ETP transcriptional program.

A–B. UMAP plots of scRNA-seq data of WT, Tcf7^−/−, and Tcf7^−/−Lef1^−/− DN1 thymocytes, all combined (A) or individually overlaid on all-cell background (B), as denoted with distinct colors. Cells from three Tcf7^−/− and three Tcf7^−/−Lef1^−/− mice were pooled for the assay. C. scRNA-seq UMAP plot showing cluster annotation of combined WT+KO single cells. D. scRNA-seq UMAP plots showing distribution of functional subsets of WT cells only (left, KO cells shown in grey) and that of WT+KO single cells (right). E. Bar graph showing relative proportion of WT and each type of KO cells in the indicated scRNA-seq-based functional subsets after normalized to counts of high-quality cells. F–G. Violin plots showing transcript levels of ETP-characteristic genes (F) and unique genes repressed by Tcf1 and Lef1 (G). H. Heat map showing the expression of differentially expressed genes between WT and Tcf7^−/− ETP_quiescent cells by the criteria of ≥1.5-fold difference and padj<0.001. Genes in red denotes those showing the same directional changes between WT and Tcf7^−/− ETP_proliferative cells (fig. S4B). The color scale denotes z-normalized transcript levels with arbitrary units.

Figure 6.. Pre-thymic ablation of Tcf1 and Lef1 perturbs ChrAcc landscape of ETPs.

A–B. UMAP plots of scATAC-seq data of WT and Tcf7^−/−Lef1^−/− DN1 thymocytes, all combined (A) or individually overlaid on all-cell background (B), as denoted with distinct colors. Cells from four Tcf7^−/−Lef1^−/− mice were pooled for the assay. C. scATAC-seq UMAP plot showing four functional subsets based on projection from scRNA-seq-defined subsets. Note that cells that were not projected to these four subsets were labeled in grey. D. Bar graph showing relative proportion of WT and Tcf7^−/−Lef1^−/− cells in the indicated scATAC-seq-based functional subsets after normalized to counts of high-quality cells. E. Detection of non-T lineage markers in DN1 thymocytes using an alternative gating strategy, where Lin⁺ cells were not excluded from the TCRβ⁻CD8⁻ DN thymocytes. The resulting DN1 and CD44^medCD25⁻ subsets were further analyzed for B220, CD19, CD11c and Mac1 expression. Representative contour plots are from ≥3 independent experiments, and cumulative data are means ± s.d. Statistical significance for multiple-group comparisons was first determined with one-way ANOVA, and Tukey’s test was used as post hoc correction for indicated pairwise comparison, *, p < 0.05; **, p < 0.01; ***, p < 0.001. F. Bar graph summarizing overlapping rates between Tcf1 binding peaks and differential ChrAcc sites in KO_specific cluster as compared to WT ETPs. G,H. scATAC-seq data presented as pseudo-bulk ATAC-seq tracks resulting from summation of WT ETP and KO_specific single cells at select ETP-characteristic (G) and Tcf1/Lef1-repressed genes (H), along with sequencing tracks of Tcf1 CUT&RUN in OP9-DL1-derived WT ETPs. In G, cyan bars denote ‘more closed’ ChrAcc sites in KO_specific cells, with filled bars highlighting overlap with Tcf1 binding peaks. In H, orange bars denote ‘more open’ ChrAcc sites in KO_specific cells. Green dotted line denotes gene TSS. Arches on the top of tracks denote coaccessibility of TSS and factor-dependent ChrAcc sites, with red ones highlighting coaccessibility with ETP marker ChrAcc sites (defined in Fig. 2G).

Figure 7.. Acute deletion of Tcf1 and Lef1 perturbs transcriptomes of emergent ETPs.

A. Experimental design for evaluating pre-thymic deletion of Tcf1 and Lef1 on ETP formation ex vivo. B. Detection of DN1 and DN2 formation from WT or Tcf7^−/−Lef1^−/− BM MPPs after ex vivo culture on OP9-DL1 monolayers for 7–9 days, with c-Kit expression analyzed on DN1 cells. Representative contour plots are from 2 independent experiments, and cumulative data are means ± s.d. in bar graphs in lower panels. ***, p < 0.001 by two-tailed Student’s t-test. C. Experimental design for acute deletion of Tcf1 and Lef1 in emergent ETPs. D. Detection of emergent ETPs from CreER⁺WT and CreER⁺dKO BM MPPs after 7-day culture on OP9-DL1 monolayers, with 4-OHT treatment during days 3–5. Representative contour plots are from 2 independent experiments, and cumulative data are means ± s.d. in bar graphs in lower panels. *, p<0.05; ***, p < 0.001 by two-tailed Student’s t-test. E. Volcano plot showing DEGs between CreER⁺WT and CreER⁺dKO emergent ETPs sorted on day 7 of culture, by the criteria of ≥1.5 fold changes and FDR<0.05. F, G. Heat maps showing Tcf1/Lef1-activated genes in glycolysis, Notch and TCR signaling pathways (F), and select differentially expressed transcriptional regulators (G). The color scales denote z-normalized transcript levels with arbitrary units. H, I. Enrichment plots of hypoxia-induced (H) and -repressed gene sets (I) in comparison of CreER⁺WT and CreER⁺dKO ETP transcriptomes as determined with GSEA, with top 20 genes in the leading edge shown in heat maps. NES, normalized enrichment score; NOM P value and nominal P values were defined in GSEA.

Figure 8.. Acute deletion of Tcf1 and Lef1 perturbs ChrAcc landscape in emergent ETPs.

A. Volcano plot showing differential ChrAcc sites between CreER⁺WT and CreER⁺dKO emergent ETPs sorted on day 7 of culture, by the criteria of ≥2 fold changes and FDR<0.05. B. Top motifs in differential ChrAcc sites in A, as determined with HOMER. C. Heatmaps showing connection of DEGs and differential ChrAcc sites between WT and dKO emergent ETPs. Heat map on the left displays grouping of Tcf1/Lef1-activated and -repressed genes based on their connection with ‘more closed’ and ‘more open’ ChrAcc sites, with gene numbers denoted in parentheses for each group and select genes marked in different colors, where color scale denotes z-normalized transcript levels with arbitrary units. In each group of genes, the associated differential ChrAcc sites were clustered according to the changes in ChrAcc (middle panel), with site numbers denoted in parentheses for each cluster, where color scale denotes normalized signal strength with arbitrary units. In the right panel, the presence of Tcf/Lef motif in the differential ChrAcc sites was detected with Jasper and marked with a horizontal red line, where the percentages denote the frequency of motif presence in each ChrAcc site cluster. D. ATAC-seq tracks at the indicated gene loci in CreER⁺WT and CreER⁺dKO emergent ETPs with both replicates shown. Cyan and orange open bars denote ‘more closed’ and ‘more open’ sites in dKO cells, respectively.