1. Developmental Biology
Download icon

Bi-fated tendon-to-bone attachment cells are regulated by shared enhancers and KLF transcription factors

  1. Shiri Kult
  2. Tsviya Olender
  3. Marco Osterwalder
  4. Svetalana Markman
  5. Dena Leshkowitz
  6. Sharon Krief
  7. Ronnie Blecher-Gonen
  8. Shani Ben-Moshe
  9. Lydia Farack
  10. Hadas Keren-Shaul
  11. Tomer-Meir Salame
  12. Terence D Capellini
  13. Shalev Itzkovitz
  14. Ido Amit
  15. Axel Visel
  16. Elazar Zelzer  Is a corresponding author
  1. Department of Molecular Genetics, Weizmann Institute of Science, Israel
  2. Environmental Genomics and Systems Biology Division, Lawrence Berkeley National, United States
  3. Department for BioMedical Research (DBMR), University of Bern, Switzerland
  4. Life Sciences Core Facilities, Weizmann Institute of Science, Israel
  5. Department of Immunology, Weizmann Institute of Science, Israel
  6. Department of Molecular Cell Biology, Weizmann Institute of Science, Israel
  7. Department of Human Evolutionary Biology, Harvard University, Department of Human Evolutionary Biology, United States; Broad Institute of Harvard and MIT, United States
  8. U.S. Department of Energy Joint Genome Institute, Lawrence Berkeley National Laboratory, United States
  9. School of Natural Sciences, University of California, Merced, United States
Research Article
  • Cited 4
  • Views 1,388
  • Annotations
Cite this article as: eLife 2021;10:e55361 doi: 10.7554/eLife.55361

Abstract

The mechanical challenge of attaching elastic tendons to stiff bones is solved by the formation of a unique transitional tissue. Here, we show that murine tendon-to-bone attachment cells are bi-fated, activating a mixture of chondrocyte and tenocyte transcriptomes, under regulation of shared regulatory elements and Krüppel-like factors (KLFs) transcription factors. High-throughput bulk and single-cell RNA sequencing of humeral attachment cells revealed expression of hundreds of chondrogenic and tenogenic genes, which was validated by in situ hybridization and single-molecule ISH. ATAC sequencing showed that attachment cells share accessible intergenic chromatin areas with either tenocytes or chondrocytes. Epigenomic analysis revealed enhancer signatures for most of these regions. Transgenic mouse enhancer reporter assays verified the shared activity of some of these enhancers. Finally, integrative chromatin and motif analyses and transcriptomic data implicated KLFs as regulators of attachment cells. Indeed, blocking expression of both Klf2 and Klf4 in developing limb mesenchyme impaired their differentiation.

Introduction

The function of the musculoskeletal system relies on the proper assemblage of its components, namely skeletal tissues (bone, cartilage, and joints), muscles, and tendons. However, the attachment of tissues composed of materials with large differences in their mechanical properties is highly challenging. In the musculoskeleton, elastic tendons, which have a Young’s modulus (a measure of stiffness) in the order of 200 megapascal, are attached to the much harder bone, with a modulus in the order of 20 gigapascal. This disparity makes the connection between these two tissues a mechanical weak point, which is subject to higher incidence of tearing by both external and internal forces acting on the musculoskeleton during movement. The evolutionary solution to this problem is the enthesis, a transitional tissue that displays a gradual shift in cellular and extracellular properties from the tendon side through to the bone side (Genin et al., 2009; Liu et al., 2011; Liu et al., 2012a; Lu and Thomopoulos, 2013; Thomopoulos et al., 2003). Yet, despite its importance, the formation of this cellular gradient as well as the underlying molecular mechanism remain largely unknown.

In recent years, the initial events that lead to the formation of the embryonic attachment unit (AU), which serves as the primordium of the enthesis, have started to be investigated. These studies identified the progenitors of the AU and showed that they express both the chondrogenic and tenogenic transcription factors Sox9 and scleraxis (Scx), respectively (Blitz et al., 2013; Sugimoto et al., 2013). The patterning of the Sox9+/Scx+ progenitors along the skeleton is regulated by a genetic program that includes several transcription factors (Eyal et al., 2019). Next, Sox9+/Scx+ progenitors differentiate to chondrocytes, which form a bone eminence on one side, or to tenocytes, which form the tendon on the other side, whereas the cells in between differentiate into Gli1+ cells that eventually will form the enthesis (Felsenthal et al., 2018; Schwartz et al., 2017; Schwartz et al., 2015).

Both molecular and mechanical signals regulate the AU. TGFβ signaling regulates the specification of AU progenitors, whereas BMP and FGF signaling as well as mechanical signals determine their fate and differentiation (Blitz et al., 2013; Blitz et al., 2009; Roberts et al., 2019). Postnatal enthesis cells have been termed fibrocartilage cells based on their histological appearance, because they display morphological features that are shared with tenocytes and chondrocytes (Thomopoulos et al., 2010). In recent years, several studies have identified some of the genes that these cells express, including collagens type I, II, and X; Indian hedgehog (Ihh); parathyroid hormone-related peptide (Pthlh); patched 1 (Ptc1); runt-related transcription factor 2 (Runx2); tenascin C (Tnc); and biglycan (Bgn) (Liu et al., 2012b; Liu et al., 2013; Liu et al., 2018; Thomopoulos et al., 2010). Interestingly, these genes are also expressed by cells in the neighboring tissues, namely by chondrocytes or tenocytes. However, despite these advances, a comprehensive molecular signature of this tissue and the mechanism that enables its formation are still missing.

In this work, we aimed to decipher the identity of the fibrocartilage cells that form the attachment tissue between tendon and bone. Bulk and single-cell transcriptomic analyses of the attachment cells, which were validated by in situ hybridization (ISH) and single-molecule fluorescent ISH (smFISH), showed that these cells express a mix of the transcriptomes of chondrocytes and tenocytes. Chromatin analysis further verified the transcriptomic results and provided a mechanistic explanation for the bi-fated behavior of attachment cells, which share enhancers with their neighboring tenocytes or chondrocytes. Finally, we identify the transcription factors KLF2 and KLF4 as regulators of attachment cell differentiation. Overall, we provide the transcriptional as well as the epigenetic mechanism that allows attachment cells to activate a combination of cartilage and tendon transcriptomes and, thereby, the formation of the unique transitional tissue.

Results

Attachment cell transcriptome is a mix of chondrocyte and tenocyte transcriptomes

To date, the transcriptome of attachment cells has not been characterized thoroughly. We therefore analyzed the transcriptome of embryonic day (E) 14.5 attachment cells from the prominent deltoid tuberosity and greater tuberosity of the humerus. With the goal to isolate these cells specifically, we generated a compound mouse by crossing three mouse lines, namely Col2a1-Cre, R26R-tdTomato, and Scx-GFP (see Materials and methods) (Blitz et al., 2013; Sugimoto et al., 2013). Thus, the fluorescent reporter tdTomato labeled Col2a1-expressing chondrocytes, whereas GFP fluorescently labeled Scx-expressing tenocytes. Unexpectedly, the two reporters failed to label the attachment cells that were located in between these two populations. This failure might be due to a missing regulatory element in one of the constructs that were used to produce each transgenic reporter. Nevertheless, the borders between tendon and attachment cells and between cartilage and attachment cells were clearly demarcated (Figure 1A, Figure 1—figure supplement 1A,B). We therefore used laser capture microdissection (LCM) to subdivide the attachment site into three cellular compartments, namely attachment cells, adjacent tenocytes, and adjacent chondrocytes. As controls, samples were also taken from two more compartments, remote tenocytes and remote chondrocytes. Initial analysis of the different transcriptomes using principal components analysis (PCA) showed that the transcriptomes of tenocytes and chondrocytes were clearly separated, whereas attachment cells were located between the two cell types, recapitulating their anatomical positions. This suggests that the attachment cell transcriptome is shared with both chondrocytes and tenocytes (Figure 1B, PC1 31.66%).

Figure 1 with 4 supplements see all
Transcriptomic analysis of tendon-to-bone attachment site domains at E14.5.

(A) A scheme of tendon-to-bone attachment site of the murine deltoid tuberosity. (B) Principal component analysis (PCA) of bulk MARS-seq data from E14.5 attachment site samples. The x-axis (PC1) shows the highest variance among the samples. Interestingly, the samples are arranged according to their anatomical locations. ‘R’ (samples 1 and 5) stands for remote and ‘A’ (samples 2 and 4) is for adjacent. The y-axis (PC2) shows that tenocytes and chondrocytes are closer to one another, while attachment cells (black circle) were found to be remote from both of them, that is with higher variance, suggesting a unique gene expression profile. (C) Heatmap of gene expression profiles at E14.5 shows 374 selected genes that exhibited differential expression between tenocytes and chondrocytes and were also expressed by attachment cells. Color bar (−1.5-0-1.5) represents the log-normalized counts standardized per gene, as yellow is higher than the mean (0) and blue is lower than the mean. Attachment cells display a gradient of gene expression profiles, reflecting their function as a transitional tissue. The upper cluster contains genes highly expressed in tenocytes (e.g. Co1a1, Col1a2, Col5a1, Scx, Bgn), whereas the lower cluster contains genes highly expressed in chondrocytes (e.g. Sox9, Col2a1, Acan). Top list on the right contains genes found to be expressed in attachment cells and in tenocytes, whereas bottom list contains genes expressed in attachment cells and chondrocytes; genes in bold type are known tenocyte or chondrocyte markers. D (a-l). Double-fluorescent ISH for mRNA of tendon (Igfbp5, biglycan, Col5a1, Col1a1) and cartilage (Col11a1 or Wwp2) genes shows that attachment cells (in yellow, shown by arrows) exhibit an in vivo gene expression profile that combines tendon and cartilage genetic programs. a-l: X20 magnification, scale bar: 50 µm; a’-l’: magnification of upper panels.

To further support our initial observation that the transcriptome of the attachment cells is a mixture of chondrocyte and tenocyte transcriptomes, we clustered the statistically significant differentially expressed genes between all samples into five clusters, using CLICK (Figure 1—figure supplement 2 and see Materials and methods). Out of 865 identified genes, 735 genes were found in two clusters. The first cluster contained mainly known tenogenic genes and the second contained chondrogenic genes. From these two clusters, 374 genes, 320 of them tenogenic and 54 chondrogenic, were also found to be expressed by attachment cells. They included major regulators and marker genes of the two tissues, such as Sox9, Col2a1, and Acan for chondrocytes and Col1a1, Col1a2, Scx, and Col5a1 for tenocytes (Figure 1C).

The third cluster comprised 54 genes that were found to be upregulated in cartilage adjacent to attachment cells alone (Figure 1—figure supplement 2). Interestingly, our analysis identified 24 and 23 genes that were found to be down- or upregulated in attachment cells, shown by the fourth and fifth clusters, respectively (Figure 1—figure supplement 2). The genes that were found to be uniquely upregulated in attachment cells included transcription factors, such as the Krüppel-like factors (KLFs), Lmo1, and Gli1, which could act as regulators of the genetic program of attachment cells. In addition, this set included differentiation markers such as Thy1, regulators of bone for example Acp5 and Alpl, protein kinases such as Mapk12 and Mast2, and signaling molecules such as Nod, Traip, Aplnr and others (Figure 1—figure supplement 3A,B). GO analysis of these genes yielded terms relating to regulation of cytokine and IL-12 production, as well as response to laminar fluid shear stress (Figure 1—figure supplement 3C). These results clearly show that the transcriptome of the attachment cells includes a mixture of tenocyte and chondrocyte genes, many of which are involved in ECM organization and developmental processes, in addition to a unique subgroup of genes that are upregulated in these cells. To determine the relative proximity of attachment cells to the other two cell types, we performed hierarchical clustering on the 865 differentially expressed genes. Results showed that these cells are closer to tenocytes (Supplementary file 1).

To validate our transcriptome analysis, we performed single- and double-fluorescent in situ hybridization (FISH) using marker genes for tenocytes and chondrocytes that were selected from the transcriptomic results, namely Igfbp5, biglycan (Bgn), Col5a1, and Col1a1 for tenocytes and Col11a1 and Wwp2 for chondrocytes. As seen in Figure 1D and Figure 1—figure supplement 4, in agreement with the transcriptome analysis, the selected markers were co-expressed by the attachment cells.

Single-cell RNA-seq reveals the bi-fated nature of attachment cells

Finding that the bulk transcriptome of attachment cells represents a mixture of tenocyte and chondrocyte genes led us to examine if this phenotype also exists at a single-cell resolution by performing scRNA-seq. To isolate E13.5 Sox9+/Scx+ attachment progenitors, we generated a compound mouse line harboring Sox9-CreER, tdTomato and Scx-GFP transgenes (Figure 1—figure supplement 1C–E). Then, FACS-sorted attachment progenitors from the proximal side of the forelimb underwent single-cell 10x Genomics RNA sequencing (Figure 2A,B and see Materials and methods). Sequenced progenitors were then filtered, normalized, and clustered (Seurat 3.1.5).

Figure 2 with 1 supplement see all
Attachment cells co-express tendon and cartilage genes at the single-cell level.

(A) scRNA-seq analysis of E13.5 Sox9+/Scx+ attachment progenitors is shown as uniform manifold approximation and projection (UMAP) embedding of jointly analyzed single-cell transcriptomes (10x Chromium platform). (B) Violin plots of distinct genes associated with tendon or cartilage. (C) Single-molecule fluorescent ISH (smFISH) of mRNA of tendon biglycan (Bgn, red) and cartilage Wwp2 (green) genes on the background of DAPI staining (blue) further validates the dFISH results. X100 magnification, scale bar: 10 µm. (D) Quantification of Bgn and Wwp2 smFISH results in cartilage, attachment, and tendon cells.

