Spatial mapping of thymic stromal microenvironments reveals unique features influencing T lymphoid differentiation

Abstract

Interaction of hematopoietic progenitors with the thymic microenvironment induces them to proliferate, adopt the T lineage fate, and asymmetrically diverge into multiple functional lineages. Progenitors at various developmental stages are stratified within the thymus, implying that the corresponding microenvironments provide distinct sets of signals to progenitors migrating between them. These differences remain largely undefined. Here we used physical and computational approaches to generate a comprehensive spatial map of stromal gene expression in the thymus. Although most stromal regions were characterized by a unique gene expression signature, the central cortex lacked distinctive features. Instead, a key function of this region appears to be the sequestration of unique microenvironments found at the cortical extremities, thus modulating the relative proximity of progenitors moving between them. Our findings compel reexamination of how cell migration, lineage specification, and proliferation are controlled by thymic architecture and provide an in-depth resource for global characterization of this control.

Figure 1

Workflow for identification of stromal genes and generation of a spatial stromal map of the thymus.

Figure 2. Spatial relationships of subregions used in this study

Typical medullary and cortical subregions selected for collection are indicated. Panel A shows tissue pre-dissection, including baseline parameters for cortical and medullary tissue depth. Panel B shows the same tissue after defining regions for microdissection; each subregion in the cortex represents approximately 15% of the total cortical depth. Panel C shows the tissue after microdissection; in panel D, the lines defining the areas of interest have been removed to show the final result of the microdissection process. Scale bars on lower right of each panel are 200 μm.

Figure 3. Unsupervised clustering reveals stromal-dominant differences among different subregions of thymic tissue

RNA from microdissected tissues or purified lymphoid cells (as indicated) was used to probe microarrays, and 25831 probesets present in at least one microdissected tissue region were subjected to two dimensional hierarchical clustering (vertical hierarchy not shown). Black = highest normalized expression. Samples clustered in predictable fashion; notably, cortical subregions were distinct from each other, and were more similar to whole cortex than to other tissue or lymphoid types. Note that, in general, genes with the largest dynamic range are characteristic of medullary stromal cells (i.e., found in whole medulla but not medullary lymphocytes).

Figure 4. Representation of predicted functional classifications in stromal genes identified by differential mapping

Genes identified as stromal-specific (Fig. 1) in the cortex or the medulla were mapped onto pathways or functions as described in the text. Mapping onto classifications that would be expected to characterize stromal or lymphoid cells in the thymus are shown. Enrichment or depletion are measured against the number of genes from the whole chip that would be predicted to fall in each category based random chance. Stromal gene lists were enriched for categories that characteristic of stromal cells, and were depleted in categories that characterize lymphoid cells.

Figure 5. Spatial mapping of stromal gene regions in the thymus

Stromal genes that changed significantly between regions (p value <0.05 AND | fold change ≥ 2) were organized into clusters using the k-means algorithm, and 8 clusters representing the dominant patterns are shown. Grey bars indicate the interquartile range (25^th–75^th percentile) and solid vertical lines indicate the interdecile range (10^th–90^th percentile) for all genes in the cluster, while outliers are indicated by dots. For the sake of visual clarity, line graphs for individual genes (probesets) are not shown, but assignment to a given cluster mandates that each gene generally follows the trend indicated by the median value (bold magenta line). The per gene normalized value across all experiments (i.e., unity, from which relative change is measured) is indicated by a dashed red line. A panel of genes known to characterize stromal cells in specific compartments are indicated on the corresponding clusters.

Figure 6. Calculation of global stromal gene signatures by proportional subtraction of lymphoid signals from tissues

Panel A shows hypothetical behavior for noiseless data. Stromal-specific genes are infinitely higher in tissue than in lymphoid cells, and thus form a portion of the curve with infinite slope. Lymphoid-specific genes (no additional signal provided by stromal cells) likewise form the section of the curve with a slope of zero. The point at which the lymphoid-specific line intersects the y-axis defines the proportion of lymphoid signal in total tissue (per chip normalized data). In real (noisy) data, the hypothetical curve is deformed into a sigmoidal function in which the y-intercept does not exist (panel B). Note that only the region of the curve with tissue:lymphoid values <1 is shown, since this is the theoretical maximum contribution of lymphoid signal to tissue). Three independent methods were used to predict the y-intercept; the tangent to the inflection point (diagonal red line), a line perpendicular to the inflection point (horizontal red line), or the average ratio for a panel of known lymphoid-specific genes (dashed red line). Panel C shows enrichment in KEGG pathways for the most-changed genes in each compartment (medulla or cortex; yellow = highest statistical significance, black = threshold of significance, and blue = lowest significance). Note that enriched pathways in cortical stroma are strongly biased towards metabolic functions, while those in medulla are heavily biased towards signaling, particularly epithelial growth factor signaling pathways.