Lgr5+ telocytes are a signaling hub at the intestinal villus tip

The intestinal epithelium is a structured organ composed of crypts harboring Lgr5+ stem cells, and villi harboring differentiated cells. Spatial transcriptomics have demonstrated profound zonation of epithelial gene expression along the villus axis, but the mechanisms shaping this spatial variability are unknown. Here, we combined laser capture micro-dissection and single cell RNA sequencing to uncover spatially zonated populations of mesenchymal cells along the crypt-villus axis. These included villus tip telocytes (VTTs) that express Lgr5, a gene previously considered a specific crypt epithelial stem cell marker. VTTs are elongated cells that line the villus tip epithelium and signal through Bmp morphogens and the non-canonical Wnt5a ligand. Their ablation strongly perturbs the zonation of enterocyte genes induced at the villus tip. Our study provides a spatially-resolved cell atlas of the small intestinal stroma and exposes Lgr5+ villus tip telocytes as regulators of the epithelial spatial expression programs along the villus axis.


Introduction
The small intestine is a structured organ composed of repeating crypt-villus units. Lgr5+ epithelial stem cells at the base of the crypts continuously divide to give rise to proliferative progenitors that migrate towards the villi 1,2 . As they approach the crypt exits, these progenitors become post-mitotic and differentiate into distinct lineages -absorptive enterocytes, mucus-secreting goblet cells, enteroendocrine cells and tuft cells. The differentiated cells operate for about three days as they continue to migrate along the villus walls until they are shed from the villi tips 3 .
Villi cells have traditionally been considered 'terminally-differentiated' in the sense that they have irreversibly committed to a functional cell state, which they carry from birth to death. Recent studies challenged this view and demonstrated that villi epithelial cells constantly change their functional state, such that around 85% of enterocyte genes are expressed in a non-uniform manner along the villus axis 4 . Enterocytes first implement antimicrobial programs at the base of the villi, then shift to sequential absorption of amino acids, carbohydrates and lipids in distinct villi zones, finally up-regulating genes associated with cell adhesion and immune modulation at the villi tips. Similarly, enteroendocrine cells change the types of hormones they produce as a function of their position along the villus axis 5,6 . The epithelial lineages along the intestinal villus therefore exhibit profound spatial heterogeneity.
What are the mechanisms that facilitate the spatial heterogeneity of the villus epithelium?
One mechanism could be an internal 'clock', whereby epithelial cells are pre-programmed to turn distinct functions on and off at different times. A second mechanism could be a zonated response to spatial gradients of nutrients and bacteria at the luminal sides of the tissue. Yet a third mechanism could entail zonated molecular cues that originate in the lamina propria, the underlying stroma at the basal sides of the epithelial layer [7][8][9][10] . Here, we use spatially-resolved transcriptomics to characterize the zonated stromal gene expression signatures along the crypt-villus axis. We identify four mesenchymal cell populations residing at distinct crypt-villus zones. These include a sub-population of telocytes localized at the villus tip that is marked by Lgr5, a gene previously considere a specific marker of epithelial crypt stem cells. We demonstrate that Lgr5+ villus tip telocytes regulate the epithelial gene expression programs at the villus tip.

A spatial expression atlas of the small intestinal stroma
To obtain a global view of spatial heterogeneity of stromal gene expression we used Laser Capture Microdissection (LCM) to isolate four stromal zones along the jejunum crypt-villus axis from five mice (Fig. 1a,b). We performed RNA sequencing of these segments, yielding a coarse zonation map (Table S1). We focused on the zonation patterns of ligands and receptors 11,12 , since these are likely to implement zonated cross-talk with the epithelium (Fig. 1c, Fig. S1a,b, Table S2). Examples of ligand/receptor genes that exhibited higher expression in the crypt stroma included Grem1, encoding the BMP pathway inhibitor Gremlin1 13 (Fig. 1c) and Il18r1, the receptor for Il18. The cytokine Il18 has been shown to be expressed in enterocytes at the lower villi zones, as part of an anti-microbial gene module 4 .