Results showed five transcriptionally distinct subpopulations (Figure 2A). Examination of the top-10 differentially expressed genes across clusters revealed that they were subdivided to three main categories. Clusters 1 and 3 were highly rich with chondrogenic genes (Sox9, Sox5, Col2a1, Wwp2), whereas clusters 2 and 4 were highly rich with tenogenic genes (Scx, Hic1, Ptn, Igfbp5). By contrast, the identity of cluster 0 was less specified, as it was highly rich with genes such as Col1a1, Alpl, Spp1, and Runx2 in addition to Bicc1, Maf, and Mef2c (Figure 2—figure supplement 1).

To study the bi-fated nature of attachment cells, we analyzed the relative expression of tenogenic (Scx, Col1a2, biglycan) and chondrogenic (Sox9, Col2a1, Wwp2) marker genes in each subpopulation. As seen in Figure 2B, all five clusters contained both chondrogenic and tenogenic marker genes. However, while tenogenic genes were relatively high in all five clusters, the expression levels of chondrogenic markers were more variable, being the lowest in cluster 2. Overall, the scRNA-seq results suggest that the Sox9+/Scx+ attachment progenitors are bi-fated. The analysis also revealed a molecular sub-classification of this progenitor population.

Finally, to confirm the co-expression of tenogenic and chondrogenic genes at the attachment site at a single-cell resolution, we performed single-molecule FISH for Wwp2 (chondrogenic marker) and Bgn (tenogenic marker). As seen and quantified in Figure 2C,D, Wwp2 and Bgn were indeed co-expressed in the attachment cells.

Overall, these results indicate that attachment cells are bi-fated, expressing in parallel both chondrogenic and tenogenic genes (referred to in the following as mixed transcriptome or mixed expression profile).

Genome-wide profiling of attachment cell-specific regulatory regions

To gain a mechanistic understanding of how attachment cells activate a combination of two transcriptomes, we compared chromatin accessibility in these cells with open chromatin signatures defining chondrocytes and tenocytes by conducting an assay for transposase-accessible chromatin with high-throughput sequencing (ATAC-seq) (Buenrostro et al., 2015). This method allows to profile open chromatin regions, some of which may act as enhancers. To isolate by FACS E13.5 humeral chondrocytes, tenocytes and attachment cells, we generated a compound mouse line harboring Sox9-CreER, tdTomato, and Scx-GFP transgenes (Figure 1—figure supplement 1C–E). However, because the number of isolated chondrocytes was insufficient, chondrocytes were FACS-sorted from E13.5 Col2a1-CreERT-tdTomato-Scx-GFP mouse. These three cell populations were then subjected to ATAC-seq (see Materials and methods).

Initial PCA analysis of accessible chromatin profiles for each FACS-sorted cell population once again revealed that tenocytes and chondrocytes were clearly separated, while attachment cells resided between these two cell types (Figure 3—figure supplement 1A, Figure 1B). Next, we compared global chromatin accessibility among the three cell types by calculating the level of overlap among the ATAC-seq peaks (Figure 3A). While the majority of the peaks were shared by all three cell types, attachment cells had a significantly lower number of unique peaks (p<1e-4 relative to both cell types, chi-square with Yates correction), and a significantly higher overlap with the other two cell types (p<2.2e-1).

Figure 3 with 1 supplement see all
Accessible chromatin reveals an epigenetic mechanism shared by attachment cells and neighboring tenocytes or chondrocytes.

(A) Venn diagram showing cell-specific or overlapping peaks of ATAC-seq among tenocytes, chondrocytes and attachment cells. (B) Heatmap of ATAC-seq peaks associated with E14.5 differentially expressed genes. Left: TSS peaks, right: intergenic or intron peaks. The peaks are sorted according to their degree of accessibility across the three cell types. (C) Percentage of common peaks (shared by three cell types) vs. differential peaks (the chromatin is open only in one or two cell types) compared between TSS and intergenic areas (p=0, chi-square test). (D) IGV snapshots of the TSS region of Mgp, Sox9, Col11a1, and Col1a2 genes, as well as potential enhancers of these genes.

Analysis of the ATAC-seq signal revealed that 13,017 peaks were located near transcription start sites (TSSs), whereas 31,856 peaks were in intergenic or intron regions. Most of the peaks that were located near TSSs were accessible in all three cell types (87%), and only 13% were accessible in one or two cell types (Figure 3—figure supplement 1C,D).

Next, we studied the ATAC-seq signal of peaks associated with the genes that were differentially expressed at E14.5, using HOMER default parameters. We found that 819 peaks were located near transcription start sites (TSSs), whereas 2340 peaks were in intergenic or intron regions. Most of the peaks that were located near TSSs (708, 86%) were accessible in all three cell types, and only 111 (13%) were accessible in one or two cell types (Figure 3B,C). This low level of differential accessibility is inconsistent with the possibility that promoter accessibility is the main mechanism regulating the bi-fated attachment cells. Interestingly, a significantly higher fraction of intergenic peaks were specific to one or two cell types (1767, 74.7%, p=0, chi-square test, Figure 3B,C).

Out of these 1767 intergenic peaks, 920 peaks were accessible in attachment cell; 672 were shared with either tenocytes or chondrocytes; and 248 were accessible only in attachment cells (Supplementary file 3). These results suggest that the 920 intergenic elements that were accessible in attachment cells may act as enhancers that drive the transcriptome of these cells. Moreover, since most of these intergenic elements were shared between attachment cells and either chondrocytes or tenocytes, they may serve as part of the mechanism that drives the mixed transcriptome of attachment cells.

To identify such dual cell-type-specific enhancers likely regulating attachment cell differentiation, we next screened for shared enhancers of 15 bona fide markers of tenocytes or chondrocytes that were found to be expressed in E14.5 attachment cells (Figure 1C). To improve the prediction of these enhancers, we selected our ATAC-seq peaks based on their proximity to genes with verified expression in attachment cells and another cell type (Figure 1—figure supplement 4, Figure 1D) and computationally intersected them with ENCODE datasets of histone modification marks associated with enhancers and promoters in mouse limbs at E13.5 (H3K27Ac, H3K9ac, H3K4me3, H3K4me1 ChIP-Seq), and other datasets (Andrey et al., 2017; Guo et al., 2017; Figure 4 and Table 1), revealing the degree of evolutionary conservation of each core sequence (Casper et al., 2018). For example, as shown in Figure 3D, we identified a region at −42 kb from the TSS of Mgp, a bona fide chondrogenic marker (Barone et al., 1991), which was accessible in chondrocytes and attachment cells, whereas in tenocytes this site was closed. Another example is Sox9, a bona fide chondrogenic marker. At +303 kb from Sox9, we identified a region that was accessible in attachment cells and chondrocytes, but not in tenocytes. The same pattern was observed for a region at −330 kb from the TSS of a third bona fide chondrogenic marker, namely Col11a1 (Li et al., 1995). The opposite pattern was observed at −17 kb from the TSS of Col1a2, a bona fide tenogenic marker, where we identified a region that was accessible in attachment cells and tenocytes, but not in chondrocytes. Similar results were obtained for additional chondrogenic markers, such as Sox6, and for tenogenic markers Tnc and Col1a1 (data not shown). Importantly, we found that the chromatin accessibility patterns of these putative enhancers were in agreement with the transcriptomic and ISH results, as shown, for example, by Sox9 and Col5a1 (Figure 1C,D, Figure 1—figure supplement 4). This suggests that the mechanism for the activation of a mixed transcriptome in attachment cells is based on sharing regulatory elements with chondrocytes or tenocytes.

In vivo analysis of enhancers identifies shared domains of activity between attachment cells and neighboring tissues.

Transgenic mouse reporter enhancer assay (lacZ) of elements positive at E14.5 (marked in light blue; for each enhancer, an E14.5 whole-mount embryo, magnification of the limb and forelimb and/or hindlimb sagittal sections are shown). Left to right: Col1a1 element (mm1995) activity is seen at the teres major insertion at the scapula (n=2/9). Klf2 element (mm1988) activity is seen in hypertrophic chondrocytes and perichondrium at the humerus and forelimb digits (n=3/3). Sox9 element (mm1989) activity is seen in hypertrophic chondrocytes of the humerus (n=4/7). Mgp element (mm1990) activity is seen in the hip, digit, and metacarpals joints in addition to the posterior distal side of the femur (n=3/3). Col11a1 element (mm1991) activity is seen in the greater tuberosity insertion and anterior distal side of the femur (n=11/11). T, tendon; C, cartilage; A, attachment; HC, hypertrophic chondrocytes.

Table 1
Criteria that were used to select candidates for the in vivo enhancer activity assay.

We first compared the ATAC-seq data to bulk MARS-seq transcriptome analysis and chose peaks that were assigned to genes with differential expression (using GREAT). According to the ATAC-seq results, the TSS of most of these genes was accessible in all three tissue types. Therefore, we searched for putative enhancers that could regulate the differential gene expression. We then compared the elements overlap with E13.5-E14.5 histone marks (H3K27Ac or H3K4me1 of ENCODE) and HiC results (Andrey et al., 2017), to increase the probability of in vivo verification, in addition to preliminary results of gene expression (i.e. the chosen elements were assigned to genes of interest, which were identified by RNA-seq and validated by ISH). The degree of evolutionary conservation of the core sequence was also taken into consideration while prioritizing the elements.

mm1995
Col1a1
mm1988
Klf2
mm1989
Sox9
mm1990
Mgp
mm1991
Col11a1
mm1992, mm1994, Elnmm1986
Igfbp5
Distance from TSS−63,010−30,343303,066−42,288−331,264−3,433
−43,487
−3,605
CoordinateChr11:94872064–94874364Chr8: 72287698–72291754Chr11:113084290–113086290Chr6:136917343–136918843Chr3: 113698476–113700076Chr5: 134749924–134751424
Chr5: 134790228–134791228
Chr1: 72877526–72879452
Encode
H3K27ac/H3K4me1
+++++++
HiC+
Core seq conservationOpossumChickenPlatypusOpossumLizardOpossumChicken
Klf2 binding site+++--

Shared regulatory elements drive expression in attachment cells and flanking cartilage or tendon cells

The identification of multiple predicted enhancer regions near genes expressed by attachment cells and tenocytes or chondrocytes suggests that attachment cells are regulated predominantly by enhancers with shared activities. To test this hypothesis in vivo, we took advantage of a recently developed, site-directed transgenic mouse enhancer reporter system (Kvon et al., 2020). Using this system, we studied the activity of eight elements, which were selected because they were associated with bona fide marker genes for tenocytes or chondrocytes, and were found to be expressed in the attachment cells. Moreover, they were predicted to drive transcription in attachment cells and one of the flanking tissues (Figure 3D). The activity of these representative elements was examined at E14.5, a stage at which chondrocytes, tenocytes, and attachment cells have already been established.

Five elements were found to be active in the mouse forelimbs, as well as in other anatomical areas (Figure 4). The Col1a1-associated element drove lacZ expression at the teres major insertion into the scapula, in agreement with the ATAC-seq results, which predicted its activity in tenocytes and attachment cells (n = 2/9). This result suggests that Col1a1 element is active in tenocytes and attachment cells. Klf2 and Sox9 elements were predicted to be active in chondrocytes and attachment cells (Figure 3D). The Klf2 element indeed drove reporter activity in hypertrophic chondrocytes and perichondrium at the humerus and forelimb digits, in addition to the skull and mandible (Figure 4, n = 3/3). The Sox9 element showed activity solely in hypertrophic chondrocytes of the humerus (n = 4/7). These results suggest a chondrocyte-specific function of these two enhancers. Mgp element was also predicted to be active in chondrocytes and attachment cells (Figure 3D). Its activity was seen in forelimb and hindlimb, specifically in hip, metacarpal joints and digits as well as in the posterior distal side of the femur, a site where ligaments (e.g. the cruciate ligaments) are inserted into the femur at the knee area and at ligament insertion into to the hip (e.g. iliofemoral ligament; Figure 4, n = 3/3), verifying its activity in chondrocytes and attachment cells. Lastly, Col11a1 element activity was predicted in chondrocytes and attachment cells (Figure 3D). Its activity verified the bioinformatic analysis, showing LacZ staining in the greater tuberosity insertion, as well as in the posterior side of the skull, the nasal bone area and the anterior distal side of the femur (Figure 4, n = 11/11).

These results suggest that Mgp and Col11a1 elements are active in chondrocytes and attachment cells, whereas Col1a1 element is active in tenocytes and attachment cells, as the chromatin analysis predicts. Overall, these results provide a proof of concept for the ability of the accessible intergenic elements we have identified to act as enhancers that drive expression in both attachment cells and chondrocytes or tenocytes. This supports our hypothesis that shared enhancers activate a mixed transcriptome in attachment cells.

Krüppel-like factors are regulators of attachment cell development

Our finding of enhancers that can drive the transcription of the mixed transcriptome of the attachment cells raised the question of the identity of the transcription factors (TFs) that can potentially bind to these elements. To identify such factors, we used Genomatix to analyze accessible elements that were associated with differentially expressed genes for over-representation of transcription-factor-binding sites (TFBS), selecting the top 50 TFBS families, and then mining our transcriptomic data for the expression of these TFBS families (see Materials and methods). Among the differentially expressed genes at E14.5 we identified NFIs (Nfia), GLIs (Gli1), KLFs (Klf2 and Klf4), ZBTBs (Zbtb48) and RUNXs (Runx3; Table 2, Figure 5A), whose expression was upregulated in attachment cells. Further support for these results was provided by HOMER motif analysis (Heinz et al., 2010), which showed significant over-representation of KLFs and RUNXs TFBSs. We therefore sought to explore the possible role of KLFs as regulators of attachment cells.

Krüppel-like factors (KLFs) are regulators of attachment cell development.

