Mapping the developing human immune system across organs

Abstract

Single-cell genomics studies have decoded the immune cell composition of several human prenatal organs but were limited in describing the developing immune system as a distributed network across tissues. We profiled nine prenatal tissues combining single-cell RNA sequencing, antigen-receptor sequencing, and spatial transcriptomics to reconstruct the developing human immune system. This revealed the late acquisition of immune-effector functions by myeloid and lymphoid cell subsets and the maturation of monocytes and T cells before peripheral tissue seeding. Moreover, we uncovered system-wide blood and immune cell development beyond primary hematopoietic organs, characterized human prenatal B1 cells, and shed light on the origin of unconventional T cells. Our atlas provides both valuable data resources and biological insights that will facilitate cell engineering, regenerative medicine, and disease understanding.

Fig. 1. Cross-tissue cellular atlas of the developing human immune system.

(A) Overview of study design and analysis pipeline. We generated scRNA-seq and scVDJ-seq data from prenatal spleen, yolk sac and skin which were integrated via scVI with a collection of publicly available single-cell RNA-seq datasets. This cell atlas was used for (i) differential abundance analysis across gestation and organs with Milo, (ii) antigen receptor repertoire analysis with scirpy and dandelion, (iii) comparison with adult immune cells and in vitro differentiated cells with scArches and CellTypist, and (iv) spatial cell type deconvolution on Visium 10X data of hematopoietic and lymphoid organs using cell2location. (B) Summary of analyzed samples by gestational stage (x-axis) and organ (y-axis). Colors denote the types of molecular assays performed for each sample. The side bar indicates the total number of cells collected for each organ (after quality control). (C) Left: UMAP embedding of scRNA-seq profiles in prenatal tissues (908,178 cells) colored by broad cellular compartments. Right: bar plot of percentage of cells assigned to each broad compartment for each of the profiled organs. Raw cell proportions are adjusted to account for FACS-based CD45 enrichment. The category “Other” denotes clusters annotated as low-quality cells. (Eo/Baso/Mast: eosinophils/basophils/mast cells; ILC: innate lymphoid cells; NK: natural killer cells). (D) Representative colocalization patterns identified with non-negative matrix factorization of spatial cell type abundances estimated with cell2location. For each annotated microenvironment, we show (top) a dot plot of relative contribution of cell types to microenvironment (dot size) and (bottom) spatial locations of microenvironments on tissue slides, with the color representing the weighted contribution of each microenvironment to each spot. Each scale bar indicates a length of 1 mm.

Fig. 2. Myeloid variation across time and tissues.

(A) Beeswarm plot of log-fold change (x-axis) in cell abundance across gestational stages in Milo neighborhoods of myeloid cells. Results from five organs are shown. Neighborhoods overlapping the same cell population are grouped together (y-axis), and colored if displaying significant differential abundance (DA) (SpatialFDR 10%). The black dot denotes the median log-fold change. The top bar denotes the range of gestational stages of the organ samples analyzed. (B) Heat map of average expression across time of a selection of markers of stage-specific macrophage neighborhoods. Mean log-normalized expression for each gene is scaled (z-score). Gestational ages are grouped in 5 age bins. Age bins where less than 30 cells of a given subset were present are not shown. The top panel shows the fraction of all macrophages belonging to the specified macrophage population in each time point and each organ (color). (C) Close-up view of monocytes on Milo neighborhood embedding of myeloid cells (subset from fig. S16). Top: neighborhoods are colored by overlapping cell population. Bottom: neighborhoods displaying significant DA (SpatialFDR 10%) are colored by log-fold change in abundance between the specified organ and all other organs. (D) Mean expression of a selection of differentially expressed genes between CCR2hi monocytes from bone marrow and other organs. Log-normalized expression for each gene is scaled (z-score). We show genes upregulated in bone marrow associated with G2/M checkpoint and genes downregulated in bone marrow associated with TNF signaling (from MSigDB Hallmark 2020 gene sets). (E) Schematic of the proposed process of monocyte egression from the bone marrow mediated by CXCR4 and CCR2 expression: CXCR4^hi monocytes are retained in the bone marrow, until they switch to a proliferative state with increased expression of CCR2, mediating tissue egression. CCR2^hi monocytes seed peripheral tissues and then mature further to the periphery-specific IL1B expressing subtype.

Fig. 3. Lymphoid variation across time and tissues.

(A) Beeswarm plot of log-fold change (x-axis) in cell abundance across gestational stages in Milo neighborhoods of lymphoid cells (as in fig. 2A). (B) Heat map showing average expression across time of a selection of genes identified as markers of early-specific and late-specific NK neighborhoods (as in fig. 2B): NK cells identified in liver and skin before 12 pcw express TNF proinflammatory genes, whereas expression of immune-effector genes such as cytokines, chemokines and granzyme genes increases after 12 pcw. Age bins where less than 30 NK cells were present in a given organ are grayed out. (C) Close-up view of single-positive T cells on Milo neighborhood embedding of lymphoid cells. Each point represents a neighborhood, the layout of points is determined by the position of the neighborhood index cell in the UMAP in fig. S4I. Top: neighborhoods are colored by the cell population they overlap. Bottom: neighborhoods are colored by their log-fold change in abundance between the specified organ and all other organs. Only neighborhoods displaying significant differential abundance (SpatialFDR 10%) are colored. (D) Mean expression of a selection of differentially expressed genes between single-positive T cells from thymus (TH) and other organs. We show genes downregulated in the thymus associated with TNF signaling (using MSigDB Hallmark 2020 gene sets) and genes upregulated in the thymus associated with an IFN-α response. (E) Schematic of the proposed mechanism of thymocyte maturation and egression from thymus mediated by type I interferon signaling and NF-κB signaling.