Lgr5 is abundantly expressed in stromal cells at the villus tip
We next focused on stromal ligands and receptors that were zonated towards the villi tips in our stromal LCM-RNAseq data. These included Bmp4, a morphogen that was shown to inhibit epithelial proliferation 15 and to control enteroendocrine cell differentiation pathways at the tip 5 (Fig. 1c). Other villus tip molecules included the non-canonical Wnt ligand Wnt4 and the epidermal growth factor Egf. Egfr, the receptor for Egf, has elevated expression levels in villus tip enterocytes 4 , constituting another example of correlated spatial expression of pairs of epithelial and stromal ligands and their matching receptors (Table S2).
We further analyzed zonated proteins constituting the extracellular matrix (ECM, also termed the 'matrisome' 16 ), identifying distinct collagens and matrix metalloproteinases that are differentially zonated between the crypt and villus stroma (Fig. S1h). These included the crypt enriched Lama2, encoding the laminin subunit alpha-2 and the villus enriched Lama5, encoding the laminin subunit alpha-5, previously shown to be differentially zonated at the protein level 17,18 .
Our zonation analysis revealed Lgr5 to be one of the most highly expressed receptors in the villus tip stroma (Fig. 1c, Table S1). Lgr5 has been shown to be a specific marker of epithelial stem cells at the crypt base 2,14 . It was therefore unexpected to observe elevated expression levels of this gene at the lamina propria of the villus tip. To validate this finding, we performed smFISH and detected abundant localized expression of Lgr5 transcripts in PDGFRa+ telocytes that co-expressed Bmp4, a classic villus tip ligand 19 (Fig. 2). Telocytes are large mesenchymal cells with elongated extensions that stain positively for the PDGFRa surface marker. They form intricate contacts with both epithelial cells and other stromal cell types 9,20 . Telocytes that surround the crypts constitute important niche cells that secrete canonical Wnt morphogens and Rspo3 to maintain stemness of crypt Lgr5+ epithelial stem cells 9,21,22 . The molecular identities of telocytes residing in the villi stroma have thus far not been characterized.
The expression levels of Lgr5 in the villus tip telocytes (VTTs) were comparable to the expression levels in the epithelial crypt base columnar stem cells (mean expression of 0.109±0.012 mRNAs/µm -3 in tip telocytes vs. 0.104±0.008 mRNAs/µm -3 in crypt stem cells, (Fig. 2f). VTTs expressed both Lgr5 and its ligand Rspo3 23 (Fig. 2g). Notably, mRNAs of Rspo3 were localized along telopodes -thin PDGFRa+ extensions of the telocytes that extended towards the lamina propria and away from the cell body, where Lgr5 mRNAs were localized ( Fig. 2g).   percentiles of the smFISH expression, horizontal red lines are medians. g) Rspo3 mRNAs are localized on telopodes that extend away from the cell bodies of the VTTs. VTTs are marked by Lgr5 mRNA (red dots), Rspo3 mRNA (green dots) is localized away from the cell body, PDGFRa antibody mark VTTs cell bodies and telopodes. Scale bar -10 µm, in inset, green arrows point to Rspo3 mRNAs (green dots) localized on PDGFRa telopodes (blue). Telocyte cell body is marked by white dashed line. inset Scale bar -5 µm.
Several mouse models have been developed to investigate intestinal Lgr5+ crypt epithelial stem cells. These include models for ablating Lgr5 cells 24 and for tracing their progenies 2 .
We asked whether the knock-in constructs in these mice were also expressed in Lgr5+ VTTs.