(A) Heatmap of selected transcription factors at E14.5. Transcriptome analysis shows upregulated expression of Klf2, Klf4, and Gli1 in attachment cells. (B) Left: E14.5 ISH validated these results, showing Klf2 expression in attachment cells (X20 magnification, scale bar: 25 µm). Right: Scheme of attachment site. (C) Single-cell RNA-seq analysis of E13.5 Sox9+Scx+ attachment cells shows Klf2 and Klf4 expression in the five cell populations. (D) KLF2 and KLF4 are regulators of attachment cell development. a-d. Histological transverse sections through the humeral deltoid tuberosity of E15.5 Prx1-Klf2-Klf4 and E18.5 Prx1-Klf2-Klf4 mutant and control embryos (×10 and ×5 magnification, scale bar: 100 µm). a’-d’. Higher magnification of upper panel (×40 for E15.5 and ×20 for E18.5, scale bar: 40 µm). e-p. ISH for Col1a1 and Bgn genes of E15.5 Prx1-Klf2-Klf4 mutant and control embryos. ISH for Gli1, Col1a1, and Bsp genes of E18.5 Prx1-Klf2-Klf4 mutant and control embryos (×20 magnification, scale bar: 40 µm).

Table 2
Genomatix analysis of the genomic regions of cis-regulatory elements identified by ATAC-seq.

Over-representation of transcription factor binding site (TFBS) families was identified. Crossing these results with E14.5 transcriptome revealed differentially expressed TFs from the KLF (Klf4 and Klf2), GLI (Gli11), NFI (Nfia), ZBTBs (Zbtb48), and RUNX families.

TFBSGene familyDescriptionZinc finger
V$E2FFMycE2F-myc activator/cell cycle regulatorno
V$KLFSKLFKrüppel-like transcription factorsyes
V$ZF02ZbtbC2H2 zinc finger transcription factorsyes
V$NF1FNfiNuclear factor 1yes
V$MAZFMaz, Patz1Myc-associated zinc fingersyes
V$GLIFGLIGLI zinc finger familyyes
V$PLAGPlag1, Plag2, Plagl1Pleomorphic adenoma geneyes
V$HAMLRunxAcute myelogenous leukemia factorsno
V$NFKBHIVEPNuclear factor kappa B/c-relyes
V$SMADSmadVertebrate SMAD family of transcription factorsno

Focusing on Klf2, which was found to be differentially expressed by the bulk transcriptome analysis of the attachment site (Figure 1—figure supplement 2, Figure 1—figure supplement 3), we first validated its expression in the forming attachment site by ISH (Figure 5B). Next, we analyzed the enhancers that were shown by the bioinformatic analysis to be active in attachment cells and either tenocytes or chondrocytes (Figure 4). Three of these enhancers had Klf2-binding sites in their sequence (Table 1), further supporting a potential role for KLFs during attachment site development.

Previous studies demonstrated that Klf2 and Klf4 are functionally redundant, as KLF4 has ~90% sequence similarity to KLF2 in its zinc finger DNA-binding domain, suggesting that these factors could have common target sequences (Chiplunkar et al., 2013). In our scRNA-seq analysis, the expression of both genes was found in Sox9+/Scx+ attachment progenitors (Figure 5C). We therefore proceeded to study attachment cell development upon blocking the expression of both Klf2 and Klf4 in limb mesenchyme, using Prx1-Cre as a deleter and focusing on E15.5 and E18.5, a period during which the attachment site of the deltoid tuberosity undergoes differentiation and consequently grows in size. Transverse histological sections through the deltoid tuberosity of E15.5 control mice showed that the attachment cells were packed together and surrounded by ECM (Figure 5Da’). In contrast, in the putative attachment site of Prx1-Klf2-Klf4 double conditional knockout (dcKO) embryos, the cells were sparse with reduced ECM (Figure 5Db’). By E18.5, this difference was more pronounced, as attachment cells failed to differentiate (Figure 5Dc’,d’). To gain a molecular understanding, we studied the expression of several genes that were previously shown to be expressed at these stages in the attachment site (Felsenthal et al., 2018; Figure 1D). Indeed, we found that the expression of Col1a1, Gli1, Bsp, Bgn, and Col5a1 was reduced in the dcKO attachment site, relative to the control (Figure 5De-p), supporting a role for KLF2/4 in attachment cells differentiation.

Finally, to further validate the involvement of KLFs in activation of gene expression in the attachment site, we searched for KLF2/4 binding sites in ATAC-seq peaks associated with the 374 genes that were shown to be expressed by attachment cells (Figure 1C). Interestingly, we found that many of these genes had KLF2/4 binding sites in their regulatory regions (72% of the 374 attachment genes relative to 53% in the whole genome, p<1e-4, chi-square). We then searched for KLF2/4 binding sites in ATAC-seq peaks that were associated with genes whose expression was reduced in the dcKO attachment site (Figure 5De-p). For Gli1, we found KLF2/4 TFBSs in peaks that reside −2.1 and −1.5 kb from its TSS (Table 3). For Col5a1, we found multiple binding sites for KLF2 or KLF4. Together, these results indicate that KLF2/4 play an essential role in regulating attachment cell gene expression.

Table 3
Gli1 and Col5a1 genes are expressed in the attachment site and have TFBS for Krüppel-like factors (KLFs).

Search for KLF2/4 binding sites in ATAC-seq peaks that were associated with genes whose expression was reduced in the dcKO attachment site (Figure 5De-p). For Gli1, we found KLF2/4 TFBSs in peaks that reside −2.1 and −1.5 kb from its TSS (indicated as + under TFBS). For Col5a1, we found multiple binding sites for KLF2 or KLF4 (indicated as + under TFBS).

HOMERGREATENCODETFBS
PeakGeneDistance to TSSAnnotationgene1Distance gene 1gene2Distance gene 2H3K27acH3K4me1H3K4me3KLF2
(Genomatix)
KLF4
(Genomatix)
KLF4
(Homer)
chr10_127339817_127340317Gli11512IntronGli11522Arhgap916340++++++
chr10_127339170_127339670Gli12159IntronGli12169Arhgap915693++++
chr10_127338776_127339276Gli12553IntronGli12563Arhgap915299+++
chr10_127341889_127342389Gli1−560TSSGli1−550NANA+++
chr10_127342228_127342728Gli1−899TSSGli1−889NANA+++
chr2_27740165_27740665Rxra29832IntronCol5a1−146010Rxra63214+++
chr2_27745961_27746569Rxra35682IntronCol5a1−140160Rxra69064+++
chr2_27934260_27934760Col5a148085IntronCol5a148085Fcnb150375+++
chr2_27886819_27887319Col5a1644IntronCol5a1644NANA+++
chr2_27802472_27802972Col5a1−83703IntergenicCol5a1−83703Rxra125521++
chr2_27848386_27848886Col5a1−37789IntergenicCol5a1−37789Rxra171435+++
chr2_27885904_27886561Col5a1−192TSSCol5a1−192NANA+++
chr2_27813321_27813821Col5a1−72854IntergenicCol5a1−72854Rxra136370++

Combined with the bioinformatic analysis of chromatin and transcriptomic data, these results suggest that KLF2/4 are major regulators of tendon-to-bone attachment, playing a central role in attachment cell differentiation.

Discussion

In this work, we describe the unique transcriptome that allows cells of the attachment between tendon and bone to act as a transitional tissue. The ability to activate a combination of chondrogenic and tenogenic transcriptomes is regulated by sharing enhancers with these cells. Finally, we identify the transcription factors KLF2/4 as regulators of these unique bi-fated cells.

The existence of borders between tissues that differ in cell type, extracellular matrix composition, structure, and function raises the question of how tissues are connected. The border can be sharp, as seen in blood vessels, where pericytes and endothelial cells are separated by a basement membrane, or between the esophagus and the stomach in the gastrointestinal tract (Bergers and Song, 2005; San Roman and Shivdasani, 2011). On the other hand, the border can be less defined histologically and molecularly, thus forming a transitional tissue. Examples for the latter are the borders between the sections of the small intestine and between tendon and bone (Liu et al., 2007; Romih et al., 2005; San Roman and Shivdasani, 2011). From a broader perspective, as all organs and systems are made of different tissues, understanding the biology of border tissues is imperative. Moreover, some of these border tissues are involved in various pathologies. For example, gastric cancers may emerge from distinct anatomical areas, such as the esophagus–stomach boundary (Chawengsaksophak, 2019; McDonald et al., 2015; San Roman and Shivdasani, 2011).

In the case of the attachment between tendon and bone, the significance of this tissue is demonstrated by enthesopathies, a collective name for injuries and pathologies of the enthesis. For example, over 30% of the population over the age of 60 will injure their shoulder’s rotator cuff (Lehman et al., 1995). Failure rates of surgical reattachment range from 20% for small tears to 94% for repair of massive tears (Galatz et al., 2004; Harryman et al., 1991). The high failure and recurrence rates of these procedures highlight the need for understanding the biology of this complex transitional tissue of the enthesis. This understanding may allow the development of new strategies to improve the treatment of enthesopathies.

There are two options to form a transitional tissue. The first strategy is by mixing cells from the two neighboring tissues, such as in the epithelia of the urinary tract (Romih et al., 2005). Alternatively, the border cells can express a mixture of the transcriptomes of the two neighboring cell types. As we show here, the attachment cells represent the latter strategy well, as they express a high number of genes that are differentially expressed by either tenocytes or chondrocytes. These cells display morphological features that are shared with tenocytes and chondrocytes (Thomopoulos et al., 2010). Our results therefore provide a molecular explanation for the age-old histological definition of enthesis cells as fibrocartilage, which was based on their morphology (Thomopoulos et al., 2010). Moreover, the finding of mixed matrix genes in the transcriptome of the attachment cells may provide a mechanism for the formation of a transitional tissue, which allows safe transfer of forces by the tendon between muscle and bone.

In addition to expression of chondrogenic and tenogenic genes, we identified genes that are uniquely expressed by attachment cells. These genes may provide another level of specificity to the regulation of the development of this unique tissue. Finally, our scRNA- seq of the attachment cells revealed heterogeneity, which was a consequence of varying levels of chondrogenic and tenogenic gene expression. This heterogeneity may represent the differentiation processes that the Sox9/Scx-positive progenitors undergo during development, ultimately leading to their terminal cell fate.

Our finding that attachment cells are bi-fated raises the question of the mechanism that underlies this fate. An immediate implication of our finding is that there must be an epigenetic mechanism that supports the bi-fated state. The observed chromatin accessibility at the sites of the promoters of most of the shared genes in all three cell types rules out the possibility of limited promoter accessibility as the main mechanism. By contrast, the high percentage of shared accessible intergenic sites between attachment cells and one group of flanking cells, that is chondrocytes or tenocytes, suggests that this is the main mechanism. Moreover, many of these shared sites correlated with putative enhancers in ENCODE datasets. Finally, we verified the ability of three different enhancers from that list to drive gene expression in attachment cells and in either tendon or cartilage. These findings strongly support our hypothesis that the regulatory mechanism is based on the ability of attachment cells to share enhancers with either chondrocytes or tenocytes in order to drive the mixed expression profile of these bi-fated cells. It is important to emphasize that we have also identified a group of intergenic elements that were accessible exclusively in attachment cells. These attachment-specific elements may act as enhancers that drive the expression of genes that are specific to attachment cells. It is, of course, possible that they participate in the regulation of shared genes as well. Overall, both sets of putative enhancers, namely shared and attachment cell-specific, may play important roles in the genetic program that regulates the development of this unique tissue.

Sharing enhancers is not the only possible strategy for the generation of a mixed transcriptome. A simple alternative would be a specific set of enhancers to be used by the attachment cells. A possible explanation for the sharing strategy is the common origin of all these cells, which is limb mesenchyme originating from lateral plate mesoderm (Johnson and Tabin, 1997). It is possible that during development, limb mesenchymal progenitors display highly accessible chromatin; yet, during differentiation, this accessibility is restricted to prevent the expression of genes from alternate lineages. In contrast to this restriction process, in the bi-fated attachment cells the shared sites are maintained accessible to allow the expression of the mixed transcriptome. A mechanism for silencing of genes of alternate lineages was previously described. For example, polycomb-repressed chromatin leads to silencing of genes of alternate lineages, leading to the commitment to a specific cell fate (Aldiri et al., 2017; Laugesen and Helin, 2014; Zhu et al., 2013). Interestingly, previous studies demonstrated the importance of chromatin repression in the developing limb, showing how deletion of Ezh2, which acts as the enzymatically active subunit of PRC2, leads to skeletal malformations (Deimling et al., 2018). This obviously raises the question of the mechanism that prevents this silencing in attachment cells.

It is clear that we cannot exclude the possibility that an active mechanism, such as the SWI/SNF remodeling complexes, opens the chromatin structure in bi-fated cells to allow attachment cell dual behavior (Hu et al., 2011). However, such a mechanism cannot explain why the strategy of shared enhancers was selected. Finally, a mixed transcriptome is one of the hallmarks of stem cell pluripotency, as progenitor cells eventually chose one fate over the other upon differentiation (Johnson et al., 2015; Soldatov et al., 2019). Our findings suggest a different strategy, where using both programs facilitates the establishment of a specific attachment tissue. Overall, our results reveal a novel function for chromatin state, which allows the activation of two sets of genes in a third cell type to create a new cell fate that forms a transitional tissue.