Fig. 4. System-wide blood and immune cell development.

(A) Boxplots of the number of progenitor cells in all donors across organs. Each point represents a donor, color-coded by organ (YS: yolk sac; LI: liver; BM: bone marrow; TH: thymus; SP: spleen; MLN: mesenteric lymph node; SK: skin; GU: gut; KI: kidney). The red dash line marks the threshold of 10 cells for potential technical artefacts. B_prog: B cell lineage progenitors; M/E_prog: megakaryocyte/erythroid progenitors; MYE_prog: myeloid progenitors; T_prog: T cell progenitors. Detailed cell types included in each lineage are shown in table S5. Boxes capture the first-to-third quartile of the cell number and whisks span a further 1.5X interquartile range on each side of the box. (B) Multiplex smFISH staining with DAPI, CDH5 for endothelial cells, and VPREB1, DNTT, RAG1 for B cell progenitors in the human prenatal intestine at 15 pcw. Left panel shows a zoomed-out view with the area of interest boxed in white (scale bar: 500 μm). Right panel shows a detailed view of the area of interest (scale bar: 50 μm). Gray arrows point to B cell progenitors associated with blood vessels and orange arrows point to B cell progenitors away from blood vessels. (C) Scaled sum of abundances of B progenitor cell types estimated with cell2location, shown on representative slides for each organ, with the corresponding H&E staining. Each scale bar represents the length of 1 mm. (D) Cell type contributions to microenvironments containing B cell progenitors in different organs identified with non-negative matrix factorization of spatial cell type abundances estimated with cell2location. The color and the size of the dots represent the relative fraction of cells of a type assigned to the microenvironment.

Fig. 5. Identification of putative prenatal B1 cells.

(A) Left: Close-up view of non-progenitor B cell populations on UMAP embedding of all lymphoid cells (fig. S4I), with marker genes listed next to each cell type. Right: expression of B1 marker genes on UMAP. (B) Top: dot plot of IGHM and IGHD expressions in B1 and mature B cells, with color of dots representing the mean expression and size representing the fraction of cells expressing the gene. Bottom: cycling cell proportions within each B cell subtype colored by organs, with dot size representing log₁₀(cell count) and only dots with at least 10 cells shown. B1 cells had significantly higher cycling proportions than mature B cells in a logistic regression controlling for donors and organs. (C) Point plots of NP-addition length, CDR3 junction length, and mutation frequency in BCR heavy chains or light chains in B1 (N=2,357) and mature B (N=7,387) cells, with points representing the mean and lines representing 95% confidence intervals. Heavy chain VD and DJ junction NP-addition lengths are only calculated for cells with high-quality D gene mapping (B1: N=615, mature B: N=2,430). Difference in characteristics were tested with linear regressions controlling for donors and organs. (D) Volcano plot summarizing results of BCR heavy and light chain V, J gene segment usage comparison between B1 and mature B cells. The y-axis is the −log₁₀(Benjamini–Hochberg adjusted P-value) and the x-axis is log(odds ratio) computed using logistic regression controlling for donors and organs. (E) Normalized proportions of antibody-secreting cells in different sorted fractions of the ELISpot experiments (raw counts in table S6), colored by donor. Each point represents a reaction well. The proportions of antibody-secreting cells were normalized against the average proportion in CCR10^hi wells for each donor to remove donor-specific effects. A representative well image for each sorted fraction is shown on the bottom. (F) Schematic illustration summarizing the features of all human prenatal B1 cells, and additional features specific to CCR10^hi prenatal B1 cells.

Fig. 6. Deep characterization of human unconventional T cells.

(A) Proportions of cells expressing paired γδTCR, paired αβTCR, or both or neither. The proportions were calculated over cells that have had both single-cell αβTCR and γδTCR sequencing. The cells that expressed neither paired αβTCR nor paired γδTCR could be due to dropouts in single-cell TCR sequencing as over 50% of these contained orphan VDJ or VJ chains of αβTCR or γδTCR (fig. S28A). (B) Heat map showing the percentage of each γδTCR gene segment present in different T cell subtypes. Differential usage between cell subtypes was computed using the chi-squared test and gene segments with Benjamini–Hochberg adjusted P-values < 0.05 are marked with *. (C) Heat map showing the proportion of each TCRα gene segment present in different T cell subtypes. The gene segment usage in unconventional T cells and conventional T cells was compared using logistic regression, controlling for donors and organs. Gene segments with Benjamini–Hochberg adjusted P-value<0.05 are marked with magenta * (for preferential usage in unconventional T cells) and green * (for preferential usage in conventional T cells). (D) PCA plot summarizing TRAV, TRAJ, TRBV, TRBJ gene segment usage proportion in different T cell subtypes. Each dot represents a sample of at least 20 cells, with its size representing the cell count. The centroid of each cell type is shown as a filled circle, and 80% confidence contours are shown around the centroids. Arrows illustrate the proposed developmental trajectories. (E) Schematic illustration showing the T–T training origin of unconventional T cells in contrast to the T–TEC training origin of conventional T cells. (F) Top: schematic showing the experimental set-up of T cell differentiation from iPSCs in ATOs. Bottom left: UMAP visualization of different cell types in the ATO. Bottom right: density plots of cells from each time point over UMAP embedding. (G) Left: predicted annotations from logistic regression overlaid on the same UMAP plot as in (F); right: ZBTB16 expression pattern overlaid onto the same UMAP plot.