Single cell transcriptomics reveals the molecular signature of VTTs
We next sought to obtain the global gene expression signatures of Lgr5+ VTTs and of other intestinal mesenchymal cell populations. LCMseq provides a spatial map of gene expression along the crypt villus axis, however each zone represents a mixture of diverse cell types, including mesenchymal cells, endothelial cells, immune cells 25 and enteric neurons 26 .
Unraveling the cellular origins of the zonated transcripts requires single cell transcriptomics.
Extracting telocytes is challenging due to their long thin extensions and their entrenchment within the ECM. We therefore used fluorescence activated cell sorting to enrich for PDGFRa+ cells that included the intestinal telocytes, as well as other mesenchymal cell types   Table S4 for the complete list of markers). Lgr5 is highly expressed in the Villus Tip Telocytes (7 out of 10 Lgr5+ cells belong to the VTT cluster, p=0.0121). f) Summed expression of the top markers for each mesenchymal cell cluster in the LCMseq data (Methods).
We focused on the mesenchymal cell clusters that were positive for PDGFRa. Reclustering these 329 cells yielded four distinct mesenchymal cell populations ( Fig. 4d-e, Table S3). We used Seurat 28 to identify distinct gene expression markers for each of the four cell populations (Fig. 4e, Fig. S3, Table S4). To identify potential zonation of these mesenchymal cell types we used our LCMseq map (Table S1) to examine the expression of the marker genes from each cluster (Fig. 4f, Methods). Using the marker gene identities and their spatial patterns, we annotated the clusters as crypt telocytes, enteric mesothelial fibroblasts, myofibroblasts and villus tip telocytes (VTTs). 10 cells out of the 329 mesenchymal cells were Lgr5+, 7 of which were in the VTT cell population (Fig. 4e, hypergeometric p=0.0121). The remaining 3 Lgr5+ cells in our scRNAseq data were interspersed with 1 cell in each of the remaining three cell populations.
We used smFISH to validate the spatial patterns of expression of the cell type markers revealed by the scRNAseq measurements. Grem1 and Sfrp1 expression was elevated in telocytes at the bottom of the crypts (Fig. S4a,c). Ly6a and Angptl4 expression was elevated in VTTs (Fig. S4b,c). Notably, although the scRNAseq data suggested higher levels of Bmp2 and Bmp4 ligands in crypt telocytes, our smFISH measurements showed significantly higher levels of Bmp2 and Bmp4 in VTTs, in line with the LCM measurements (Fig. 1c, Fig. S4c).
Enteric mesothelial fibroblasts 26 were marked by Thbs2, Fibin and Rgs5, as well as Angptl2, a gene previously shown to inhibit Bmp 29 . These were interspersed throughout the cryptvillus axis and were localized at the core of the villus, away from the epithelial layer ( Fig.   S4d), as were myofibroblasts, marked by Acta2, Tagln and Aldh1a1 (Fig. 4e).
Differential gene expression between crypt telocytes and VTTs further revealed elevated levels of the non-canonical Wnt ligand Wnt5a in VTTs (Fig. 5a), suggesting that VTTs implement a switch from canonical to non-canonical Wnt signaling 30,31 . Indeed, we observed broad expression of Axin2, a transcriptional readout of canonical Wnt signaling, along the villi epithelial cells, with a decrease at the villus tip (Fig. 5c), the zone of stromal Wnt5a expression.
Fig. 5 -VTTs implement a spatial switch from canonical to non-canonical Wnt signaling. a) Volcano plot demonstrating differential gene expression between VTTs and crypt telocytes. Black dots have qval<0.2 and expression fold change larger than 2 or smaller than ½. Wnt5a and Lgr5 are marked in red. b) Spatial shift in the stroma from the expression of canonical Wnt2b (red dots, marked by red arrows) to non-canonical Wnt5a (green dots, marked by green arrows). Scale bar -20 µm. c) Axin2 (red dots) is expressed broadly along the villus axis and repressed at the villus tip. Fos (green dots) is highly expressed in villus tip enterocytes. scale bar -20 µm. d) E-Cadherin protein (gray), encoded by Cdh1 gene, is induced in villus tip enterocytes. Large blobs in B-D are autofluorescent signals originating in immune cells.