KLF2 and KLF4 are known to regulate several biological processes, such as promoting the differentiation of gut and skin (KLF4, [Katz et al., 2002; Segre et al., 1999]) as well as the immune system (KLF2, [Kuo et al., 1997]), maintaining pluripotency of embryonic stem cells (KLF2 and KLF4, [Jiang et al., 2008]), and, together with other factors, inducing pluripotency to generate iPSC by reprogramming (KLF4, [Takahashi and Yamanaka, 2006]). Several works describe the involvement of KLF2 and KLF4 in the musculoskeletal system. In bones, Klf4 over-expression in osteoblasts caused delayed bone development, in addition to impaired blood vessel invasion and osteoclast recruitment (Michikami et al., 2012). Another study showed that KLF2/4 are expressed during chick limb development in tendons and ligaments as part of the genetic program that regulates connective tissues (Orgeur et al., 2018). Previous studies showed that KLF2 and KLF4 display high similarity in protein sequences (Dang et al., 2000; Shields and Yang, 1997), suggesting that these factors could have common target sequences and may be functionally redundant. Indeed, loss of both KLF2 and KLF4 during embryogenesis led to abnormal blood vessel development and early lethality. This phenotype was more severe than what was observed in embryos that lost only KLF2 or KLF4 (Chiplunkar et al., 2013). Furthermore, previous work identified 313 target genes shared between KLF2 and KLf4, suggesting that they overlap in regulating gene expression (Orgeur et al., 2018). In this work, we show that KLF2/4 are central regulators of the attachment site. While the attachment did form initially in their absence, the subsequent differentiation failed, suggesting that KLF2/4 play a role at this stage. While we concentrated in this study on the attachment site, it is most likely that KLF2/4 play a role also in other musculoskeletal tissues such as the skeleton, tendon, and muscle. This possibility is supported by previous studies, where KLF2/4 were shown to be expressed in osteoblasts, chondrocytes, tenocytes, and muscle connective tissues (Michikami et al., 2012; Orgeur et al., 2018). Previous studies demonstrate the role of muscle-induced mechanical load in the development of attachment site (Blitz et al., 2009). In that context, our finding that KLF2/4 regulate the differentiation of attachment cells is interesting, because previous works have shown that these factors are mechanically regulated. It was shown in mice that shear stress on the vessels induced by blood flow leads to upregulation of Klf2 expression (Lee et al., 2006). Additional in vitro studies showed that KLF2 and KLF4 are influenced by shear stress (Chiplunkar et al., 2013; Dekker et al., 2005; Villarreal et al., 2010). It is therefore possible that these factors are regulated by muscle forces, leading to the proper differentiation and maturation of the attachment site.

The ability of KLF2/4 to regulate gene expression in the attachment site is supported by our finding that many of the genes that were expressed by attachment cells had in their regulatory region KLF2/4 binding sites. Yet, it is clear that not all of them share this property, suggesting that these two factors are part of a larger transcriptional network. For example, our bioinformatic analysis identified other TF families such as GLI’s, RUNX’s, and NFI’s as regulators of the attachment sites. Gli1 was previously reported as a marker for enthesis cells (Dyment et al., 2015; Felsenthal et al., 2018; Liu et al., 2012a; Schwartz et al., 2015). Since gene expression by attachment cells is regulated by sharing enhancers with chondrocytes or tenocytes, it is reasonable to assume that some regulators of these cells might be part of the network that regulates the attachment cells. Indeed, loss of the tendon regulator Scx in mice led to failure of attachment cells to differentiate. Additionally, loss of the chondrogenic regulator Sox9 in Scx-expressing cells led to failure in attachment site formation (Blitz et al., 2013; Blitz et al., 2009).

To conclude, by characterizing the transcriptome and chromatin landscape of tendon-to-bone attachment cells, we provide a molecular understanding of the bi-fated identity of these cells. Moreover, by identifying the transcription factors KLF2/4 as central regulators and the strategy of sharing enhancers with either tenocytes or chondrocytes, we provide a mechanism that regulates these bi-fated cells (Figure 6). These findings present a new concept for the formation of a border tissue, which is based on the simultaneous expression of a mixed transcriptome of the two flanking cell types by the intermediate cells. This strategy allows the formation of a unique transitional tissue without developing de novo a dedicated genetic program that regulates a third, new cell fate.

Proposed model of bi-fated tendon-to-bone attachment cells that are regulated by shared enhancers and KLF transcription factors.

Tenocytes (green, left) and chondrocytes (red, right) express tenogenic (i.e. Col1a1) or chondrogenic (i.e. Col11a1) genes, respectively, whereas attachment cells (yellow, middle) express both chondrogenic and tenogenic genes to form the attachment site. Attachment cells duality of gene expression is regulated epigenetically by intergenic chromatin areas, which are accessible in these cells and in either tenocytes or chondrocytes. Additionally, at the transcriptional level, the transcription factors KLF2/4 are expressed by attachment cells and regulate their differentiation.

Materials and methods

For Key Resources Table, see Appendix.

Animals

The generation of floxed Klf2 (Lee et al., 2006), floxed Klf4 (Katz et al., 2002), Prx1-Cre (Logan et al., 2002), Sox9-CreER (Soeda et al., 2010), Col2-CreERT (Nakamura et al., 2006), Col2a1-Cre (Ovchinnikov et al., 2000), R26R-tdTomato (Madisen et al., 2010), and Scx-GFP (Pryce et al., 2007) mice have been described previously.

To create Col2a1-CreER, R26R-tdTomato, Sox9-CreER, R26R-tdTomato, and Col2a1-Cre, R26R-tdTomato reporter mice, all on the background of Scx-GFP, floxed R26R-tdTomato mice were mated with Col2a1-CreER, Sox9-CreER or Col2a1-Cre mice, respectively. These strains were mated on a mixed background of C57BL/6 and B6.129 (ICR) mice and used for LCM and FACS experiments. To create Prx1-Klf2-Klf4 mutant mice, floxed Klf2-Klf4 mice were mated with Prx1-Klf2-Klf4 mice. As a control, we used embryos that lack Cre alleles. E14.5 wild-type C57BL/6 mice were used for ISH experiments as well.

For FACS experiments, Sox9-CreER or Col2a1-CreER mice were crossed with Rosa26-tdTomato reporter mice, all on the background of Scx-GFP. Induction of Cre recombinase was performed at various pregnancy stages by administration of 0.03 mg/gr tamoxifen/body weight in corn oil by oral gavage (stock concentration was 5 mg/ml). In all timed pregnancies, plug date was defined as E0.5. For harvesting of embryos, timed-pregnant females were sacrificed by cervical dislocation. Tail genomic DNA was used for genotyping by PCR.

Laser capture microdissection

Request a detailed protocol

E14.5 Col2a1-tdTomato-Scx-GFP mouse forelimbs were dissected, shortly fixated in 4% PFA, washed with PBS and cryo-embedded (as described by Bhattacherjee et al., 2004 Bhattacherjee et al., 2004). Next, samples were cryo-sectioned and mounted on LCM slides (PET, Zeiss), washed with RNase-free water and EtOH (Arcturus dehydration component kit) according to an altered protocol of Pazin et al., 2012. LCM (PALM MicroBeam C system, Zeiss) was calibrated for refined tissue cutting. Isolated cells were collected to LCM caps (Adhesive Cap 500 clear, Zeiss) and RNA was purified using RNeasy FFPE Kit (Qiagen). The resulting RNA was the input for the bulk MARS-seq protocol (RNA sequencing).

Preparation of single-cell suspension for fluorescence-activated cell sorting (FACS)

Request a detailed protocol

Flow cytometry analysis and sorting were performed at the Weizmann Institute of Science Flow Cytometry Core Facility on a BD FACS AriaIII instrument (BD Immunocytometry Systems) equipped with 488, 407, 561, and 633 nm lasers, using a 70 μm nozzle, controlled by BD FACS Diva software v8.0.1 (BD Biosciences). Further analysis was performed using FlowJo software v10.2 (Tree Star). For collection of cells, Sox9-CreERT2-tdTomato;ScxGFP or Col2a1-CreER-tdTomato;ScxGFP mice were crossed with Rosa26-tdTomato;ScxGFP reporter mice. Embryos were harvested at E13.5 following tamoxifen administration at E12.0, as described above. Forelimbs were dissected and suspended in cold PBS. To extract cells from tissues, PBS was replaced with 1 ml heated 0.05% trypsin and collagenase type V (dissolved in DMEM, Sigma) and incubated for 15 min at 37°C, gently agitated every 5 min. Tissues were then dissociated by vigorous pipetting using 1 ml tips. Next, 4 ml of DMEM supplemented with 10% FBS and 1% Pen-Strep was added and cell suspensions were filtered with 40 μm filter net. Finally, tubes were centrifuged at 1000 rpm for 7 min, supernatant was removed and cells were resuspended in 0.5–1 ml of cold PBS and used immediately for FACS. Single-stained GFP and tdTomato control cells were used for configuration and determining gate boundaries. Live cells were gated by size and granularity using FSC-A versus SSC-A and according to DAPI staining (1 μg/ml). FSC-W versus FSC-A was used to further distinguish single cells. In addition, unstained, GFP-stained only and tdTomato-stained only cells were mixed in various combinations to verify that the analysis excluded false-positive doublets. GFP was detected by excitation at 488 nm and collection of emission using 502 longpass (LP) and 530/30 bandpass (BP) filters. tdTomato was detected by excitation at 561 nm and collection of emission using a 582/15 BP filter. DAPI was detected by excitation at 407 nm and collection of emission using a 450/40 BP filter.

Real-time PCR (RT-PCR)

Request a detailed protocol

Total RNA was purified from LCM-isolated samples of E14.5 mouse forelimbs using RNeasy FFPE Kit (Qiagen). Reverse transcription was performed with High Capacity Reverse Transcription Kit (Applied Biosystems) according to the manufacturer's protocol. Analysis of Col2a1 and Scx was performed to monitor RNA quality during LCM calibrations, whereas RNA quantity was monitored by analysis of β–actin. RT-PCR was performed using Fast SYBR Green master mix (Applied Biosystems) on the StepOnePlus machine (Applied Biosystems). Values were calculated using the StepOne software. Data were normalized to 18S rRNA or β-actin in all cases.

In situ hybridization

Request a detailed protocol

Section ISH were performed as described previously (Riddle et al., 1993). Single-and double-fluorescent ISH on paraffin sections were performed using DIG- and/or FITC-labeled probes (listed in Table 4; Shwartz and Zelzer, 2014). After hybridization, slides were washed, quenched, and blocked. Probes were detected by incubation with anti-DIG-POD (Roche; 1:300) and anti-FITC-POD (Roche, 1:200), followed by Cy2-tyramide- and Cy3-tyramide-labeled fluorescent dyes according to the instructions of the TSA Plus Fluorescent Systems Kit (Perkin Elmer).

Table 4
List of probes used for in situ hybridization.
Probe nameGenomic positionRefseq templateSize (bp)
Igfbp5792–1264NM_010518.2502
Col11a1831–1224NM_007729.3393
Biglycan147–716NM_007542.4569
Wwp22488–3056NM_025830.3568
Col5a1626–1298NM_015734.2672
Col1a14295–4475NM_007742.4180
Klf2895–1512NM_008452.2617
Gli11810–2447NM_010296.2638
Bsp145–1058NM_008318.31955
Scx273–1129NM_198885.3856
Mfap4276–966NM_029568.2690
Ptn1247–1932NM_008973.2685
Col3a1707–1388NM_009930.2681
Mgp69–553NM_008597.4485
Eln154–693NM_007925.4540
Mmp14708–1283NM_008608.3575
Col2a14474–4879NM_001113515.2406
Col9a12553–3100NM_007740.3547
Sox950–797NM_011448.4748
Snorc32–437NM_028473.1405

Single-molecule fluorescent in situ hybridization (smFISH)

Request a detailed protocol

Harvested E14.5 forelimbs were fixed with cold 4% formaldehyde (FA) in PBS and incubated first in 4% FA/PBS for 3 hr, then in 30% sucrose in 4% FA/PBS overnight at 4°C with constant agitation. Fixed tissues were embedded in OCT and sectioned at a thickness of 10 µm. The preparation of the probe library, hybridization procedure, and imaging conditions were previously described (Itzkovitz et al., 2012; Lyubimova et al., 2013; Raj et al., 2008). In brief, probe libraries were designed against biglycan (Bgn) and Wwp2 mRNA sequences using the Stellaris FISH Probe Designer (Biosearch Technologies, Inc, Petaluma, CA) coupled to Quasar 670 and CAL Fluor Red 610, respectively. Libraries consisted of 17–96 probes each of length 20 bps, complementary to the coding sequence of each gene (Supplementary file 2). Nuclei were stained with DAPI. To detect cell borders, Alexa Fluor 488 conjugated phalloidin (Thermo Fisher, A12379) was added to the GLOX buffer, which was wash for 15 min. Slides were mounted using ProLong Gold (Molecular Probes, P36934).

Image acquisition and analysis

Request a detailed protocol

For smFISH image acquisition, we used a Nikon-Ti-E inverted fluorescence microscope equipped with a Photometrics Pixis 1024 CCD camera to image 10-µm-thick cryosections. For image analysis, we used ImageM, a custom MATLAB program (Lyubimova et al., 2013), which was used to compute single-cell mRNA concentrations by segmenting each cell manually according to the cell borders and the nucleus. The size of the nucleus was detected automatically by the program according to the DAPI signal. For each cell, the concentration of cytoplasmic mRNA of each gene was calculated by measuring the number of dots per volume. Images were visualized and processed using ImageJ 1.51 hr (Schindelin et al., 2012) and Adobe Illustrator CC2018.

Bulk RNA sequencing

Request a detailed protocol

For this analysis, we adapted the MARS-seq protocol (Jaitin et al., 2014; Keren-Shaul et al., 2019) to generate RNA-seq libraries for expression profiling of the purified RNA from E14.5 LCM-isolated samples. Briefly, RNA from each sample was barcoded during reverse transcription and pooled. Following Agencourt Ampure XP beads cleanup (Beckman Coulter), the pooled samples underwent second strand synthesis and were linearly amplified by T7 in vitro transcription. The resulting RNA was fragmented and converted into a sequencing-ready library by tagging the samples with Illumina sequences during ligation, RT, and PCR. Libraries were quantified by Qubit and TapeStation as well as by qPCR for actb gene as previously described (Jaitin et al., 2014; Keren-Shaul et al., 2019). Sequencing was done on a Hiseq 2500 SR50 cycles kit (Illumina).

