- Research article
- Open Access
Heterogeneous combinatorial expression of Hoxd genes in single cells during limb development
© Duboule et al. 2018
- Received: 25 May 2018
- Accepted: 29 August 2018
- Published: 18 September 2018
Global analyses of gene expression during development reveal specific transcription patterns associated with the emergence of various cell types, tissues, and organs. These heterogeneous patterns are instrumental to ensure the proper formation of the different parts of our body, as shown by the phenotypic effects generated by functional genetic approaches. However, variations at the cellular level can be observed within each structure or organ. In the developing mammalian limbs, expression of Hox genes from the HoxD cluster is differentially controlled in space and time, in cells that will pattern the digits and the forearms. While the Hoxd genes broadly share a common regulatory landscape and large-scale analyses have suggested a homogenous Hox gene transcriptional program, it has not previously been clear whether Hoxd genes are expressed together at the same levels in the same cells.
We report a high degree of heterogeneity in the expression of the Hoxd11 and Hoxd13 genes. We analyzed single-limb bud cell transcriptomes and show that Hox genes are expressed in specific combinations that appear to match particular cell types. In cells giving rise to digits, we find that the expression of the five relevant Hoxd genes (Hoxd9 to Hoxd13) is unbalanced, despite their control by known global enhancers. We also report that specific combinatorial expression follows a pseudo-time sequence, which is established based on the transcriptional diversity of limb progenitors.
Our observations reveal the existence of distinct combinations of Hoxd genes at the single-cell level during limb development. In addition, we document that the increasing combinatorial expression of Hoxd genes in this developing structure is associated with specific transcriptional signatures and that these signatures illustrate a temporal progression in the differentiation of these cells.
- Hox genes
- Gene expression
Chromatin interactions between Hox genes and their enhancers in single cells showed variability [4, 7, 8], and super-resolution microscopy confirmed that the HoxD gene cluster can display a variety of structural conformations in future autopod cells . This heterogeneity is difficult to integrate with chromosome conformation datasets produced at this locus, since the latter approach reflects the averaged behaviors of a cellular population. Consequently, a higher variability can be expected in cell-specific Hox gene transcriptions, when compared to the apparently rather homogenous expression profiles previously reported [6, 8].
Moreover, while genetic approaches have revealed the critical function of these genes during limb outgrowth and patterning, the homogeneous or heterogeneous impact of mutations at the cellular level is more difficult to evaluate. The ablation of Hoxd13 alone leads to a morphological effect in digits weaker than when a simultaneous deletion of Hoxd11, Hoxd12, and Hoxd13 is achieved [10–12], suggesting that Hoxd11, Hoxd12, and Hoxd13 functionally cooperate during digit development. However, how this cooperation occurs at the cellular level is unknown.
One potential cause for transcriptional heterogeneity may involve a competition between the various promoters located in-cis and the global enhancers driving transcription in digits . It was indeed recently reported that the HoxD cluster lies between two large topologically associating domains (TADs) [14, 15] each of them containing series of enhancer elements with distinct specificities [3, 7]. The TAD located centromeric to HoxD (C-DOM) contains several enhancers specific for autopod (digit) cells, whereas T-DOM, the TAD located telomeric to HoxD, hosts a series of enhancers specific for the future arm and forearm cells. While genes located at either extremity of the cluster respond to their neighboring TAD, those genes located at a central position in the cluster such as Hoxd9, Hoxd10, or Hoxd11 are targeted successively by enhancers belonging to the two different TADs. Initially, in future forearm cells, they respond to T-DOM regulation, whereas in a subsequent phase, in future digit cells, they respond to C-DOM enhancers , which may lead to an even greater heterogeneity in transcript distribution. In order to try and evaluate the heterogeneity in Hoxd transcript distribution during limb development, we produced single-limb cell transcriptomes of different origins, to see whether the apparently homogenous Hox gene transcriptional program as observed upon large-scale analyses could be observed at the cellular level as well. Here, we report that Hoxd gene transcripts are present in various combinations in different limb cells. We discuss the impact of these results upon our understanding of how Hoxd genes are regulated and how their global functions are achieved in these structures.
Heterogenous distribution of posterior Hoxd gene transcripts in single cells
In order to document the expression pattern of Hoxd13 at the single-cell level, embryonic day (E) 12.5 limb sections were use in RNA-FISH experiments (Fig. 1b). As expected, we observed a high expression specificity in presumptive digit cells in the distal part of the forelimb, with the highest transcript levels in cells located at the boundary between the digital and the interdigital compartments, while lower levels were scored in interdigital mesenchyme. Signal was detected neither within the digital compartment nor in more proximal parts of the limb  (Fig. 1b). However, a high heterogeneity in gene expression was recorded, with stippled signal pattern contrasting with the broader expression domain previously described by whole-mount in situ hybridization (WISH). As a consequence, we asked whether all cells expressing Hoxd13 would also contain Hoxd11 transcripts, knowing that both genes are under the same regulatory control in these distal cells [2, 17]. We micro-dissected autopod tissue to obtain a single-cell suspension and performed double fluorescent RNA labelling. The single-cell preparation was then analyzed by fluorescence-activated cell sorting (FACS) and revealed that only a minority of cells was in fact expressing Hoxd11 and/or Hoxd13 (Fig. 1c). Amongst positive cells, the largest fraction was Hoxd13 positive and negative for Hoxd11 (d13+d11−; 53%), whereas double positive cells (d13+d11+) represented 38% only and 9% of the cells contained Hoxd11 mRNAs alone (d11+) (Fig. 1c).
Because a substantial number of cells did not express any Hoxd genes, we enriched for the positive fraction using a mouse line containing a GFP reporter sequence knocked in Hoxd11. In these mice, GFP was produced in those cells where Hoxd11 had been transcribed (Additional file 1: Fig. S1). We monitored the fluorescence at E12.5 and observed a pattern recapitulating Hoxd11 endogenous expression (Fig. 1d). E12.5 limb cells from these animals were FACS-sorted using the GFP (Fig. 1e) and, under these conditions, the double labelling of GFP-positive cells increased to more than a third of the cells (Fig. 1f). However, amongst the positive cells, the ratio between the three Hoxd-positive populations (Hoxd13 only, Hoxd11 only, and double positive) was roughly the same as before (37%, 7%, and 55%, respectively). To confirm the presence of these different populations, we performed Hoxd13 RNA-FISH on sections from Hoxd11::GFP E12.5 forelimbs (Fig. 1g) and observed a high variability in GFP levels (Fig. 1h). We found that high levels of Hoxd13 were observed in cells with either little or no Hoxd11 activity (Fig. 1i), yet the majority of cells displayed high signals for both Hoxd11 and Hoxd13, suggesting that in these cells the two genes were regulated in a similar manner.
To quantify a potential correlation between Hoxd11 and Hoxd13 expression levels in these GFP-positive cells, we binned Hoxd11-positive cells in three categories: negative cells (d11neg, orange), cells expressing at low levels (d11low, red), and cells expressing at high levels (d11hi, gray; Fig. 1j, left panel). Flow cytometry analysis revealed that higher Hoxd13 levels were clearly observed in the d11hi population, indicating that in single cells, whenever both genes are expressed, they tend to respond to enhancers with the same efficiency (Fig. 1j). To relate these latter results with the level of GFP observed by microscopy (Fig. 1i), we monitored the levels of GFP in single cells and found a correlation between abundant Hoxd11 mRNAs, on the one hand, and higher levels of the GFP protein, on the other hand (Fig. 1j, right panel). Altogether, these results suggested that some cellular heterogeneity exists with respect to Hoxd gene transcription in presumptive digit cells, with the possibility for sub-populations of cells to selectively express either one or two genes. Overall, these observations contrasted with the view that all limb cells transcribe all posterior Hoxd genes, a view conveyed by the global analysis of expression patterns by whole-mount in situ hybridization (WISH) and accentuated by schematics published to summarize these expression domains (e.g., ).
Single-limb cell transcriptomics
To have a wider view of this cellular heterogeneity by expanding the analysis to all Hox genes, as well as to see whether it depends on the position and fate of various limb cells, we performed single-cell RNA-seq. Because of its potential to detect as little as single-digit input spike-in molecules, we used the Fluidigm microfluidics C1 captures to obtain the maximal intensity of transcript detection . We enriched for cells expressing at least one Hoxd gene by using only the GFP-positive cells sorted by flow cytometry from the Hoxd1I::GFP mouse E12.5 forelimbs (see Fig. 1d–f). After capture, the cells were sequenced at very high depth to reach the finest sensitivity of gene detection, with an average of about 8.7 M reads per cell (Additional file 2: Fig. S2).
To visualize the relative mRNA contributions from all Hoxd genes, we plotted their cumulative expressions with color-coded single cell (Fig. 2e and Additional file 5: Fig. S4). While the distribution of absolute levels mirrored quite well the pattern previously established using other approaches [4, 6], we observed again a selectivity in expression, which also applied to Hoxd12, Hoxd11, and Hoxd10. Of note, amongst autopod cells positive either for Hoxd13 and/or for Hoxd11, we identified similar proportions as before, with the majority of cells expressing both Hoxd11 and Hoxd13, 40% containing Hoxd13 mRNAs only and 10% with Hoxd11 mRNAs only.
To assess the potential covariances between the five Hoxd genes important for limb development (from Hoxd9 to Hoxd13), we classified by Spearman’s rank correlation the genes that covaried with at least one of the Hoxd genes. A hierarchical clustering from these 76 genes showed a clear segregation between Hoxd11/Hoxd13, on the one hand, and Hoxd9, Hoxd10, and Hoxd12, on the other hand (Fig. 2f). While Hoxd9, Hoxd10, and Hoxd12 were closely associated in either the presence or the absence of their mRNAs, Hoxd11 and Hoxd13 were part of two different sub-clusters associated with different set of genes, suggesting that the cell-specific expression of combinations of Hoxd genes may have some biological relevance.
Combinatorial Hoxd gene expression observed in single limb bud cells
We then performed separate tSNE for autopod and zeugopod cells by clustering cells according to their Hoxd combinatorial patterns (Fig. 3d) and observed that some combinations tend to cluster together. This effect was particularly clear in autopod cells whenever a sufficient number of cells (> 5) were plotted and we noticed that the transcriptional diversity increased along the second dimension of the tSNE, when a higher diversity of Hoxd mRNAs was scored in the same cells. In zeugopod cells, groups of cells also segregated, though not as distinctly, suggesting a more homogeneous distribution of Hoxd mRNAs. These results suggested that sub-populations of autopod cells transcribe various combinations of Hoxd genes.
Analysis of Hoxd cellular clusters
Noteworthy, the clustering of expressed transcripts showed a hierarchical organization with a progression from those cells expressing Hoxd13 only to two, then three, and finally four Hoxd genes (Fig. 4c). As some of these genes were previously identified either as HOX protein targets (e.g., Ppp2ca ) or being part of a Hox functional pathways (e.g., Uty [22, 23] or Hoxa11os [24, 25]), we assessed whether specific target genes could be associated with particular combinations of Hoxd mRNAs. We generated a supervised clustering showing the covariance of known target genes in a Spearman correlation matrix (Fig. 4d, e). When the 199 cells originating from both the autopod and the zeugopod were considered, we found a clear partition of target gene mRNAs into two groups corresponding to the nature of Hoxd mRNAs present (Fig. 4d). The presence of Hoxd9 and Hoxd10 mRNAs aggregated with targets genes such as Hand2 and Sfrp1, whereas Hoxd11, Hoxd12, and Hoxd13 were co-expressed with different target genes such as Ppp2ca and Bmp2/4. Finally, the highest clustering across all cells was observed between Hoxd12, Hoxd13, Dach1, and Lhx9, thus revealing a robust link between these genes.
When only autopod cells were considered, we observed two groups, with Hoxd9 and Hoxd10 transcripts in one cluster, while the more centromeric genes Hoxd11, Hoxd12, and Hoxd13 were transcribed in the other (Fig. 4e). As the former group did not express any of those genes typically upregulated in distal cells (Hoxd12 and Hoxd13), we wondered whether such differences in transcript distribution may reflect various stages in the progression of distal limb cells towards their final fates. We thus implemented a measure of cellular pseudo-age, a strategy that evaluates a temporal hierarchy amongst single cells based on their respective transcriptomes. This approach allows to plot cells along a linearized axis to infer whether the combination alignments observed in the tSNE may correlate with a modulation of the time component [26–28].
During the early stages of limb growth and patterning, limb bud cells absolutely require the expression of Hox genes originating from two distinct clusters, HoxA and HoxD [29, 30]. Here, we describe the single-cell combinatorial expression of Hoxd genes found in cells sorted out by using a Hoxd11::GFP mouse strain. Albeit some cells tend to show higher level of Hoxa genes whenever the levels of Hoxd transcripts were low, this was not the general rule. The fact that we did not score many Hoxa mRNA-positive cells after the enrichment for Hoxd gene expression (Additional file 5: Fig. S4) may however reflect a compensatory mechanism whereby a strong global expression of one cluster would result in the weak transcription of the other. Distinct cellular content for either Hoxd or Hoxa mRNAs could account for the different phenotypic effects of inactivating these genes upon limb morphology, as exemplified by Hoxd13 and Hoxa13 [18, 29, 30]. A more comprehensive single-limb cell sequencing strategy will likely fix this issue.
Our data show that Hoxd quantitative (or “reverse”) collinearity [5, 6]  ought to be considered at a global level since it results in fact from a sum of combinatorial expression of various genes in different cells. We emphasized that in autopod cells, the most frequently expressed gene is Hoxd13, as was expected from previous studies where it was described that this gene is expressed at the highest level, due to its position at the extremity of the gene cluster, i.e., on the centromeric side of a strong TAD boundary [6, 7]. Apart from Hoxd13, other Hoxd genes were more sparsely activated, indicating either a stochastic process or a functional requirement for specific mRNA combinations in different cell types. This heterogeneous cellular situation raises two separate questions; the first concerns the underlying regulatory mechanism, whereas the second has to do with a potential functional significance of these different mRNA combinations.
Different regulatory conformations?
We had previously shown that the regulation of posterior Hoxd genes in the distal limb bud was not implemented exactly in the same manner for all genes. In particular, Hoxd13 is the only gene to be expressed in presumptive thumb cells, the other Hoxd mRNAs being excluded from this very digit [6, 32]. Also, the deletion of the Hoxd13 locus lead to the upregulation of Hoxd12 in thumb cells, yet not of the other remaining genes, suggesting that this thumb-specific expression was associated with the final and most 5′position of the gene on the cluster . The recent identification, in the posterior part of the HoxD cluster, of an unusually high density of bound CTCF molecules may cause this transcriptional selectivity through the use of various sites, taking advantage of their distinct orientations [4, 7, 33].
In this view, the particular orientations of CTCF-binding sites may allow for the transient stabilization of various loop conformations, for example after extrusion driven by the cohesin complex [34, 35]. Accordingly, distinct combinations of posterior Hoxd mRNAs could reflect the formation of specific loop extrusion patterns, in any single cell, as a choice between a fixed number of possibilities determined by the presence of bound CTCF, with some conformations being favored over others. Of note, we found that the cohesin-loading factor Nipbl is strongly downregulated in cells from the Hoxd13 group when compared to cells expressing the full combination (Hoxd10 to Hoxd13). Mutations in this gene have been found in patients with Cornelia de Lange syndrome who have notably a clinodactyly of the fifth finger. Also, a recent report showed that mice heterozygotes for Nipbl display polydactyly and that lower dose of Hoxd11 to Hoxd13 in these mice could further enhance this phenotype .
Other chromatin regulators were found enriched in this list of genes including Jag1, Brd7, Jmjd6, Phf8, Ddb1, Hdac1, Swi5, Smarcc1, Smarce1, Hmgb3, Dnm3os, Cbx1, and Lmnb1 (Additional file 3: Table S1). The product of the latter gene has been associated with the architecture of large domains of inactive chromatin (LADs; ), where the HoxD cluster is not located . Since reduced levels of Lmnb1 gene product have been shown to be associated with reduced expression of polycomb target genes, including the posterior Hoxd genes , its increased expression in those cells containing mRNAs from Hoxd10 to Hoxd13 may reflect a global change in chromatin configuration [37, 40–42]. How would this change relate to a more permissive expression of Hoxd genes, as a cause or as a consequence, remains to be established.
The analysis of single-cell transcriptomes revealed an unexpected hierarchical progression of Hoxd gene expression, from cells expressing a single posterior gene (Hoxd13) to the full combination, from Hoxd10 to Hoxd13. This global transcriptional sequence was inferred from a pseudo-time approach, a method whereby a temporal progression of cells is deduced based on their transcript patterns [26, 28]. We tested this hypothesis using diffusion pseudo-time and found that autopod cells are much more subject to align along a developmental trajectory. This specificity may be associated to the particular way Hoxd genes are regulated in distal limb buds, with a rapid and strong activation of Hoxd13 due to its leading position in the cluster favoring privileged contacts with the various enhancers [2, 43]. It is possible that the recruitment of additional Hoxd genes located nearby may be more progressive, along with local epigenetic modifications, which could be inherited from one cell to its daughter cells. In this view, the number of Hoxd genes expressed would increase along with mitotic divisions leading to the hierarchical progression observed.
Additive cellular or emerging functions?
The second question relates to the potential different functions that limb bud cells may display by carrying distinct combinations of Hoxd mRNAs. The question here is to discriminate between two views of the genotype-phenotype relationship during limb bud development; in a first scenario, each cell would express a determined combination of Hoxd mRNAs, for example in response either to its topological position within the growing limb or to its own “regulatory history,” i.e., the regulations at work in its ancestor cells. In a second scenario, a balanced distribution of cells expressing various Hoxd mRNAs could result from a stochastic distribution of a fixed set of various chromatin architectures . In the former context, the resulting limb phenotype would derive from the additive effect of every single cell, providing one out of the possible sets of information delivered by the various transcriptomes associated. In the second framework, the phenotype would derive from the random mixture of multiple cells expressing distinct transcriptomes with a given balance fixed by the choice of one possible chromatin conformations. Under physiological conditions, these various combinations may allow for differential responses to signaling molecules to generate cellular diversity at the time digit patterns are being established across the limb. It remains to be seen how such potential modulations in the responses to signaling pathways may integrate a theoretical framework whereby digit formation may be an emerging property of the this cellular system [45, 46].
Genetic approaches cannot easily discriminate between these alternatives. In previous studies where the functions of Hoxd genes during limb development were aimed to be assessed separately, various combinations of multiple gene inactivation were used. In most cases however, this consistently led to limited phenotypes due to a fair level of redundancy, particularly amongst Hoxd and Hoxa genes, preventing precise functions to be attributed to specific (groups of) Hox genes (see refs in ). However, the use of multiple gene inactivation revealed that the transcription of Hoxd11 and Hoxd12 contributed functionally and thus added to the mere presence of Hoxd13 transcripts, even though autopods double mutant for Hoxd13 and Hoxa13 would no longer grow and develop digits [12, 29, 43]. This is coherent with our data suggesting that the specific presence of Hoxd11 or Hoxd12 mRNAs is associated with distinct transcriptomes containing additional key regulators of cell fate and chromatin remodeling genes.
Therefore, part of the limb phenotypes observed in Hoxd multiple mutant alleles may result from the different response of a sub-group of cells, which would be differentially impacted by the loss of a given gene. For example, cells that express only Hoxd13 or a combination of Hoxd13 and Hoxd10 mRNAs may not be sensitive to the absence of Hoxd11 transcripts in the corresponding mutant stock. Our results thus stress the necessity to keep in mind the cellular heterogeneity of transcriptional programs even in instances where WISH patterns seem to reveal homogenous distributions of transcripts. In this context, transcript patterns at the single-cell level can help solve the interpretation of genetically deficient phenotypes, even though the co-regulation of Hoxd genes and the functional redundancy of their products make this statement difficult to apply to the present work.
Our results reveal the existence of distinct combinations of Hoxd genes at the single-cell level during limb development. In addition, we document that the increasing combinatorial expression of Hoxd genes in this tissue is associated with specific transcriptional signatures and that these signatures illustrate a time progression in the differentiation of these cells. While this cellular heterogeneity in the combinations of Hox mRNAs may help in understanding the complex transcriptional regulation of these neighbor genes, it will have to be considered when considering the phenotypic outcome of functional studies where one or several such genes were inactivated. Also, further analysis at different developmental stages may enable the reconstruction of the cell fate trajectories and the state transitions that causes the cellular heterogeneity of the early limb bud tissue.
Forelimb tissue samples were isolated from Hoxd11::GFP heterozygous animals at embryonic day 12.5 (E12.5) with day E0.5 being noon on the day of the vaginal plug. The cloning steps for the generation of the Hoxd11 transgenic mice is described in (Additional file 1: Fig. S1). Briefly, the knock-in was done by introducing a bi-cistronic cassette along with an IRES sequence. Hoxd11 was inactivated by the insertion of a TauGFP sequence in a frame into the coding sequence. The BamH1 site was used for insertion of the IRES cassette. The cassette was introduced as a single-copy knock-in (Additional file 1: Fig. S1). The GFP signal detected in this mouse stock reflects the endogenous distribution of Hoxd11 transcription.
E12.5 forelimbs were micro-dissected and fixed with 4% paraformaldehyde for 3 h. Then, the limbs were treated with sucrose at 5, 10, and 15% and then frozen in OCT. Twenty-five-micrometer cryostat sections were dried for 30 min, post-fixed in 4% paraformaldehyde for 10 min, and quenched with 0.6% H2O2 in methanol for 20 min. Slides were then processed using the Ventana Discovery xT with the RiboMap kit. The pre-treatment was performed with mild heating in CC2 for 12 min, followed by protease3 (Ventana, Roche) for 20 min at room temperature. Finally, the sections were hybridized using automated system (Ventana) with a Hoxd13 probe diluted 1:1000 in ribohyde at 64 °C for 6 h. Three washes of 8 min in 2× SSC followed at hybridization temperature (64 °C). Slides were incubated with anti-DIG POD (Roche Diagnostics) for 1 h at 37 °C in BSA 1% followed by a 10-min revelation with TSA substrate (Perkin Elmer) and 10 min DAPI. Slides were mounted in ProLong fluorogold. Images were acquired using a B/W CCD ORCA ER B7W Hamamatsu camera associated with an inverted Olympus IX81 microscope. The image stacks with a 2-μm step were saved as TIFF stacks. Image reconstruction and deconvolution were performed using FIJI (NIH, ImageJ v1.47q) and Huygens Remote Manager (Scientific Volume Imaging, version 3.0.3).
RNA flow cytometry
Double in situ hybridization in single cells for RNA flow cytometry was performed using PrimeFlow RNA (Affymetrix, Santa Clara, CA) reagents following the manufacturer’s protocols. Cell viability was assessed by live/dead fixable dead cell, violet (ThermoFischer; L34955). Hsp90ab RNA probe (a gene expressed ubiquitously) served as a positive control. Hoxd11 and Hoxd13 RNA probes were used for the actual analysis. Cell staining was analyzed on a FACS Astrios located at the EPFL flow cytometry platform. Data analysis was performed by using FlowJoX (Treestar, Ashland, OR). The labelling and flow cytometry were performed on dissociated cells from eight forelimbs obtained from four different animals pooled together.
Single-cell dissociation and fluorescence-activated cell sorting
Pools of embryonic forelimbs obtained from eight embryos at stage E12.5 were dissociated into a single-cell suspension using collagenase from Sigma (collagenase type XI) at 37 °C for 15 min with 10 s trituration. Cells were then filtered on a cell strainer to get rid of clumps. Single cells were then resuspended in FACS solution (10% FCS in PBS with 2 mM EDTA). Fluorescence-activated cell sorting was performed using the MoFlow ASTRIOS EQ cell sorter with a 100-μm nozzle. Through flow cytometry analysis performed using FlowJo (FlowJo LLC ©), we detected 1,602,844 cells positive for GFP in the autopod tissue and 235,000 simply negative. In the zeugopod tissue, 1,527,167 cells were positive, whereas 1,296,068 were negative giving thus a total of 87% GFP-positive autopod cells and 54% positive zeugopod cells.
Control of Gfp expression levels in single cells using RNAscope
To control for the differences between cells positive for the GFP protein and the absence—or low levels—of Hoxd11, we monitored the expression of Gfp mRNA using RNAscope technology (ACD, 320851). Cells were fixed and placed directly on slides following the manufacturer instructions. We used a probe against Gfp mRNA (#400281, C1, as designed by ACD) to assess the number of cells positive for the mRNA after being sorted, based on their fluorescence that reflects only the protein levels of GFP. Images were acquired as five Z stacks on an Axiocam (Zeiss) microscope using a 100X Plan-Neofluar × 100/1.30 Oil objective. 2D projections of the multiple planes were then transformed in mask to count Gfp levels per cells. Automated counting using MATLAB scored 90 positive cells for Gfp mRNA out of 115 cells analyzed (78%).
Single-cell RNA sequencing, library preparation, and mapping
Dissociated single cells were obtained from eight Hoxd11::GFP forelimbs micro-dissected at E12.5 from four littermate embryos. Cells with the highest level of GFP fluorescence (top 20%) were sorted using an Astrios cell sorter with a 100-μm nozzle. Seventy-five-base pair large reads were uniquely mapped to the latest Mus musculus reference genome (mm10) and the ERCC sequences using bowtie2  in a local mode. Raw counts for the annotated ENSEMBL mouse genes (GRCm38) and the ERCC were obtained using the RNA-seq module of the HTSstation portal . The raw counts are summarized in (Additional file 10: Table S3). All single-cell RNA-seq data can be found in the Gene Expression Omnibus (GEO) repository under accession number GSE114748.
Filtering low-quality cells and genes expressed at low levels
Those counts were used to filter out some low-quality cells based on the following criteria: total number of reads mapped > 250, number of genes “expressed” > 2000 (“expressed” = with count > 0), and percent of reads mapped to spike-in sequences < 25%. A total of 199 cells was retained (123 zeugopods and 76 autopods cells). As a control, the positive correlation of expression (Loess regression curve) between the Hoxd11 and the Gfp RNAs is shown (Additional file 11: Fig. S8), tested by Spearman correlation with r = 0.69 in autopod cells and r = 0.49 in zeugopod cells. The Gfp was detected > 4 uniquely mapped reads in 84% of the cells (183 out of 199). Genes expressed at low levels were also removed from the rest of the analysis, and only genes present (raw count > 0) in at least 10% of either the 76 autopods or the 123 zeugopods cells were retained. Hox genes were manually added if they did not satisfy these criteria. A total of 10,948 genes remained. ERCC with null counts through the remaining cells were also excluded from the rest of the analysis. Additional file 12: Table S4 summarizes those criteria (see also Additional file 2: Fig. S2).
Raw counts were normalized with spike-in counts using the R package scran (methods used computeSpikeFactors and normalize version 1.0.4) (http://bioconductor.org/packages/scran/). Prior to normalization, size factors were mean-centered to their batch of origin. An additional normalization step was also applied in order to correct for a potential gene length bias. Additional file 13: Table S5 compiles all the normalized values.
Grouping of Hoxd gene combinations for differential gene expression analyses
HoxD groups were defined per cell and were composed by Hoxd genes with a minimum normalized expression of 5 when count represented at least 5% of the most expressed Hoxd genes in the cell. The differential gene expression analysis was performed with the R package limma (version 3.28.21) . Genes with a minimum absolute log fold change of 2 and a BH-adjusted p value less than 0.01 (false discovery rate (FDR) of 1%) were considered differentially expressed.
The tSNE were computed using the package Rtsne (version 0.13) with the following parameters: two dimensions and a perplexity of 30, a maximum of iterations of 3000, and a seed set at 42. The top highly variable genes (HVG) that were used to plot the tSNE in Fig. 3c, d were selected using the trendVar & decomposeVar methods of the R package scran (version 1.0.4) (http://bioconductor.org/packages/scran/).
Diffusion maps  are tools to analyze single-cell differentiation data. It implements a distance metric relevant to how differentiation data is generated biologically, as cells follow noisy diffusion-like dynamics in the course of taking several differentiation lineage paths . The distances between cells reflect the transition probability based on several paths of random walks between the cells . The analysis was performed using the R package destiny (http://bioconductor.org/packages/destiny).
Network visualization and Venn diagram
The network shown in Fig. 2 was built using weighted interaction networks from various sources of data and is able to process user data into such networks using a system that distinguishes between three different types of user-defined data in its import procedures: real- and binary-valued interaction networks, e.g., physical interaction networks; real-valued gene profile datasets, e.g., multi-sample microarray expression datasets; and binary-valued gene profile datasets . The network shown in Fig. 4 is a summary network of differentially expressed genes that was made with the R package Igraph (version 1.1.2; http://igraph.org).
We thank B. Mangeat for his help with the Fluidigm C1 capture, Jessica Dessimoz for her help with the RNA-FISH procedure, and J. Faget and R. Colisson for their help with the RNA prime flow labelling and flow cytometry analysis. We thank J. Zakany for the help with the production of the knocked in mouse strain and E. Klingler, L. Beccari, and other members of the Duboule laboratories for the discussions. We also thank P. Schwalie and A. Necsulea for their help in the early steps of the single-cell RNA-seq analysis. Cell sorting was performed at the EPFL Flow Cytometry Core Facility.
This work was supported by funds from the Ecole Polytechnique Fédérale (EPFL, Lausanne), the University of Geneva, the Swiss National Research Fund (No. 310030B_138662), and the European Research Council grants SystemHox (No 232790) and RegulHox (No 588029) (to D.D.), as well as the SNF Ambizione grant PZ00P3_174032 (to P.J.F.). Funding bodies had no role in the design of the study and collection, analysis and interpretation of data, and in writing the manuscript.
Availability of data and materials
The datasets generated and analyzed for this study are available in the GEO repository under accession number GSE114748.
PJF and DD conceived and designed the study. PJF carried out the experiments. PJF and ML analyzed the single-cell RNA-seq data. PJF and QLG performed the diffusion pseudo-time analysis. BM and JC generated and maintained the Hoxd11:GFP mouse line. PJF and DD wrote the paper. All authors read and approved the final manuscript.
Ethics approval and consent to participate
All experiments were performed in agreement with the Swiss law on animal protection (LPA), under license No GE 81/14 (to DD).
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Dolle P, Izpisua-Belmonte JC, Falkenstein H, Renucci A, Duboule D. Coordinate expression of the murine Hox-5 complex homoeobox-containing genes during limb pattern formation. Nature. 1989;342:767–72.View ArticlePubMedPubMed CentralGoogle Scholar
- Montavon T, Soshnikova N, Mascrez B, Joye E, Thevenet L, Splinter E, et al. A regulatory archipelago controls Hox genes transcription in digits. Cell. 2011;147:1132–45.View ArticlePubMedPubMed CentralGoogle Scholar
- Andrey G, Montavon T, Mascrez B, Gonzalez F, Noordermeer D, Leleu M, et al. A switch between topological domains underlies HoxD genes collinearity in mouse limbs. Science. 2013;340:1234167.View ArticlePubMedPubMed CentralGoogle Scholar
- Fabre PJ, Leleu M, Mormann BH, Lopez-Delisle L, Noordermeer D, Beccari L, et al. Large scale genomic reorganization of topological domains at the HoxD locus. Genome Biol. 2017;18:149.View ArticlePubMedPubMed CentralGoogle Scholar
- Dolle P, Izpisua-Belmonte JC, Brown JM, Tickle C, Duboule D. HOX-4 genes and the morphogenesis of mammalian genitalia. Genes Dev. 1991;5:1767.View ArticlePubMedPubMed CentralGoogle Scholar
- Montavon T, Le Garrec J-F, Kerszberg M, Duboule D. Modeling Hox gene regulation in digits: reverse collinearity and the molecular origin of thumbness. Genes Dev. 2008;22:346–59.View ArticlePubMedPubMed CentralGoogle Scholar
- Rodriguez-Carballo E, Lopez-Delisle L, Zhan Y, Fabre PJ, Beccari L, El-Idrissi I, et al. The HoxD cluster is a dynamic and resilient TAD boundary controlling the segregation of antagonistic regulatory landscapes. Genes Dev. 2017;31:2264–81.View ArticlePubMedPubMed CentralGoogle Scholar
- Fabre PJ, Benke A, Manley S, Duboule D. Visualizing the HoxD gene cluster at the nanoscale level. Cold Spring Harb Symp Quant Biol. 2015;80:9–16.View ArticlePubMedPubMed CentralGoogle Scholar
- Fabre PJ, Benke A, Joye E, Nguyen Huynh TH, Manley S, Duboule D. Nanoscale spatial organization of the HoxD gene cluster in distinct transcriptional states. Proc Natl Acad Sci U S A. 2015;112:13964–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Zakany J, Fromental-Ramain C, Warot X, Duboule D. Regulation of number and size of digits by posterior Hox genes: a dose-dependent mechanism with potential evolutionary implications. Proc Natl Acad Sci U S A. 1997;94:13695–700.View ArticlePubMedPubMed CentralGoogle Scholar
- Dolle P, Dierich A, LeMeur M, Schimmang T, Schuhbaur B, Chambon P, et al. Disruption of the Hoxd-13 gene induces localized heterochrony leading to mice with neotenic limbs. Cell. 1993;75:431–41.View ArticlePubMedGoogle Scholar
- Delpretti S, Zakany J, Duboule D. A function for all posterior Hoxd genes during digit development? Dev Dyn. 2012;241:792–802.View ArticlePubMedGoogle Scholar
- Kmita M, Fraudeau N, Herault Y, Duboule D. Serial deletions and duplications suggest a mechanism for the collinearity of Hoxd genes in limbs. Nature. 2002;420:145–50.View ArticlePubMedGoogle Scholar
- Dixon JR, Selvaraj S, Yue F, Kim A, Li Y, Shen Y, et al. Topological domains in mammalian genomes identified by analysis of chromatin interactions. Nature. 2012;485:376–80.View ArticlePubMedPubMed CentralGoogle Scholar
- Nora EP, Lajoie BR, Schulz EG, Giorgetti L, Okamoto I, Servant N, et al. Spatial partitioning of the regulatory landscape of the X-inactivation centre. Nature. 2012;485:381–5.View ArticlePubMedPubMed CentralGoogle Scholar
- Andrey G, Duboule D. SnapShot: Hox gene regulation. Cell. 2014;156:856. e1View ArticlePubMedGoogle Scholar
- Spitz F, Gonzalez F, Duboule D. A global control region defines a chromosomal regulatory landscape containing the HoxD cluster. Cell. 2003;113:405–17.View ArticlePubMedGoogle Scholar
- Zakany J, Duboule D. The role of Hox genes during vertebrate limb development. Curr Opin Genet Dev. 2007;17:359–66.View ArticlePubMedGoogle Scholar
- Svensson V, Natarajan KN, Ly L-H, Miragaia RJ, Labalette C, Macaulay IC, et al. Power analysis of single-cell RNA-sequencing experiments. Nat Methods. 2017;14:381–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Sheth R, Barozzi I, Langlais D, Osterwalder M, Nemec S, Carlson HL, et al. Distal limb patterning requires modulation of cis-regulatory activities by HOX13. Cell Rep. 2016;17:2913–26.View ArticlePubMedPubMed CentralGoogle Scholar
- Salsi V, Vigano MA, Cocchiarella F, Mantovani R, Zappavigna V. Hoxd13 binds in vivo and regulates the expression of genes acting in key pathways for early limb and skeletal patterning. Dev Biol. 2008;317:497–507.View ArticlePubMedGoogle Scholar
- Walport LJ, Hopkinson RJ, Vollmar M, Madden SK, Gileadi C, Oppermann U, et al. Human UTY(KDM6C) is a male-specific Nϵ-methyl lysyl demethylase. J Biol Chem. 2014;289:18302–13.View ArticlePubMedPubMed CentralGoogle Scholar
- Copur Ö, Müller J. Histone demethylase activity of Utx is essential for viability and regulation of HOX gene expression in drosophila. Genetics. 2018;208:633–7.View ArticlePubMedGoogle Scholar
- Kherdjemil Y, Lalonde RL, Sheth R, Dumouchel A, de Martino G, Pineault KM, et al. Evolution of Hoxa11 regulation in vertebrates is linked to the pentadactyl state. Nature. 2016;539:89–92.View ArticlePubMedPubMed CentralGoogle Scholar
- Villavicencio-Lorini P, Kuss P, Friedrich J, Haupt J, Farooq M, Turkmen S, et al. Homeobox genes d11-d13 and a13 control mouse autopod cortical bone and joint formation. J Clin Invest. 2010;120:1994–2004.View ArticlePubMedPubMed CentralGoogle Scholar
- Haghverdi L, Büttner M, Wolf FA, Buettner F, Theis FJ. Diffusion pseudotime robustly reconstructs lineage branching. Nat Methods. 2016;13:845–8.View ArticlePubMedGoogle Scholar
- Haghverdi L, Buettner F, Theis FJ. Diffusion maps for high-dimensional single-cell analysis of differentiation data. Bioinformatics. 2015;31:2989–98.View ArticlePubMedGoogle Scholar
- Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 2014;32:381–6.View ArticlePubMedPubMed CentralGoogle Scholar
- Fromental-Ramain C, Warot X, Messadecq N, LeMeur M, Dolle P, Chambon P. Hoxa-13 and Hoxd-13 play a crucial role in the patterning of the limb autopod. Development. 1996;122:2997–3011.PubMedGoogle Scholar
- Kmita M, Tarchini B, Zakany J, Logan M, Tabin CJ, Duboule D. Early developmental arrest of mammalian limbs lacking HoxA/HoxD gene function. Nature. 2005;435:1113–6.View ArticlePubMedPubMed CentralGoogle Scholar
- Nelson CE, Morgan BA, Burke AC, Laufer E, DiMambro E, Murtaugh LC, et al. Analysis of Hox gene expression in the chick limb bud. Development. 1996;122:1449–66.PubMedPubMed CentralGoogle Scholar
- Wagner GP, Vargas AO. On the nature of thumbs. Genome Biol. 2008;9:213.View ArticlePubMedPubMed CentralGoogle Scholar
- Soshnikova N, Montavon T, Leleu M, Galjart N, Duboule D. Functional analysis of CTCF during mammalian limb development. Dev Cell. 2010;19:819–30.View ArticlePubMedPubMed CentralGoogle Scholar
- Fudenberg G, Imakaev M, Lu C, Goloborodko A, Abdennur N, Mirny LA. Formation of chromosomal domains by loop extrusion. Cell Rep. 2016;15:2038–49.View ArticlePubMedPubMed CentralGoogle Scholar
- Haarhuis JHI, van der Weide RH, Blomen VA, Yáñez-Cuna JO, Amendola M, van Ruiten MS, et al. The cohesin release factor WAPL restricts chromatin loop extension. Cell. 2017;169:693–707.e14.View ArticlePubMedPubMed CentralGoogle Scholar
- Lopez-Burks ME, Santos R, Kawauchi S, Calof AL, Lander AD. Genetic enhancement of limb defects in a mouse model of Cornelia de Lange syndrome. Am J Med Genet C Semin Med Genet. 2016;172:146–54.View ArticlePubMedPubMed CentralGoogle Scholar
- Kind J, Pagie L, Ortabozkoyun H, Boyle S, de Vries SS, Janssen H, et al. Single-cell dynamics of genome-nuclear lamina interactions. Cell. 2013;153:178–92.View ArticlePubMedPubMed CentralGoogle Scholar
- Vieux-Rochas M, Fabre PJ, Leleu M, Duboule D, Noordermeer D. Clustering of mammalian Hox genes with other H3K27me3 targets within an active nuclear domain. Proc Natl Acad Sci U S A. 2015;112:4672–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Sadaie M, Salama R, Carroll T, Tomimatsu K, Chandra T, Young ARJ, et al. Redistribution of the Lamin B1 genomic binding profile affects rearrangement of heterochromatic domains and SAHF formation during senescence. Genes Dev. 2013;27:1800–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Reddy KL, Zullo JM, Bertolino E, Singh H. Transcriptional repression mediated by repositioning of genes to the nuclear lamina. Nature. 2008;452:243–7.View ArticlePubMedGoogle Scholar
- Peric-Hupkes D, Meuleman W, Pagie L, Bruggeman SW, Solovei I, Brugman W, et al. Molecular maps of the reorganization of genome-nuclear lamina interactions during differentiation. Mol Cell. 2010;38:603–13.View ArticlePubMedPubMed CentralGoogle Scholar
- Gonzalez-Sandoval A, Gasser SM. On TADs and LADs: spatial control over gene expression. Trends Genet. 2016;32:485–95.View ArticlePubMedGoogle Scholar
- Beccari L, Yakushiji-Kaminatsui N, Woltering JM, Necsulea A, Lonfat N, Rodriguez-Carballo E, et al. A role for HOX13 proteins in the regulatory switch between TADs at the HoxD locus. Genes Dev. 2016;30:1172–86.PubMedPubMed CentralGoogle Scholar
- Giorgetti L, Galupa R, Nora EP, Piolot T, Lam F, Dekker J, et al. Predictive polymer modeling reveals coupled fluctuations in chromosome conformation and transcription. Cell. 2014;157:950–63.View ArticlePubMedPubMed CentralGoogle Scholar
- Sheth R, Marcon L, Bastida MF, Junco M, Quintana L, Dahn R, et al. Hox genes regulate digit patterning by controlling the wavelength of a Turing-type mechanism. Science. 2012;338:1476–80.View ArticlePubMedPubMed CentralGoogle Scholar
- Raspopovic J, Marcon L, Russo L, Sharpe J. Modeling digits. Digit patterning is controlled by a Bmp-Sox9-Wnt Turing network modulated by morphogen gradients. Science. 2014;345:566–70.View ArticlePubMedPubMed CentralGoogle Scholar
- Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–9.View ArticlePubMedPubMed CentralGoogle Scholar
- David FP, Delafontaine J, Carat S, Ross FJ, Lefebvre G, Jarosz Y, et al. HTSstation: a web application and open-access libraries for high-throughput sequencing data analysis. PLoS One. 2014;9:e85879.View ArticlePubMedPubMed CentralGoogle Scholar
- Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47.View ArticlePubMedPubMed CentralGoogle Scholar
- Coifman RR, Lafon S, Lee AB, Maggioni M, Nadler B, Warner F, et al. Geometric diffusions as a tool for harmonic analysis and structure definition of data: diffusion maps. Proc Natl Acad Sci U S A. 2005;102:7426–31.View ArticlePubMedPubMed CentralGoogle Scholar
- Angerer P, Haghverdi L, Büttner M, Theis FJ, Marr C, Buettner F. destiny: diffusion maps for large-scale single-cell data in R. Bioinformatics. 2016;32:1241–3.View ArticlePubMedGoogle Scholar
- Mostafavi S, Ray D, Warde-Farley D, Grouios C, Morris Q. GeneMANIA: a real-time multiple association network integration algorithm for predicting gene function. Genome Biol. 2008;9(Suppl 1):S4.View ArticlePubMedPubMed CentralGoogle Scholar