VTT ablation perturbs epithelial gene expression at the villus tip
Enterocytes undergo profound gene expression changes as they approach the villus tip, a few hours before they are shed off into the lumen 4   b) Both VTTs and CBCs are ablated 24hr after DT administration, as evident by the lack of GFP fluorescence. c) GFP fluorescence re-appears in the crypt (white arrowheads) but not at the villus tip stroma 48hr after DT administration, indicating stable loss of VTTs. Inset on right shows a blow up of villus tip with no GFP+ VTTs. Green blobs are autofluorescent elements, also marked by green arrows. Scale bar in A-C -100 µm, insets scale bar -10 µm. d) Volcano plot demonstrating the changes in enterocyte gene expression 48hr following VTT ablation. Enterocyte villus tip genes that are reduced include Ada, Nt5e and Slc28a2 (red), composing the purine metabolism immune-modullatory tip module 4 . Black genes have q-values lower than 0.2 and max expression higher than 5*10 -6 (Methods). e) Enterocyte genes normally induced at the villus tip are reduced in expression 48hr following VTT ablation. Correlation between change in expression and expression zone -R=-0.53, p<10 -60 . Enterocytes were classified into six villus zones as in Moor et. al. 4 . f-h) smFISH validations of enterocyte villus tip genes that are changed 48hr following VTT ablation -in each panel Ctrl-left, DT-right. f) Cdh1. g) Egfr. h) Ada. Scale bar in f-i -20 µm i) Quantification of smFISH experiments, demonstrating that key epithelial villus tip genes are repressed 48hr after VTT ablation (Egfr, Cdh1, Klf4, Fos , Ada, Nt5e, Slc28a2 **p<10 -3 *p<0.05), whereas others remain unchanged (Creb3l3, p=0.5). Analysis was performed on two mice. Boxes show 25-75 percentiles of the smFISH expression, horizonatal red lines are medians.
Ablation of Lgr5+ crypt base columnar cells has been shown to initiate a regenerative program in the remaining crypt cells 32 . Importantly, migration of epithelial cells from the crypt to the villus tip takes 3-5 days 1,3 . We therefore argued that at short times following ablation of Lgr5+ cells, reduction in the expression of genes normally confined to the villus tip epithelium would be a result of loss of signals from the Lgr5+ VTTs, rather than incoming enterocytes, the state of which was perturbed in the crypt. To obtain a global view of the impact of VTTs on villus tip enterocyte gene expression, we therefore performed RNA sequencing of the epithelial layer 48 hours after DT administration (Table S5, Fig. 6d,e, Methods). We identified a significant repression in enterocyte genes that were normally zonated towards the tip (Spearman correlation of R=-0.53, p<10 -60 48hr post ablation between expression change following ablation and expression zone, Fig. 6e, Table S6). Villus tip enterocyte genes that were strongly reduced in expression upon VTT ablation included the purine metabolism module genes Slc28a2, Nt5e and Ada (Fig. 6d-i).

Additional prominent villus tip enterocyte genes include the transcription factors Klf4 and
Fos, Cdh1, encoding E-cadherin and Egfr 4 . These genes have a relatively high basal expression level either throughout the villi axis (Cdh1, Klf4 and Fos) or at the crypt (Egfr).
Expression changes in enterocyte at the villus tip could thus be masked for such genes in bulk RNA sequencing (Fig. 6d). To assess their expression changes, we therefore performed smFISH and measured expression specifically at the tip enterocytes, identifying a significant reduction of Egfr, Cdh1, Klf4 and Fos, in addition to Ada, Nt5e and Slc28a2 (Fig. 6g-i).
Creb3l3, an enterocyte gene that is elevated at the villi tip, did not exhibit changes in expression upon VTT ablation. These results demonstrate that VTTs are important regulators of the spatial expression programs of enterocytes at the villus tip, instructing the epithelial expression of key genes such as Egfr, Cdh1, Klf4, Fos, Nt5e, Ada and Slc28a2.
To assess the long-term consequences of VTT ablation we examined small intestinal tissue three weeks after VTT ablation. At this time point, VTTs re-appeared in 65% of the villi tips