The data were analyzed using the Pipeline Pilot-designed pipeline for transSeq (by INCPM, https://incpmpm.atlassian.net/wiki/spaces/PUB/pages/36405284/tranSeq+on+Pipeline-Pilot). Briefly, the analysis included adapter trimming, mapping to the mm9 genome, collapsing of reads with the same unique molecular identifiers (UMI) of 4 bases (R2) and counting of the number of reads per gene with HTseq-count (Anders et al., 2015), using the most 3’ 1000 bp of each RefSeq transcript. DESeq2 (version 1.4.5, Love et al., 2014) was used for normalization and differential expression analysis with betaPrior set to true, cooksCutoff=FALSE, independentFiltering=FALSE. Benjamini-Hochberg method was used to adjust the raw p-values for multiple testing. Genes with adjusted p-value≤0.05 and fold change ≥ 2 between every two conditions were considered as differential. PCA analysis was done using log-transformed normalized data (DESeq2 function rlogTransformation with parameter blind=TRUE) with a modified plotPCA function of DESeq2. Clustering of the log-normalized read count of differentially expressed genes was done using CLICK algorithm (Expander package version 7.1, Ulitsky et al., 2010), followed by visualization by R (R Development Core Team, 2013). Further analysis was performed using GSEA (Broad institute) and Gorilla (Eden et al., 2009; Subramanian et al., 2005). The crude data of the work has been deposited on NCBI GEO (GSE144306).

Single-cell library preparation using chromium 10x genomics platform

Request a detailed protocol

Cells were counted and diluted to a final concentration in PBS supplemented with 0.04% BSA. Cellular suspension was loaded onto Next GEM Chip G targeting 4000 cells and then ran on a Chromium Controller instrument to generate GEM emulsion (10x Genomics). Single-cell 3’ RNA-seq libraries were generated according to the manufacturer’s protocol (10x Genomics Chromium Single Cell 3’ Reagent Kit User Guide v3/v3.1 Chemistry).

Next-generation sequencing of single-cell libraries

Request a detailed protocol

Single-cell 3’ RNA-seq libraries were quantified using NEBNext Library Quant Kit for Illumina (NEB) and high sensitivity D1000 TapeStation (Agilent). Libraries were pooled according to targeted cell number, aiming for ~50,000,000 reads per cell. Pooled libraries were sequenced on a NovaSeq 6000 instrument using an SP 100 cycles reagent kit (Illumina) (R1 28 bases, R2 82 bases, and I1 eight bases), specifically 290M fragments of the relevant library were sequenced.

scRNA-seq bioinformatic analysis

Request a detailed protocol

Demultiplexing and alignment were performed using cellranger (10x Genomics version 3.0.2) bioinformatics pipeline using mm10 genome, followed by detecting swapped barcodes between libraries, since the NovaSeq run contained a mix of three libraries. For this, we used R (3.6.3) and the package DropletUtiles (Lun et al., 2019) (1.6.1) with the function swappedDrops (parameter min.frac = 0.9) and the molecule_info.h5 input. 1.5% of the molecules were detected as swapped. Detecting empty droplets was done with the function emptyDrops (parameter lower = 300), the number of cells detected was 2213.

For additional analysis, the R package Seurat (3.1.5) was used (Butler et al., 2018). The analysis included filtering genes (must be expressed in at least three cells), filtering cells with over 20% expression of mitochondria genes and high and low 3% percentiles of total number of genes and UMIs; this resulted in 1218 cells. In addition, one of the clusters had significantly lower RNA counts as compare to other collected cells, which might have indicated that it contained mainly low-quality cells. Indeed, examination of the 10 most highly expressed genes using heatmap demonstrated this; therefore, this cluster was removed from further analysis. A total of 1076 cells were clustered using 2000 variable genes and 15 principal components (PCs) (resolution = 0.4). UMAP plot was made with the Seurat functions RunUMAP and DimPlot and violin plots were made with the VpnPlot function. The crude data of the work have been deposited on NCBI GEO (GSE160090).

Assay for transposase-accessible chromatin with high-throughput sequencing

Request a detailed protocol

ATAC-seq data were trimmed from their adaptors and filtered from low quality reads using Cutadapt followed by alignment to the mm10 genome (GRCm38.p5) using Bowtie2 (version 2.3.4.1) (Langdon, 2015). PCR-duplicate reads were removed with Picard ‘MarkDuplicates’ (http://broadinstitute.github.io/picard/). Mitochondrial reads were removed from the alignment, and the data were further filtered to contain only reads with a unique mapping with SAMtools (-F 4 f 0 × 2). Read pairs with inner distance of up to 120 bp were selected as representing the accessible chromatin region. MACS2 (version 2.1.1.20160309) (Zhang et al., 2008) was applied for peak calling using the setting: callpeak -f BAMPE--nomodel. Peaks from all samples were combined and merged with BEDTools (Quinlan, 2014), followed by extension to a minimum length of 500 bp. For every tissue, a set of reproducible peaks was obtained by voting, which means that a normalized read count ≥ 30 was detected in at least 50% of the replicates. Peaks that were not reproducible in any tissue were removed. Peaks that reside in the ENCODE ‘Blacklist’ regions, that is regions that were previously found by ENCODE (ENCODE Project Consortium, 2012) to produce artificial signal (http://mitra.stanford.edu/kundaje/akundaje/release/blacklists/), were also eliminated. Peak quantification was done with BedTools Quinlan, 2014 following by DESeq2 (Love et al., 2014) normalization. Peaks with an averaged normalized read count ≥ 30 in at least one of the studied tissues were selected for the downstream analyses. The crude data of the work have been deposited on NCBI GEO (GSE144306).

Annotation and genomic feature enrichment analysis

Request a detailed protocol

Annotation of ATAC-seq peaks was performed using HOMER (Heinz et al., 2010) and GREAT (McLean et al., 2010). When a peak was associated by GREAT to multiple genes, the two closest genes were selected for further analysis. ATAC-seq peaks that were at a distance of up to −2 kb down or +0.5 kb up from a TSS of their annotated gene (HOMER) were considered as promoter peaks; otherwise, peaks were considered as distal. To rank distal ATAC-seq peaks as putative cis-regulatory elements, we calculated the overlap between the peaks and relevant histone modification datasets (ChIP-seq) performed by the ENCODE project (Yue et al., 2014) on E13.5 C57BL/6 mouse embryo limb. The overlap was calculated using BEDTools intersect (Quinlan, 2014). The following datasets were used: ENCSR905FFU (H3K27ac) and ENCSR426EZM (H3K4me1) as markers of enhancers, ENCSR416OYH (H3K4me3) as a marker of promoters and ENCSR022DED (H3K9me3). Overlap with the phastConsElements60wayPlacental track downloaded from the UCSC site (Casper et al., 2018) was calculated to account for evolutionary conservation. Enrichment analysis of over-representation of TFBSs in the ATAC-seq peaks was performed with the RegionMiner tool of Genomatix and HOMER.

Enhancer reporter assays in mouse embryos

Request a detailed protocol

Candidate enhancers were PCR-amplified and cloned upstream of a Shh-promoter-LacZ-reporter cassette. We used a mouse enhancer-reporter assay that relies on site-specific integration of a transgene into the mouse genome (Kvon et al., 2020). In this assay, the reporter cassette is flanked by homology arms targeting the safe harbor locus (Tasic et al., 2011). Cas9 protein and a sgRNA targeting H11 were co-injected into the pronucleus of FVB single cell-stage mouse embryos (E0.5) together with the reporter vector (Kvon et al., 2020). Embryos were sampled and stained at E14.5. Embryos were excluded from further analysis if they did not carry the reporter transgene. All mouse works were reviewed and approved by the Lawrence Berkeley National Laboratory Animal Welfare and Research Committee.

Appendix 1

Appendix 1—key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Strain, strain background (Mus musculus)Sox9-CreER, C57BL/6DOI:10.1002/dvg.20667FACS
Strain, strain background (M. musculus)R26R-tdTomato, C57BL/6DOI:10.1038/nn.2467FACS
Strain, strain background (M. musculus)Scx-GFP, C57BL/6DOI:10.1002/dvdy.21179RRID:MGI:3717422FACS
Strain, strain background (M. musculus)Col2a1-Cre, C57BL/6 and B6.129 (ICR)PMID:10686612Used as mixed background (C57BL/6 and B6.129 (ICR)) for LCM
Strain, strain background (M. musculus)Col2-CreERT,
C57BL/6 and B6.129 (ICR)
DOI:10.1002/dvdy.20892RRID:IMSR_JAX:006774Used as mixed background (C57BL/6 and B6.129 (ICR)) for FACS
Strain, strain background (M. musculus)Floxed Klf2, C57BL/6Eric SebzdaDOI:10.1016/j.devcel.2006.09.006
Strain, strain background (M. musculus)Floxed Klf4, C57BL/6MMRRC*RRID:MMRRC_029877-MUPMID:12015290 *Mutant Mouse Regional Resource Center at UC Davis
Strain, strain background (M. musculus)Prx1-Cre, C57BL/6DOI:10.1002/gene.10092
Strain, strain background (M. musculus)C57BL/6The Jackson LaboratoryIn situ hybridization and single-molecule fluorescent in situ hybridization
Strain, strain background (M. musculus)FVB/Shh-ZRSem7Axvi (396C>T variant knock-in), FVBDOI:10.1016/j.cell.2020.02.031Enhancer reporter assays in mouse embryos
Strain, strain background (M. musculus)FVBCharles Riverhttps://www.criver.com/Enhancer reporter assays in mouse embryos
Antibodyanti-DIG-POD (sheep polyclonal)RocheCat# 11207733910ISH (1:300), DOI:10.1007/978-1-62703-989-5_15
Antibodyanti-FITC-POD (sheep polyclonal)RocheCat# 11426346910ISH (1:200),
DOI:10.1007/978-1-62703-989-5_15
Commercial assay or kitTSA Plus Fluorescent Systems KitPerkin ElmerCat# NEL753001KTISH (Cy3 + fluorescein), DOI:10.1007/978-1-62703-989-5_15
AntibodyAnti-SOX9 (rabbit polyclonal)MilliporeAB5535(1:200)
Commercial assay or kitHistogene LCM Frozen Section Staining KitThermoFisher ScientificCat# KIT0401LCM
Commercial assay or kitMembraneSlide 1.0 PETCarl Zeiss MicroscopyCat# 415190-9051-000LCM
Commercial assay or kitAdhesiveCap 500 clearCarl Zeiss MicroscopyCat# 415190-9211-000LCM
Commercial assay or kitRNeasy FFPE KitQiagenCat# 73504LCM
Sequence-based reagentsmFISH probesStellaris FISH Probe Designer (Biosearch Technologies, Inc, Petaluma, CA)See Supplementary file 2
Software, algorithmPipeline Pilot-designed pipeline for transSeqINCPMhttps://incpmpm.atlassian.net/wiki/spaces/PUB/pages/36405284/tranSeq+on+Pipeline-PilotBulk RNA sequencing
Software, algorithmHTseq-countDOI: 10.1093/bioinformatics/btu638
Software, algorithmDESeq2DOI:10.1186/s13059-014-0550-8Version 1.4.5
Software, algorithmExpander packageDOI: 10.1038/nprot.2009.230Version 7.1, CLICK algorithm
Software, algorithmcellranger10x GenomicsVersion 3.0.2scRNA-seq bioinformatic analysis
Software, algorithmDropletUtilesDOI:10.1186/s13059-019-1662-yVersion 1.6.1scRNA-seq bioinformatic analysis
Software, algorithmSeuratDOI:10.1038/nbt.4096Version 3.1.5scRNA-seq bioinformatic analysis
OtherBulk RNA sequencing datasetThis paperNCBI GEO (GSE144306)
OtherSingle-cell RNA sequencing datasetThis paperNCBI GEO (GSE160090)
Commercial assay or kitHiseq 2500 SR50 cycles kitIlluminaBulk RNA sequencing
Commercial assay or kit10x Genomics Chromium Single Cell 3’ Reagent Kit10x GenomicsUser Guide v3/v3.1 Chemistry
Commercial assay or kitNEBNext Library Quant Kit for IlluminaNEB
Commercial assay or kitSP 100 cycles reagent kitIllumina
OtherAlexa Fluor 488 conjugated phalloidinThermoFisher ScientificA12379smFISH
OtherProLong GoldMolecular ProbesP36934smFISH
Commercial assay or kitNextSeq 500 High Output v2 Kit (75 cycles)IlluminaFC-404–2005
Commercial assay or kitNextera Index Kit 24ind, 96smpIlluminaFC-121–1011
Commercial assay or kitNextera DNA Sample Prep Kit (24 sam)IlluminaFC-121–1030
Software, algorithmBowtie2DOI:10.1186/s13040-014-0034-0Version 2.3.4.1
Software, algorithmPicard ‘MarkDuplicates’http://broadinstitute.github.io/picard/Version 1.119
Software, algorithmMACS2DOI:10.1186/gb-2008-9-9-r137Version 2.1.1.20160309
Software, algorithmBEDtoolsDOI:10.1002/0471250953.bi1112s47Version 2.27.1
Software, algorithmSAMtoolsDOI:10.1093/bioinformatics/btp352 Version 1.19
Software, algorithmGREATDOI: 10.1038/nbt.1630Version 4.0.4
Software, algorithmHOMERPMID:20513432Version 1.9
Software, algorithmGenomatix
OtherATAC-seq datasetThis paperNCBI GEO (GSE144306)
Otherhistone modification datasets
H3K27ac, H3K4me1, H3K4me3, H3K9me3
ENCODE Consortium and the ENCODE production laboratoryENCSR905FFU
ENCSR426EZM
ENCSR416OYH
ENCSR022DED
ChIP-seq
OtherphastConsElements60wayPlacentalUCSC site; Casper et al., 2018

Data availability

The data was submitted to GEO (GSE144306). Additional data of scRNA-seq was submitted to GEO (GSE160090).

The following data sets were generated
The following previously published data sets were used
    1. Ren B
    2. UCSD
    (2015) ENCODE
    ID ENCSR905FFU. H3K27ac ChIP-seq on embryonic 13.5 day mouse limb.
    1. Ren B
    2. UCSD
    (2015) ENCODE
    ID ENCSR426EZM. H3K4me1 ChIP-seq on embryonic 13.5 day mouse limb.
    1. Ren B
    2. UCSD
    (2015) ENCODE
    ID ENCSR416OYH. H3K4me3 ChIP-seq on embryonic 13.5 day mouse limb.
    1. Ren B
    2. UCSD
    (2015) ENCODE
    ID ENCSR022DED. H3K9me3 ChIP-seq on embryonic 13.5 day mouse limb.

References

    1. Katz JP
    2. Perreault N
    3. Goldstein BG
    4. Lee CS
    5. Labosky PA
    6. Yang VW
    7. Kaestner KH
    (2002)
    The zinc-finger transcription factor Klf4 is required for terminal differentiation of goblet cells in the Colon
    Development 129:2619–2628.
    1. Lehman C
    2. Cuomo F
    3. Kummer FJ
    4. Zuckerman JD
    (1995)
    The incidence of full thickness rotator cuff tears in a large cadaveric population
    Bulletin (Hospital for Joint Diseases 54:30–31.
  1. Software
    1. R Development Core Team
    (2013) R: A language and environment for statistical computing
    R Foundation for Statistical Computing, Vienna, Austria.
    1. Thomopoulos S
    2. Genin GM
    3. Galatz LM
    (2010)
    The development and morphogenesis of the tendon-to-bone insertion - what development can teach us about healing
    Journal of Musculoskeletal & Neuronal Interactions 10:35–45.

Decision letter

  1. Shahragim Tajbakhsh
    Reviewing Editor; Institut Pasteur, France
  2. Mone Zaidi
    Senior Editor; Icahn School of Medicine at Mount Sinai, United States

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

Your study describes bi-fated cells at the murine tendon-to-bone attachment, namely with activation of a mixture of chondrocyte and tenocyte transcriptomes, under the regulation of shared regulatory elements and KLF transcription factors, notably KLF2 and KLF4. The report provides novel insights into the yet unknown molecular and cellular architecture of the tendon-to-bone attachments, and hence is seen to be both novel and medically relevant.

Decision letter after peer review:

Thank you for submitting your article "Bi-fated tendon-to-bone attachment cells are regulated by shared enhancers and KLF transcription factors" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Clifford Rosen as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Senior Editor has drafted this decision to help you prepare a revised submission.

The manuscript from Kult and colleagues focuses on the transcriptional regulation of a unique set of cells located between bone and tendon known as the attachment unit. Using RNA-seq and ATAC-seq, the authors investigate shared and unique transcriptional programs and accessible genomic regions. The ATACseq and epigenome profiling reveals transcriptional enhancers, with overlapping intergenic areas between bifated and both fates. Transgenic enhancer reporters expose common enhancers for tenocytes and chondrocytes. Klf2 and Klf4 were identified as critically required for differentiation. These findings lead the authors to propose that bi-fated attachment cells that connect tendon to bone. Previous studies showed that these cells express cartilage and tendon markers so the present study would need to clearly highlight the advances made compared to previous work. As a general note for example, the co-option of existing enhancers does not rule out the existence of de novo ones. This needs to be addressed clearly in the study.

Essential revisions:

1) The test of enhancers by transgenesis is somewhat limited in scope (only 8 tested) and yields some surprising results that can be explored further. An important request would be to make sure of the reproducibility of experiments. It is said that 5 out of 8 selected elements drove expression in the forelimb, but not how many attempts were done for the negative ones. Moreover, as reported in the transparent report form, N=1 for Col1a1 element, N=3 for Klf2 element, N=2 for Sox9 element, N=3 for Mgp element, N=4 for Col11a1. The reviewers understand that a site-directed approach is likely to be more reproducible that random insertion, it is recommended to examine at least 3 instances per element, to be certain of some surprising results. For example, although Sox9 is not expressed in hypertrophic chondrocytes, the Sox9 element drives expression in hypertrophic chondrocytes. Moreover, Sox9 and Klf2 elements drive expression in hypertrophic chondrocytes but not in AU cells. If these results are confirmed, they may cast doubt on the conclusion that co-option of enhancers is the (only) mechanism that regulates expression in AU cells.

2) The authors have not systemically looked for AU enhancers that are not shared with tenocytes or chondrocytes. Combined ATAC-seq dataset and published ChIP-seq of histone marks can potentially identify new enhancers. The authors could speculate or assess if those enhancers were acquired de novo and exclusive of AU cells.

3) Figure 1A, the authors present a PCA biplot : Can they be more specific on how the data were transformed prior the dimension reduction (FPKM, VST, Log transformed, CPM…?). How many genes were taken into account in the PCA; All or as it is more commonly done on the 500 most variant ones ?