Discussion
Our study exposed the spatial diversity of mesenchymal cells along the crypt-villus axis. We identified a cell population of Lgr5+ VTTs that form a highly localized niche for the villus tip epithelium, facilitating the massive transcriptional changes that enterocytes undergo before they are shed from the tissue 4 . VTTs are a source of Bmp ligands (Fig. 2d,e, Fig. S4c, Table   S1) and may implement a switch from canonical to non-canonical Wnt signaling, through the expression of the non-canonical Wnt ligand Wnt5a (Fig. 5b). Similar spatially antagonistic expression of stromal Wnt5a expression and canonical Wnt activity has been recently demonstrated in the mouse prostate 30 . This switch may be important for the induction of Ecadherin, a major component of the adherens junction that is strongly induced at the villus tip enterocytes (Fig. 5d) 33 . Wnt5a is also an activator of the AP1 transcription factors 34  Unlike the short-lived epithelial cells and the often mobile immune cells, VTTs are static, quiescent and long-lived, and can therefore integrate long-term information related to diets and pathologies. The intestine exhibits a remarkable ability to adapt to such perturbations 35,36 . Future work will resolve the roles of VTTs in shaping villus epithelial programs in such perturbed states. The re-appearance of VTTs 3 weeks following their ablation (Fig. 7) supports their importance for preserving the normal epithelial function. It will be important to identify the cellular origins of these re-appearing VTTs. Our study exposed VTTs as major regulators of the expression of villus tip enterocyte genes (Fig. 6e), yet enterocyte gene expression changes throughout the villus axis 4 . It will be interesting to explore the additional effects of zonated luminal signals and potential enterocyte 'internal clocks' in shaping enterocyte zonation.
Lgr5+ fibroblasts have previously been identified in the lung 37 Figure S1 -Zonated expression of stromal ligands (a) and receptors (b). Il18r1 marked in red. c) tSNE plot of immune cells in the small intestine taken from Biton et. al. 25

. Il18r1 is expressed in a subset of a cluster annotated as CD4 (red cells). d) Differential gene expression between the Il18r1+ cells (red cells in c)) and Il18r1-cells of the same cluster (green cells in c)
). These genes are markers of innate lymphoid cells type 3. The y axis was truncated at qval=10 -20 , leaving out Il18r1 (having expression ratio of 4,240, qval=10 -79 ). Vertical lines denote log2 ratios of 1 and -1. Horizontal line denotes qvalue of 0.01. e) smFISH for Il18 (green) and Il18r1 (red) demonstrating spatial colocalization at the lower villus zones. Scale bar -50µm. f) blow-up of a region boxed in e) demonstrating spatial adjacency of stromal immune cells with Il18r1 transcripts (red dots, cells marked by red arrows and outlined in red) and Il18+ villus bottom epithelial cells. Scale bar -10 µm. g) Quantification of the fraction of Il18r1+ cells (left y-axis) along the villus axis showing colocalization to the region where Il18 is highly expressed (right y-axis). h) Zonated expression of proteins constituting the extracellular matrix (ECM, also termed the 'matrisome' 16 ). Red arrows highlight the crypt-zonated Lama2 and villus tip-zonated Lama5, previously described to be inversely zonated 17,18 . Each row in (a,b,h) was scaled by its max expression across the four zones.

Supplementary Tables
Supplementary Table 1 -Average gene expression in the laser capture microdissected stromal zones. For each sample, UMI counts were normalized to the sum for all nonenterocyte genes that individually take up less than 1% of the total sample counts (Methods). TPM values presents the complete non-filtered data. Supplementary

Mice and tissues
All animal studies were approved by the Institutional Animal Care and Use Committee of