4) Figure 1 C – GO terms found are very generic, this information does not really seem to be useful. Can the authors can be more specific on the parameters they used in their GSEA analysis : test used and p-value correction (FDR q-value suggests a Benjamin and Hochberg correction, it that right ?

5) In general, one has to give the detail on the software version (including package version) and OS type used for the bioinformatic analysis, these informations are missing from the manuscript.

6) In the sentence: "This suggests that the attachment cell transcriptome is largely shared with both chondrocytes and tenocytes (Figure 1A, PC1 52.47%)", the word largely is misleading.

7) The authors do not discuss the variation both on the first and the second PC and of the attachment samples. This is a big issue because there are only 2 samples for this category of cells in which the intra-group variability is very high. This leads to a poor statistical parameter estimation giving rise to poor statistical test outcome.

8) Legend of Figure 1: The term MARS-Seq is slightly misleading as it is usually associated with single cell RNAseq analysis. For clarity, please write instead bulk-MARS-Seq.

9) "To further support our initial observation that the transcriptome of the attachment cells is a mixture of chondrocyte and tenocyte transcriptomes, we clustered the statistically significant differentially expressed genes between all samples into 5 clusters, using CLICK". Same remark as the use of contrasts in DESeq2.

10) In the sentence : "From these two clusters, 374 genes, 320 of them tenogenic and 54 chondrogenic, were also found to be expressed by attachment cells." It is unclear what "so found to be expressed by attachment cells" mean? For instance, for the tenogenic markers, does this mean that in the attachment vs. chondrocytes comparison these genes are up-regulated in the attachment cells? In that case, how are the tenogenic markers defined, using the tenocytes vs chondrocytes comparison? Would it be possible to have a Venn diagram to help follow the process to define the different marker identifications? Has a simpler method using contrasts in DESeq2 been tested? If yes, do the results converge with the ones presented here? Are these "statistically significant differentially expressed genes between all samples" coming from a pairwise wald-test or a likelihood ratio test? the Materials and methods suggest that the wald-test was used. Please clarify.

11) In a previous study the authors investigated the emergence of the attachment unit (AU) with focus on bone eminence progenitors (co-expressing Sox9 and SCX up to E12.5 and expressed Col2 after E12.5 according to Col2a1CreERT2 lineage tracing). Here, they focus on the transcriptome of E14.5 attachment cells from the deltoid tuberosity, however these appear to be different from tuberosity progenitors (adjacent chondrocytes) as described in Figure 1—figure supplement 1. A better definition of what is defined as attachment unit in this paper vs previous papers and/or the AU subcompartments would help clarify the populations that are being examined.

12) Moreover, as the constitutive Col2a1Cre did not label the AU, but in the previous study did label the AU/bone eminence progenitors, it is unclear what the exact definition of AU is.

13) When using the Col2a1-Cre, R26R-tdTomato and Scx-GFP, the authors mention : Unexpectedly, the two reporters failed to label the attachment cells that were located in between these two populations. This failure might be due to a missing regulatory element in one of the constructs that was used to produce each transgenic reporter. However, in Figure 1—figure supplement 1, subpopulation 5 seems to have SCX+ cells. Is this an error of labelling? What is the orientation of this section?

14) For FACS and ATACseq analysis the authors use Sox9CreERT2;tdTomato;SCX-GFP and Col2CreERT2;tdTomato;SCX-GFP. It is not clear why for FACS the Col2CREERT2 line is used while for LCM the constitutive one is used. Moreover, as previously reported, they isolate attachment cells as double positive SOX9/SCX cells. Here again, do the cells taken for analysis include those of the tuberosity itself? Col2CreERT2 with Tamox at E12.5 should be labeling the tuberosity too. It is unclear which cells from which Cre/reporter combination have been used for the ATACseq experiments of Figure 3.

15) For Figure 2, a scheme showing where exactly in the bone we are located and how it has been sectioned would be helpful. Also, it would be nice to perform single molecule FISH on top of Col2Cre:R26TOM:SCX lineage tracing to show the specificity of the colocalization in the "double reporter-negative" area. Also, including the KLF2/4 FISH at this point would help visualize distinctions between genes belonging to cluster 5 (unique to AU) vs genes referred as mixed transcriptome (Wwp2, Bgn).

16) The authors propose a role of Klf2/4 in attachment differentiation. What is the temporality in expression of Klf2/4 vs the putative downstream factors such as Gli1, Col5a1?

17) In ISH figures, some cells in the cartilage compartment also seem to coexpress tenocyte/cartilage markers. Can the authors comment on that?

18) How did the authors adapt MARS-Seq (a single cell RNA seq pipeline taking advantage of cell sorting) to a bulk analysis? More specifically, it isn't clear how laser capture technique was combined with the MARS-seq protocol.

19) The resolution on the single molecule FISH does not allow to really appreciate a large coexpression of the presented markers in the area.

20) Figure 4: Could the authors indicate more clearly the demarcation between cartilage and connective tissue where double labeling is found?

21) It would be interesting to know if the loss of the gene expression in Figure 5 results in a morphological abnormal attachment at later postnatal stages. If the authors have looked, it would be helpful to comment in the Discussion or include the data. How much of the intermediate gene expression program in the attachment site is dependent upon Klf regulation?

22) Is this transcriptional state-sharing permanent or transitional? Their work could be nicely contrasted and compared with some studies examining transcriptional heterogeneity/the co-expression of multiple cell fates as a mechanism cells used to transition from (multipotent) progenitor states to committed fates. Enthesis tissue would be an interesting and unique situation where possibly this intermediate shared transcriptional state is maintained to generate a new cell type. Possible references for transcriptional heterogeneity in progenitors include: Soldatov et al., 2019 and Johnson et al., 2015.

23) The section on the AEG/esophagus-stomach boundary should be better integrated with their own data or removed from the Discussion. It was not clearly stated how these two tissues are similar other than being border tissues. It is recommended to expand this section to include more specific examples how these regions (enthesis and esophagus) are related. Perhaps this esophageal boundary has also been shown to have a shared transcriptional/epigenetic state with neighboring tissues?

https://doi.org/10.7554/eLife.55361.sa1

Author response

Essential revisions

1) The test of enhancers by transgenesis is somewhat limited in scope (only 8 tested) and yields some surprising results that can be explored further. An important request would be to make sure of the reproducibility of experiments. It is said that 5 out of 8 selected elements drove expression in the forelimb, but not how many attempts were done for the negative ones. Moreover, as reported in the transparent report form, N=1 for Col1a1 element, N=3 for Klf2 element, N=2 for Sox9 element, N=3 for Mgp element, N=4 for Col11a1. The reviewers understand that a site-directed approach is likely to be more reproducible that random insertion, it is recommended to examine at least 3 instances per element, to be certain of some surprising results. For example, although Sox9 is not expressed in hypertrophic chondrocytes, the Sox9 element drives expression in hypertrophic chondrocytes. Moreover, Sox9 and Klf2 elements drive expression in hypertrophic chondrocytes but not in AU cells. If these results are confirmed, they may cast doubt on the conclusion that co-option of enhancers is the (only) mechanism that regulates expression in AU cells.

We appreciate the reviewer’s comments about the results of the site-directed LacZ enhancer-reporter assay. The paper that describes the new site-directed transgenic analysis method that we used (termed “enSERT”, for “enhancer INSERTION”) has been published very recently (Kvon et al., 2020). This publication informs about the superior specificity and reproducibility of this new approach, as verified by extensive analyses in mouse embryos. This method, which uses CRISPR/Cas9 for site-directed integration (knock-in) of LacZ reporter transgenes at the H11 safe-harbor locus, allows in vivo validation of enhancer activity with higher accuracy and specificity than previous, mostly random integration-based methods. Specifically, enSERT allows to achieve “an average transgenic rate of 50% for transgenes as large as 11 kb (based on >1,200 transgenic mice resulting from the injection of >150 independent transgenic constructs), compared to a 12% transgenic rate observed in conventional random transgenesis (based on >3,000 independent transgenic constructs)” (Kvon et al., 2020).

In Kvon et al. (Figure S2), the reliable reproducibility of the method was confirmed at E11.5 in a wide range of embryonic tissues, e.g. limbs, brain and heart. Typically, site-specific multi-copy transgenic insertion events resulted in the most sensitive (and reproducible) patterns of enhancer activity (see Materials and methods subsection “Tandem integration at H11” in Kvon et al.).

In our study, we leveraged the reliable nature of enSERT to validate in vivo multiple enhancer regions at E14.5 that had been predicted by our ATAC-seq experiments. The majority of tested regions indeed revealed in vivo enhancer activities with superior reproducibility in multiple tissues. For example, while mm1988 (Klf2 element) displayed reproducible patterns in craniofacial tissues and limbs, mm1990 (Mgp element) had reproducible patterns specifically in limbs (n=3/3 for both).

Regarding the number of examined embryos, we updated and increased the number of repetitions. We examined at least 3 embryos per enhancer and, in some cases, between 6 and 11 repetitions were made. As the reviewers suggested, we have also added the transgenic reproducibility numbers (X/Y) to Figure 4. As is the standard in the field, we compare the number of embryos with sub-regionally reproducible LacZ staining (X) to the total number of transgenic embryos obtained per tested genomic element.