Antibodies used in this study
The following antibodies were used for FACS cell isolation: CD31 (PE-Cy7 102418,

Laser capture microdissection (LCM)
Tissue blocks for laser-capture microdissection were freshly prepared at the morning of the collection day. The jejunum part was transected on top of Wattman paper soaked in PBS to support the opened structure in petri dish on ice. The tissue was washed with cold PBS and embedded in OCT on dry ice. LCM protocol was applied as previously described 4 (Fig. 1). An area of ~30,000 μm 2 per zone was collected for each mouse by pooling 9-11 distinct villi from a single tissue section. For mice 1-3 (Table S1) we additionally catapulted crypt stroma and extracted material from two sequential tissue sections for each zone, with the exception of the villus center zone of mouse 3, where one section failed to amplify (Table S1, TPM value tab).

SMART-Seq for bulk LCM samples
RNA libraries from the bulk tissues were prepared using SMART-Seq v4 Ultra Low Input RNA Kit (Clontech, 634888), with a 15 PCR cycles for the amplification step. Subsequent steps were applied as mentioned in the protocol. Nextera XT DNA Library Prep Kit (Illumina) was use to finalize the libraries. Library concentration and quality control were determined using NEBNext Library Quant Kit (New England Biolabs) and Agilent High Sensitivity D1000 ScreenTape System (Agilent, 5067-5584). Library final concentration of 2.4pM was loaded on NextSeq550 (Illumina) sequencing machine aiming for 20M reads per sample with the following cycle distribution: 38bp read1, 8bp index1, 8bp index2, 38bp read2.

Stroma LCM data processing
We removed 31 enterocyte genes, determined as genes the summed fraction of which makes up the top 50% of enterocyte expression. To this end we considered the maximal expression of all genes in enterocytes over the zones as measured in Moor et al. 4 . We normalized the TPM values upon removing these genes by the sum of all remaining TPM.
We next divided the gene expression values in each sample by the sum for all genes that individually make up less than 1% of the sample's summed expression values. For each of the four zones, we computed the means and standard errors of the means over the different samples. Ligand receptor pairs were extracted from Ramilowsky et al. 11 , Matrisome components were extracted from Naba et al. 16 . Only genes with maximal zonation value larger than 5*10 -6 were considered.

Analysis of Il18r1 stromal cells
We used the single cell RNAseq data of Biton et al. 25 (6558 cells annotated as WT control, Fig. S1c). Il18r1 positive cells belonged to a cluster annotated as 'CD4'. We performed differential gene expression using Wilcoxon rank-sum tests between the Il18r1+ and Il18rcells in that cluster (Fig. S1d). We used Immgen My Geneset tab (https://www.immgen.org) to annotate the 200 genes that had the highest fold-change between Il18r1+ cells and Il18r1-cells within the CD4-annotated cluster (only genes with expression above 5*10 -5 were used for this analysis). The gene set was enriched in innate lymphoid cells type 3.

Hybridization and imaging
Probe libraries were designed using the Stellaris FISH Probe Designer Software (Biosearch Technologies, Inc., Petaluma, CA). 7µm thick sections of fixed Jejunum were sectioned onto poly L-lysine coated coverslips and used for smFISH staining. The intestinal sections were hybridized with smFISH probe sets according to a previously published protocol 14 . SmFISH probe libraries (Table S7)

Single cell isolation
Due to the difficulties of isolating intestinal stromal cells we applied several dissociation protocols and surface marker staining to obtain a broad sampling of cells. For all mice, cells were isolated from the jejunum. The jejunum was extracted and rinsed in cold PBS. The tissue was opened longitudinally and sliced into small fragments roughly 2 cm long and incubated in 10mM EDTA-PBS on ice for 10 min. The tissue was then moved to warm 5ml 10mM EDTA-PBS containing Liberase TM (100 μg/mL, Sigma) and DNaseI (2 U/mL,Sigma) and incubated at 37°C for 20 min while shaken vigorously every few minutes. At the end of the incubation time, 5ml of cold PBS was added to the cell suspension. The supernatant was filtered through a 100μm filter and centrifuged at 300 g for 5 min. the pellet was resuspended in FACS buffer (2mM EDTA, 0.5% BSA in PBS) and stained with the required antibodies for flow cytometry sorting.
Isolation of telocytes was performed as previously decribed in 9

Single-cell sorting
Single cells were sorted with SORP-FACSAriaII machine (BD Biosciences) using a 100 μm nozzle. Dead cells were excluded on the basis of 500 ng/ml Dapi incorporation. Sorte cells were negative for CD45 and EpCAM and positive for PDGFRa. Cells were sorted into 384well cell capture plates containing 2 μl of lysis solution and barcoded poly(T) reversetranscription (RT) primers for single-cell RNA-seq 27 . Barcoded single cell capture plates were prepared with a Bravo automated liquid handling platform (Agilent) as described previously 27 . Four empty wells were kept in each 384-well plate as a no-cell control during data analysis. Immediately after sorting, each plate was spun down to ensure cell immersion into the lysis solution, snap frozen on dry ice and stored at -80°C until processed.

Massively Parallel Single Cell RNA-Seq (MARS-Seq) library preparation
Single cell libraries were prepared, as described in 27 . Briefly, mRNA from cells sorted into MARS-Seq capture plates were barcoded and converted into cDNA and pooled using an automated pipeline. The pooled sample was then linearly amplified by T7 in vitro transcription and the resulting RNA was fragmented and converted into sequencing ready library by tagging the samples with pool barcodes and Illumina i7 barcode sequences during ligation, reverse transcription and PCR. Each pool of cells was tested for library quality and concentration was assessed as described in 27 . Machine raw files were converted to fastq files using bcl2fastq package, to obtain the UMI counts, reads were aligned to the mouse reference genome (GRCm38.84) using zUMI packge 44

scRNAseq data processing
For each single cell and for each gene we first subtracted the estimated background expression. Background was calculated for each 384-well plate separately, as the mean gene expression in the four empty wells. After subtraction, negative values were set to zero.
Next, cells with total UMI counts lower than 400 or higher than 8,000 and total gene counts lower than 250 were removed. We used Seurat v2.3.4 package in R 28 v3.5.3 to visualize and cluster the single cell RNAseq data (Fig. 4, Fig. S3). Gene expression measurements (UMIs per gene) were normalized for each cell by the summed UMI, multiplied by a scale factor 10,000, and then log-transformed. To avoid undesired sources of variation in gene expression, we used Seurat to regress out cell-cell variation driven by total number of UMIs, and mitochondrial genes fraction. For detection of variable genes, we set a bottom cutoff of 0.25 and a top cutoff of 4 on the regressed log-transformed average gene expression, as well as a bottom cutoff of 0.5 on the dispersion. Cell clustering was based on PCA dimensionality reduction using the first 11 PCs, and a resolution value of 2.
We used cell type-specific markers to interpret the resulting 3 main clusters, Epcam was highly expressed in the epithelial cells clusters, Ptprc was highly expressed in the immune cells clusters and Pdgfra was highly expressed in the fibroblast cells clusters.
We next focused on the cell clusters that had an average expression of Pdgfra that was higher than 10 -4 per cell (fraction of UMI counts). These included 329 cells. For this new set of cells, we re-ran Seurat with the following parameters; For detection of variable genes, we set a bottom cutoff of 0.25 and a top cutoff of 4 on the regressed log-transformed average gene expression, as well as a bottom cutoff of 1 on the dispersion. Cell clustering was based on PCA dimensionality reduction using the first 5 PCs, and a resolution value of 0.7. Marker genes were detected with Seurat FindAllMarkers function with parameters min.pct=0.2, logfc.threshold=0.25 (Table S4).
To identify zonation of the four mesenchymal cell types we examined the summed expression of concise sets of markers that were specific to each cell type. To select these markers, we included genes that had an average expression higher than 5*10 -6 and were expressed in at least 10 cells in the respective cell type. We sorted these genes by the fold change of their expression levels compared to the maximal average expression in the other three mesenchymal cell types. We included the top 20 genes and further excluded genes for which the fold change ratio was lower than 5-fold. Differential gene expression between crypt telocytes and VTTs (Fig. 5a) considered all genes with average expression larger than 10 -5 of cellular UMIs and more than 5 cells positive for the gene in at least one of the two populations. Storey's method was used to compute q-values.

Lgr5+ cell ablation
For Lgr5+ cell ablation, male Lgr5 DTRGFP mice aged 4-5 months were administered 50 g/kg diphtheria toxin (322326 Sigma) or saline vehicle intraperitoneally. Sample were collected 24 or 48 hr post injection. For the 24 hr time points, four DT-injected and two mock-injected mice were used for histological examination (Fig. 6). For the 48 hr time point, four DTinjected and four mock-injected mice were used for histological examination, smFISH and epithelial cell sequencing. Additionally, two DT-injected and two mock-injected mice were sacrificed 20 days after administration and used for histological examination and smFISH.
For epithelial cell isolation, 3 cm from the proximal jejunum were extracted and rinsed in cold PBS. The tissue was opened longitudinally and sliced into small fragments roughly 2 cm long and incubated in 10mM EDTA-PBS on ice for 20 min. The tissue was then moved to warm 5ml 10mM EDTA-PBS and incubated at 37°C for 5 min while shaken vigorously every few minutes. At the end of the incubation time, 5ml of cold PBS was added to the cell suspension. The supernatant was filtered through a 100μm filter and centrifuged at 300 g for 5 min. the pellet was snap frozen in liquid nitrogen and was processed for bulk RNA sequencing. For smFISH blocks, 5 cm from the distal jejunum was fixed immediately in 4% FA and processed as mentioned before.

Bulk RNA sequencing of ablation experiment samples
Snap frozen cells were thawed into TRI reagent (sigma), RNA was isolated by Direct-zol RNA MiniPrep kit (Zymo research) according to the manufacturer instructions. RNA was processed by the mcSCRBseq protocol 45 with minor modifications. RT reaction was applied on 10 ng of total RNA with a final volume of 10µl (1x Maxima H Buffer, 1mM dNTPs, 2µM TSO* E5V6NEXT, 7.5% PEG8000, 20U Maxima H enzyme, 1µl barcoded RT primer).
Subsequent steps were applied as mentioned in the protocol. Library preparation was done using Nextera XT kit (Illumina) on 0.6ng amplified cDNA. Library final concentration of 2.2pM was loaded on NextSeq 500 (Illumina) sequencing machine aiming for 20M reads per sample. Raw files were converted to FASTQ files using bcl2fastq package, to obtain the UMI counts, fastq reads were aligned to the mouse reference genome (GRCm38.84) using zUMI packge 46 with the following parameters RD1 16bp, RD2 66bp with a barcode (i7) length of 8bp. UMI counts were processed with the edgeR package 47 . Our analysis included four intestinal samples of DT-injected mice and four from mock-injected, 48 hrs after injection (Table S5). We first filtered the genes by expression to maintain only genes with expression greater than 10 -5 of the summed UMI counts of the sample in at least 2 samples and next ran the calcNormFactors function with default settings. We used the "robust" option in the glmQLFit function to robustly estimate the QL dispersion and glmQLFTest with default parameters to compute differential gene expression with Benjamini-Hochberg correction for multiple hypotheses (Table S6). Fig. 6d shows the volcano plots of all genes with with maximal expression higher than 5*10 -6 in at least one of the samples, as well as maximal zonation higher than 5*10 -6 in Moor et al. 4 .

Data availability
All data has been deposited in GEO and will be publically available upon publication.

Code availability
All codes used in this study will be available upon request.