As shown in Author response table 1, while some of the elements, such as the Klf2 element mm1988, displayed a high and extremely convincing reproducibly rate (N=3/3), limb-specific reproducible staining for the Col1 element mm1995 was observed in only 2 out of 9 embryos. Nevertheless, because we confirmed the limb-subregional reproducibility (on sections) in two independent biological replicates and because this element was selected based on our limb-specific ATAC-seq and ENCODE datasets of histone modification marks associated with enhancers and promoters (H3K27Ac and H3K4me1 ChIP-seq), we have reasons to believe that this result is real and represents the activity of a Col1 enhancer in attachment and tendon cells.

Author response table 1
Vista IDNReproducible stainingLimbNotes:
mm19863/6CraniofacialNegative
mm19952/9LimbPositive** Scored as “random integration” (Kvon et al., Cell 2020).
mm19883/3VariousPositive
mm19894/7LimbPositive
mm19903/3LimbPositive
mm199111/11Limb, other structures with less penetrancePositive
mm19920/6NoneNegative2 embryos displayed non-reproducible or background LacZ staining
mm19940/9NoneNegative5 embryos displayed non-reproducible or background LacZ staining

One explanation for the reduced reproducibility of Col1 element may be the temporal in vivo changes in enhancer usage during development. Previous studies have shown that enhancers exhibited tightly temporally restricted predicted activity windows (Nord et al., Cell 2013). It is therefore possible that mm1995 element activity at E14.5 exemplifies the start or end of this element activity window. We therefore agree that in order to establish mm1995 as an enhancer, additional repetitions are required. Unfortunately, due to the COVID-19 situation, we do not know when this will be possible. Therefore, if the reviewers approve, we would like to keep this result in our manuscript. Obviously, we clearly state its limitation.

We fully agree with the reviewers that de novo enhancers may play a role in the regulation of attachment cells in addition to co-option of enhancers (i.e. sharing of elements with tendon or cartilage cells). This issue is discussed in detail in our next answer.

2) The authors have not systemically looked for AU enhancers that are not shared with tenocytes or chondrocytes. Combined ATAC-seq dataset and published ChIP-seq of histone marks can potentially identify new enhancers. The authors could speculate or assess if those enhancers were acquired de novo and exclusive of AU cells.

We thank the reviewer for this comment. It is important to emphasize that we do not exclude the possibility that attachment cells express a unique set of genes. However, our main purpose in this manuscript is to describe the dual fate of attachment cells. Yet, to accommodate the reviewer’s concern, we have included a list of genes that we identified to be expressed specifically in the AU. In addition, we also identified a group of intergenic elements that were accessible specifically in attachment cells. These elements may act as enhancers that drive gene expression in the attachment cells, thereby providing another level of specificity to the regulation of the development of this unique tissue. We have added a table listing these intergenic regions (newly added Supplementary file 1). Among the genes that are associated with AU-specific elements, we identified several transcription factors such as Klf4, Runx3 and Nfia, in addition to ECM-associated genes such as Col11a1 and Col1a1. As the reviewers suggested, we have added text to the Results and Discussion sections describing these findings.

3) Figure 1A, the authors present a PCA biplot : Can they be more specific on how the data were transformed prior the dimension reduction (FPKM, VST, Log transformed, CPM…?). How many genes were taken into account in the PCA; All or as it is more commonly done on the 500 most variant ones ?

The PCA shows log transformed normalized data of the 500 most variant genes. The analysis was done using a modified plotPCA function of DESeq2 (rlogTransformation with parameter blind=TRUE).

4) Figure 1 C – GO terms found are very generic, this information does not really seem to be useful. Can the authors can be more specific on the parameters they used in their GSEA analysis : test used and p-value correction (FDR q-value suggests a Benjamin and Hochberg correction, it that right ?

As the reviewer suggested, we excluded the GO analysis from the revised version of the paper. We used GOrilla to find GO annotations by loading the list of 374 genes on the background of all genes and received the GO processes that were shown in Figure 1C (please find the report with all results attached as GOPROCESS_374.txt).

5) In general, one has to give the detail on the software version (including package version) and OS type used for the bioinformatic analysis, these informations are missing from the manuscript.

As requested by the reviewer, the following information has been added to the Materials and methods section:

“The data were analyzed using the Pipeline Pilot-designed pipeline for transSeq (by INCPM, https://incpmpm.atlassian.net/wiki/spaces/PUB/pages/36405284/tranSeq+on+Pipeline-Pilot). Briefly, the analysis included adapter trimming, mapping to the mm9 genome, collapsing of reads with the same unique molecular identifiers (UMI) of 4 bases (R2) and counting of the number of reads per gene with HTseq-count (Anders et al., 2015), using the most 3’ 1000 bp of each RefSeq’s transcript. DESeq2 (version 1.4.5, Love, Huber and Anders, 2014) was used for normalization and differential expression analysis with betaPrior set to true, cooksCutoff=FALSE, independentFiltering=FALSE. Benjamini-Hochberg method was used to adjust the raw p-values for multiple testing. Genes with adjusted p-value ≤ 0.05 and fold change ≥ 2 between every two conditions were considered as differential. PCA analysis was done using log transformed normalized data (DESeq2 function rlogTransformation with parameter blind=TRUE) was done with a modified plotPCA function of DESeq2.”

6) In the sentence: "This suggests that the attachment cell transcriptome is largely shared with both chondrocytes and tenocytes (Figure 1A, PC1 52.47%)", the word largely is misleading.

The word “largely” was removed from the PCA description.

7) The authors do not discuss the variation both on the first and the second PC and of the attachment samples. This is a big issue because there are only 2 samples for this category of cells in which the intra-group variability is very high. This leads to a poor statistical parameter estimation giving rise to poor statistical test outcome.

We thank the reviewer for this comment. Indeed, we fully agree that more samples would have increased our statistical power and possibly identify more AU-specific genes. However, obtaining this tissue was very challenging technically. Nevertheless, applying adjusted p-values to these samples provided us with a list of statistically significant genes that are differentially expressed in attachment cells. Importantly, the PCA results were further validated by an extensive set of experiments, such as in situ hybridization and scRNA-Seq.

8) Legend of Figure 1: The term MARS-Seq is slightly misleading as it is usually associated with single cell RNAseq analysis. For clarity, please write instead bulk-MARS-Seq.

We have changed the term MARS-Seq to bulk-MARS-Seq in the manuscript.

9) "To further support our initial observation that the transcriptome of the attachment cells is a mixture of chondrocyte and tenocyte transcriptomes, we clustered the statistically significant differentially expressed genes between all samples into 5 clusters, using CLICK". Same remark as the use of contrasts in DESeq2.

As requested by the reviewer, this text have been added to the Materials and methods section: “Clustering of the log normalized read count of differentially expressed genes was done using click algorithm (Expander package version 7.1, Ulitsky et al., 2010), followed by visualization by R (R Core Team, 2013).”

10) In the sentence : "From these two clusters, 374 genes, 320 of them tenogenic and 54 chondrogenic, were also found to be expressed by attachment cells." It is unclear what "so found to be expressed by attachment cells" mean? For instance, for the tenogenic markers, does this mean that in the attachment vs. chondrocytes comparison these genes are up-regulated in the attachment cells? In that case, how are the tenogenic markers defined, using the tenocytes vs chondrocytes comparison? Would it be possible to have a Venn diagram to help follow the process to define the different marker identifications? Has a simpler method using contrasts in DESeq2 been tested? If yes, do the results converge with the ones presented here? Are these "statistically significant differentially expressed genes between all samples" coming from a pairwise wald-test or a likelihood ratio test? the Materials and methods suggest that the wald-test was used. Please clarify.

Genes were selected as follows: First, DESeq2 yielded 865 differential genes (i.e., between all samples, Figure 1—figure supplement 2), based on criteria that are described in Materials and methods (“Genes with adjusted p-value ≤ 0.05 and fold change ≥ 2 between every two conditions were considered as differential”). Second, to demonstrate how attachment cells function as an intermediate tissue, from the 865 differential genes we selected the genes in cluster 1 (up-regulated in tenocytes vs. chondrocytes samples, hence “tenogenic”) and cluster 2 (up-regulated in chondrocytes vs. tenocytes samples, hence “chondrogenic”). Third, we averaged the normalized number of reads in attachment samples of each gene on this list and selected genes with 30 normalized reads or more. Thus, the 374 genes that are shown in Figure 1 represent the activation of both tenogenic and chondrogenic gene expression in E14.5 attachment cells. We have added this description to the Materials and methods section.

11) In a previous study the authors investigated the emergence of the attachment unit (AU) with focus on bone eminence progenitors (co-expressing Sox9 and SCX up to E12.5 and expressed Col2 after E12.5 according to Col2a1CreERT2 lineage tracing). Here, they focus on the transcriptome of E14.5 attachment cells from the deltoid tuberosity, however these appear to be different from tuberosity progenitors (adjacent chondrocytes) as described in Figure 1—figure supplement 1. A better definition of what is defined as attachment unit in this paper vs previous papers and/or the AU subcompartments would help clarify the populations that are being examined.

Indeed, the reviewer accurately cites our previous paper and we understand how the use of AU in both papers can lead to some confusion regarding its definition. We have added to the Introduction a paragraph that explains the development of the AU, serving as the enthesis primordium, whose cells co-express Sox9 and Scx. As development proceeds, this field is divided into sub compartments, namely the cartilaginous bone eminence (Col2-positive) on one side, the tendon (Scx-positive) on the other side and, connecting between them, the attachment cells, which will eventually give rise to the enthesis. The objective of this paper is to address the understudied differentiation sequence of this third population.

12) Moreover, as the constitutive Col2a1Cre did not label the AU, but in the previous study did label the AU/bone eminence progenitors, it is unclear what the exact definition of AU is.

In the previous study, as in this work, the constitutive Col2a1Cre labeled the bone side of the AU, namely the bone eminence. In this work, the term AU refers to the tissue that forms between tendon and cartilage, as explained in the previous response.

13) When using the Col2a1-Cre, R26R-tdTomato and Scx-GFP, the authors mention : Unexpectedly, the two reporters failed to label the attachment cells that were located in between these two populations. This failure might be due to a missing regulatory element in one of the constructs that was used to produce each transgenic reporter. However, in Figure 1—figure supplement 1, subpopulation 5 seems to have SCX+ cells. Is this an error of labelling? What is the orientation of this section?

The orientation of this section is transverse. When observing a series of transverse sections from E14.5 forelimbs, we identified in each section Col2+ labeled cells, Scx-GFP+ labeled cells and, between them, cells that were unlabeled. Figure 1—figure supplement 1 is one of many collected sections in which subpopulation 5 is of unlabeled cells. Regarding GFP-labelled cells at the margins of the laser-captured area, to maximize the size of extracted tissue, we stretched the borders of the AU all the way to the flanking Scx-GFP+ cells. When the laser cuts the piece of tissue it burns the boundaries, so the RNA in cells located at the boundaries is degraded.

14) For FACS and ATACseq analysis the authors use Sox9CreERT2;tdTomato;SCX-GFP and Col2CreERT2;tdTomato;SCX-GFP. It is not clear why for FACS the Col2CREERT2 line is used while for LCM the constitutive one is used. Moreover, as previously reported, they isolate attachment cells as double positive SOX9/SCX cells. Here again, do the cells taken for analysis include those of the tuberosity itself? Col2CreERT2 with Tamox at E12.5 should be labeling the tuberosity too. It is unclear which cells from which Cre/reporter combination have been used for the ATACseq experiments of Figure 3.

When we initiated the LCM experiment, we assumed that Col2 would label the attachment cells. However, to our surprise, attachment cells were Col2-negative. Nevertheless, as the Scx and Col2 reporters clearly defined the boundaries between cell types, we could still use this line for LCM. This part is explained in detail in the Materials and methods section.

Because for the ATAC-seq experiment we needed to isolate the cells using FACS, we had to change our labeling strategy. Therefore, we used a combination of Sox9-CreER and Scx-GFP mice. This indeed worked well for tendon and attachment cells, but the number of isolated chondrocytes was insufficient. Therefore, we used the Col2-CreERT line to label chondrocytes and isolate them. To clarify this issue, we have added text to the Results section.

As the reviewer noted, this labeling should also label the chondrocytes that compose the tuberosity. However, due to the large difference in cell numbers between the tuberosity and the rest of the cartilage that was collected, we think that this had a limited (if any) effect on our results.

15) For Figure 2, a scheme showing where exactly in the bone we are located and how it has been sectioned would be helpful. Also, it would be nice to perform single molecule FISH on top of Col2Cre:R26TOM:SCX lineage tracing to show the specificity of the colocalization in the "double reporter-negative" area. Also, including the KLF2/4 FISH at this point would help visualize distinctions between genes belonging to cluster 5 (unique to AU) vs genes referred as mixed transcriptome (Wwp2, Bgn).

We have added to Figure 1 a scheme showing how we sectioned the deltoid and great tuberosities for FISH. Using E14.5 ScxGFP-Col2Cre-tdTomato transverse sections for Bgn and Wwp2 double smFISH (conjugated to Quasar 670 and CAL Fluor Red 610, respectively) is problematic, since tdTomato (540nm-581nm) is a very strong reporter that may mask the signal of the Wwp2 probe coupled to Quazar 610 (590nm-610nm), due to a "bleed through" of the tdTomato into the Quazar610 filter.

As an alternative, to address the reviewer’s comment we performed a 10x genomics experiment on Sox9+Scx+ cells, where we could search for the transcripts of Bgn and Wwp2 in addition to expression of Klf2/Klf4 in Sox9+Scx+ cells. Figure 2A,B summarizes this analysis. In short, this experiment exemplifies that Sox9+Scx+ cells indeed express both chondrogenic and tenogenic markers, in addition to KLF’s (Figure 5B).

16) The authors propose a role of Klf2/4 in attachment differentiation. What is the temporality in expression of Klf2/4 vs the putative downstream factors such as Gli1, Col5a1?

To address the temporality in gene expression, we performed transcriptome analysis (bulk MARS-Seq) of FACS-sorted attachment cells E13.5, when it is first possible anatomically to identify the attachment site (Sox9+Scx+ cells, n=3 for each gene; see Author response table 2). Results showed that at that time point, both Klf2 and Klf4, as well as Col5a1 and Gli1, are already expressed. Since our analysis of the Klf2/4 double mutants clearly showed a reduction (Figure 5) in expression of Gli1 and Col5a1, bona fide markers of the differentiating enthesis (Felsenthal et al., 2018; Schwartz et al., 2015), we concluded that Klf2 and Klf4 are necessary for attachment cell differentiation.

Author response table 2
GeneBulk MARS-Seq normalized reads (mean)
Klf2203.67
Klf4384.75
Gli1349.15
Col5a1391.62

17) In ISH figures, some cells in the cartilage compartment also seem to coexpress tenocyte/cartilage markers. Can the authors comment on that?

We agree with the reviewer that in some cases, at the borders the expression patterns of tendon and cartilage markers are not sharply segregated. This is mainly because during development, the borders between the different tissues that compose the attachment site are not sharp and clear, as the process of compartmentalization is still underway. Nevertheless, we believe that in general, expression of ISH markers followed the expected spatial patterns. The scRNA-seq results further support this conclusion.

18) How did the authors adapt MARS-Seq (a single cell RNA seq pipeline taking advantage of cell sorting) to a bulk analysis? More specifically, it isn't clear how laser capture technique was combined with the MARS-seq protocol.

MARS-Seq was indeed originally developed as a FACS-based method for single-cell RNA sequencing [69] (Jaitin et al., 2014; Keren-Shaul, 2019). Due to the low RNA input in scRNA-seq, the authors further developed a bulk MARS-seq protocol, which is an adaptation of the single-cell MARS-seq. The input for bulk MARS-seq is clean RNA and hence it can be used for expression profiling by purifying RNA from LCM-isolated samples, FACS-sorted cells and more.

As explained in detail in the Materials and methods subsection “Laser capture microdissection”, following LCM of E14.5 cryo-sections of the attachment site, we purified the RNA using an RNA purification kit (RNeasy FFPE Kit, Qiagen). The resulting RNA was the input for the bulk MARS-seq protocol. Briefly, RNA from each sample was barcoded during reverse transcription and pooled. Following Agencourct Ampure XP beads cleanup (Beckman Coulter), the pooled samples underwent second strand synthesis and were linearly amplified by T7 in vitro transcription. The resulting RNA was fragmented and converted into a sequencing-ready library by tagging the samples with Illumina sequences during ligation, RT, and PCR. Libraries were quantified by Qubit and TapeStation as well as by qPCR for actb gene as previously described (Jaitin et al., 2014; Keren-Shaul et al., 2019). Sequencing was done on a Hiseq 2500 SR50 cycles kit (Illumina).

19) The resolution on the single molecule FISH does not allow to really appreciate a large coexpression of the presented markers in the area.

Figure 2C, showing the smFISH of biglycan and Wwp2 in E14.5 attachment cells, was increased in resolution to easily appreciate the coexpression of the markers.

20) Figure 4: Could the authors indicate more clearly the demarcation between cartilage and connective tissue where double labeling is found?

We have added clearer demarcation between cartilage and connective tissue.

21) It would be interesting to know if the loss of the gene expression in Figure 5 results in a morphological abnormal attachment at later postnatal stages. If the authors have looked, it would be helpful to comment in the Discussion or include the data. How much of the intermediate gene expression program in the attachment site is dependent upon Klf regulation?

We agree with the reviewers. However, due to postnatal lethality, our analysis of E18.5 is the latest time point that can be studied in this mouse line.

22) Is this transcriptional state-sharing permanent or transitional? Their work could be nicely contrasted and compared with some studies examining transcriptional heterogeneity/the co-expression of multiple cell fates as a mechanism cells used to transition from (multipotent) progenitor states to committed fates. Enthesis tissue would be an interesting and unique situation where possibly this intermediate shared transcriptional state is maintained to generate a new cell type. Possible references for transcriptional heterogeneity in progenitors include: Soldatov et al., 2019 and Johnson et al., 2015.

As suggested, we have added a paragraph to the Discussion regarding this issue.

23) The section on the AEG/esophagus-stomach boundary should be better integrated with their own data or removed from the Discussion. It was not clearly stated how these two tissues are similar other than being border tissues. It is recommended to expand this section to include more specific examples how these regions (enthesis and esophagus) are related. Perhaps this esophageal boundary has also been shown to have a shared transcriptional/epigenetic state with neighboring tissues?

This part of the Discussion was edited to fit with the reviewer requirements.

https://doi.org/10.7554/eLife.55361.sa2

Article and author information

Author details

  1. Shiri Kult

    Department of Molecular Genetics, Weizmann Institute of Science, Rehovot, Israel
    Contribution
    Conceptualization, Validation, Investigation, Visualization, Methodology, Writing - original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8216-2908
  2. Tsviya Olender

    Department of Molecular Genetics, Weizmann Institute of Science, Rehovot, Israel
    Contribution
    Data curation, Software, Formal analysis, Validation, Visualization, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  3. Marco Osterwalder

    1. Environmental Genomics and Systems Biology Division, Lawrence Berkeley National, Berkeley, United States
    2. Department for BioMedical Research (DBMR), University of Bern, Bern, Switzerland
    Contribution
    Validation, Investigation, Visualization, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1969-2313
  4. Svetalana Markman

    Department of Molecular Genetics, Weizmann Institute of Science, Rehovot, Israel
    Contribution
    Visualization, Writing - review and editing
    Competing interests
    No competing interests declared
  5. Dena Leshkowitz

    Life Sciences Core Facilities, Weizmann Institute of Science, Rehovot, Israel
    Contribution
    Data curation, Software, Formal analysis, Validation, Visualization, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  6. Sharon Krief

    Department of Molecular Genetics, Weizmann Institute of Science, Rehovot, Israel
    Contribution
    Validation, Investigation, Visualization, Methodology
    Competing interests
    No competing interests declared
  7. Ronnie Blecher-Gonen

    Department of Immunology, Weizmann Institute of Science, Rehovot, Israel
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
  8. Shani Ben-Moshe

    Department of Molecular Cell Biology, Weizmann Institute of Science, Rehovot, Israel
    Contribution
    Investigation, Visualization, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  9. Lydia Farack

    Department of Molecular Cell Biology, Weizmann Institute of Science, Rehovot, Israel
    Contribution
    Software, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9597-5078
  10. Hadas Keren-Shaul

    Life Sciences Core Facilities, Weizmann Institute of Science, Rehovot, Israel
    Contribution
    Investigation, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  11. Tomer-Meir Salame

    Life Sciences Core Facilities, Weizmann Institute of Science, Rehovot, Israel
    Contribution
    Data curation, Formal analysis, Validation, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  12. Terence D Capellini

    Department of Human Evolutionary Biology, Harvard University, Department of Human Evolutionary Biology, United States; Broad Institute of Harvard and MIT, Cambridge, United States
    Contribution
    Resources, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3842-8478
  13. Shalev Itzkovitz

    Department of Molecular Cell Biology, Weizmann Institute of Science, Rehovot, Israel
    Contribution
    Resources, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0685-2522
  14. Ido Amit

    Department of Immunology, Weizmann Institute of Science, Rehovot, Israel
    Contribution
    Resources, Methodology
    Competing interests
    No competing interests declared
  15. Axel Visel

    1. Environmental Genomics and Systems Biology Division, Lawrence Berkeley National, Berkeley, United States
    2. U.S. Department of Energy Joint Genome Institute, Lawrence Berkeley National Laboratory, Berkeley, United States
    3. School of Natural Sciences, University of California, Merced, Merced, United States
    Contribution
    Resources, Funding acquisition, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-4130-7784
  16. Elazar Zelzer

    Department of Molecular Genetics, Weizmann Institute of Science, Rehovot, Israel
    Contribution
    Conceptualization, Supervision, Funding acquisition, Methodology, Writing - original draft, Project administration
    For correspondence
    eli.zelzer@weizmann.ac.il
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1584-6602

Funding

National Institutes of Health (R01 AR055580)

  • Elazar Zelzer

National Institutes of Health (R01HG003988)

  • Axel Visel

Israel Science Foundation (grant No. 345/16)

  • Elazar Zelzer

Minerva Foundation (grant No. 713533)

  • Elazar Zelzer

David and Fela Shapell Family Center for Genetic Disorders

  • Elazar Zelzer

Estate of Mr. and Mrs. van Adelsberge

  • Elazar Zelzer

University of California (DEAC02-05CH11231)

  • Marco Osterwalder

Swiss National Science Foundation (PCEFP3_186993)

  • Marco Osterwalder

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Nitzan Konstantin for expert editorial assistance, Dr. Douglas Lutz and service engineer Tal Alon (Getter Bio-Med, Zeiss) for LCM calibration, Neria Sharabi from the Department of Veterinary Resources, Weizmann Institute, and all Zelzer lab members for suggestions and advice. We thank Drs. Merav Kedmi, David Pilzer and Hadas Keren-Shaul from the Genomics Sandbox unit at the Life Science Core Facility, Weizmann Institute of Science, for critical advice. We thank E Sebzda for providing floxed Klf2 mice and the Mutant Mouse Regional Resource Center (MMRRC) at UC Davis for providing floxed Klf4 mice. We thank the ENCODE Consortium and the ENCODE production laboratory for generating the described datasets. This study was supported by grants from the Israel Science Foundation (ISF, grant No. 345/16), National Institute of Health (grant No. R01 AR055580), Minerva Foundation (grant No. 713533), the Estate of Mr. and Mrs. van Adelsbergen, and the David and Fela Shapell Family Center for Genetic Disorders (to EZ). Work at Lawrence Berkeley National Lab was supported by National Institute of Health grant R01HG003988 (to AV) and performed under Department of Energy Contract DE-AC02-05CH11231, University of California. MO was supported by Swiss National Science Foundation grant PCEFP3_186993.

Ethics

Animal experimentation: All mice were maintained and used in accordance with protocols approved by the Weizmann Institutional Animal Care and Use Committee (IACUC; permission no. 10550119-2). All animal work was reviewed and approved by the Lawrence Berkeley National Laboratory (LBNL) Animal Welfare Committee. All mice used in this study were housed at the Animal Care Facility (ACF) at LBNL. Mice were monitored daily for food and water intake, and animals were inspected weekly by the Chair of the Animal Welfare and Research Committee and the head of the animal facility in consultation with the veterinary staff. The LBNL ACF is accredited by the American Association for the Accreditation of Laboratory Animal Care (AAALAC). Transgenic mouse assays were performed in Mus musculus FVB background mice.

Senior Editor

  1. Mone Zaidi, Icahn School of Medicine at Mount Sinai, United States

Reviewing Editor

  1. Shahragim Tajbakhsh, Institut Pasteur, France

Publication history

  1. Received: January 21, 2020
  2. Accepted: November 30, 2020
  3. Version of Record published: January 15, 2021 (version 1)

Copyright

© 2021, Kult et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,388
    Page views
  • 166
    Downloads
  • 4
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Cell Biology
    2. Developmental Biology
    Sun-Hee Hwang et al.
    Research Article

    The role of compartmentalized signaling in primary cilia during tissue morphogenesis is not well understood. The cilia-localized G-protein-coupled receptor—Gpr161 represses hedgehog pathway via cAMP signaling. We engineered a knock-in at Gpr161 locus in mice to generate a variant (Gpr161mut1), which was ciliary localization defective but cAMP signaling competent. Tissue phenotypes from hedgehog signaling depend on downstream bifunctional Gli transcriptional factors functioning as activators/repressors. Compared to knockout (ko), Gpr161mut1/ko had delayed embryonic lethality, moderately increased hedgehog targets and partially down-regulated Gli3-repressor. Unlike ko, the Gpr161mut1/ko neural tube did not show Gli2-activator-dependent expansion of ventral-most progenitors. Instead, the intermediate neural tube showed progenitor expansion that depends on loss of Gli3-repressor. Increased extraciliary receptor (Gpr161mut1/mut1) prevented ventralization. Morphogenesis in limb buds and midface requires Gli-repressor; these tissues in Gpr161mut1/mut1 manifested hedgehog hyperactivation phenotypes—polydactyly and midfacial widening. Thus, ciliary and extraciliary Gpr161 pools likely establish tissue-specific Gli-repressor thresholds in determining morpho-phenotypic outcomes.

    1. Cell Biology
    2. Developmental Biology
    Evelien Eenjes et al.
    Research Article Updated

    SOX2 expression levels are crucial for the balance between maintenance and differentiation of airway progenitor cells during development and regeneration. Here, we describe patterning of the mouse proximal airway epithelium by SOX21, which coincides with high levels of SOX2 during development. Airway progenitor cells in this SOX2+/SOX21+ zone show differentiation to basal cells, specifying cells for the extrapulmonary airways. Loss of SOX21 showed an increased differentiation of SOX2+ progenitor cells to basal and ciliated cells during mouse lung development. We propose a mechanism where SOX21 inhibits differentiation of airway progenitors by antagonizing SOX2-induced expression of specific genes involved in airway differentiation. Additionally, in the adult tracheal epithelium, SOX21 inhibits basal to ciliated cell differentiation. This suppressing function of SOX21 on differentiation contrasts SOX2, which mainly drives differentiation of epithelial cells during development and regeneration after injury. Furthermore, using human fetal lung organoids and adult bronchial epithelial cells, we show that SOX2+/SOX21+ regionalization is conserved. Lastly, we show that the interplay between SOX2 and SOX21 is context and concentration dependent leading to regulation of differentiation of the airway epithelium.