A multistep computational approach reveals a neuro-mesenchymal cell population in the embryonic hematopoietic stem cell niche Free

The first adult-type hematopoietic stem cells (HSCs), defined by their capability to repopulate the hematopoietic system of adult irradiated recipients, are generated in the Aorta-Gonad-Mesonephros (AGM) region of mid-gestation mouse embryo (Medvinsky and Dzierzak, 1996; Müller et al., 1994). Intra-aortic hematopoietic clusters (IAHCs), containing hematopoietic stem and progenitor cells (HSPCs) closely attached to the ventral side of the dorsal aorta (Tavian et al., 1996; Yokomizo and Dzierzak, 2010; Jaffredo et al., 1998), originate from hemogenic endothelial cells through a process of endothelial-to-hematopoietic transition (EHT) (Jaffredo et al., 1998; Oberlin et al., 2002; de Bruijn et al., 2002; Chen et al., 2009; Zovein et al., 2008; Bertrand et al., 2010; Boisset et al., 2010; Kissa and Herbomel, 2010). In the mouse, IAHCs are also found on the dorsal side of the aorta (Taoudi and Medvinsky, 2007). However, HSCs are preferentially located in the ventral aortic domain (Taoudi and Medvinsky, 2007; Souilhol et al., 2016). Moreover, experiments in the chick embryo revealed that the subaortic tissue originated from splanchnic mesoderm and is essential for IAHC production (Richard et al., 2013). Collectively, these data indicate that the aortic hematopoietic activity is polarized along the dorsoventral axis and that the ventral domain constitutes a supportive HSPC microenvironment. Recent studies have provided novel insight on the AGM microenvironment in terms of phenotype, gene expression and lineage history (Chandrakanthan et al., 2022; Kapeni et al., 2022; Sa da Bandeira et al., 2022; Fadlullah et al., 2022; Crosse et al., 2020; Yvernogeau et al., 2020). At the molecular level, extrinsic factors including cytokines (Robin et al., 2006; Souilhol et al., 2016; Petit-Cocault et al., 2007), signaling pathways (Marshall et al., 2000; Durand et al., 2007; Peeters et al., 2009; Wilkinson et al., 2009; Robert-Moreno et al., 2005; Pouget et al., 2014; Lee et al., 2014; Souilhol et al., 2016; McGarvey et al., 2017), shear stress mediators, proinflammatory and sympathetic nervous system mediators (Adamo et al., 2009; North et al., 2009; Fitch et al., 2012; Espín-Palazón et al., 2014; Li et al., 2014; Mariani et al., 2019; Sawamiphak et al., 2014) participate in HSPC development.

The mesenchyme is primarily responsible for extracellular matrix (ECM) production, organization and remodeling. The AGM mesenchyme remains less well-defined than the hematopoietic component forming the IAHCs or the endothelial lining at the origin of the HSPCs through EHT (Calvanese et al., 2022; Fadlullah et al., 2022; McKinney-Freeman et al., 2012; Solaimani Kartalaei et al., 2015; Baron et al., 2018). Capturing the molecular signature of the AGM mesenchyme in vivo indeed constitutes a technical and biological challenge, as an array of diverse cell types (endothelial cells forming the vascular system, disseminated macrophages and neural crest cells forming the para-aortic lateral primary sympathetic ganglia and the more ventrally located adrenal primordia) are enmeshed within that ECM-embedded mesenchyme. We hypothesize that discrete mesenchymal cell subtypes may coexist in the AGM, depending on their topological association with dorsal (e.g. notochord and somites) and ventral (e.g. gut and gonads) developing tissues.

To investigate this issue, and in line with our previous reports aiming at extracting the molecular core of HSPC niches (Desterke et al., 2020; Charbord et al., 2014), we present here an integrated strategy combining precise laser capture, bulk and single cell transcriptomics with computational analyses including correlation- and information-based gene network analyses. Our analyses indicated the presence of a ventrally located neuro-mesenchymal subset at embryonic day (E) 11.5 within a large mesenchymal cell population. Morphological and lineage tracing studies confirmed the existence of this cell subset as bona fide yet-unexplored population and demonstrated that it was not derived from the neural crest. Moreover, loss-of-function experiments in the zebrafish validated the essential role of the characteristic ECM component decorin (Dcn) for HSPC development. Together, this study provides an unexpected insight onto the transcriptomic landscape of the AGM mesenchyme.

We used a laser microdissection approach to isolate the ventral (VT) and dorsal (DT) aortic tissues from E10.5, E11.5 and E12.5 mouse embryos (Fig. 1A, area delineated in red; Fig. S1A; Table S1; Movie 1) and to compare their gene expression profiles.

The lists of differentially expressed genes (DEGs) between VT and DT were extracted by ANOVA and DEGs were characterized by different gene ontology (GO) profiles. Some of the VT profiles (‘ECM’, ‘metalloproteinase activity’, ‘cell junction’, ‘synapse’, ‘Wnt signaling’, ‘innate immunity’) were of most interest as they peaked at E11.5 (Fig. 1B), a developmental stage after which AGM HSC activity decreases (Kumaravelu et al., 2002). This profile was confirmed using Gene Set Enrichment Analysis (GSEA) (Subramanian et al., 2005) (Fig. 1C,D). The gene sets ‘Complement’, ‘IL6-Jak-Stat3 signaling’ and ‘Interferon alpha response’ were significantly enriched in VT, in agreement with recent studies demonstrating the role of proinflammatory signals on AGM HSCs (Espín-Palazón et al., 2014; Sawamiphak et al., 2014; Li et al., 2014; Mariani et al., 2019). Importantly, the ‘Epithelial-to-Mesenchymal Transition’ gene set was also significantly enriched in VT, including numerous highly ranked genes encoding ECM proteins, such as the interstitial collagen I (Col1a1 and Col1a2), Dcn, a member of the small leucine-rich proteoglycan family with an essential role in collagen fibril assembly, tenascin C (Tnc), an interstitial fibronectin-binding molecule implicated in neuronal guidance, lysyl oxidase homolog 1 (Loxl1), an oxidoreductase essential for collagen crosslinking, and transforming growth factor-beta-induced protein ig-h3 (Tgfbi), a molecule implicated in cell-collagen interactions. Taken together, these analyses revealed a strong enrichment for ECM-associated molecules in the subaortic tissue at E11.5.

Our further study focused on E11.5, as the expression profiles were most contrasted at this time point. To enrich our database in mesenchyme-specific genes, we merged our E11.5 transcriptome with two other expression sets. The first expression set corresponded to flow-cytometry-sorted AGM cells, endothelial cells (CD144^posCD45^neg), differentiated hematopoietic cells (CD34^lowKit^lowCD45^pos) and HSPCs (CD34^highKit^high). Gating strategies and molecular identities of these cell fractions are presented in Fig. S2. The second expression set corresponded to AGM stromal lines derived from the aorta region (AM14-1C4, AM30-3F4) or the urogenital ridges (UG26-1B6 and UG26-3B5) (Oostendorp et al., 2002; Charbord et al., 2014). Using ANOVA and then intersecting the resulting DT and VT DEGs led to a list of 2401 genes: 1199 upregulated in DT and 1202 in VT (Table S2). Using Principal Component Analysis (PCA) (Fig. 2A) allowed us to discriminate DT from VT according to PC1, and endothelial cells (EC) from differentiated hematopoietic cells (HC) according to PC2 (34.2 and 19.1% of the variance, respectively).

To find out which genes were specific to each of the four discriminated tissue/cell populations (VT, DT, EC, HC) and to reveal their relationships, Weighed Gene Correlation Network Analysis (WGCNA) combined with PCA was applied to the set of 2401 genes (Langfelder and Horvath, 2008). PC1 and PC2 corresponded to the external traits (DT-to-VT and EC-to-HC contrasts), and the merged dynamic tree cut method was used to identify the modules (Fig. 2B, Fig. S1B). The distribution of module eigenvectors indicated that two modules (designated blue and yellow) corresponded to VT, another two (brown and turquoise) to DT, one (cyan) to EC and two (magenta and salmon) to HC (Fig. 2C).

Gene analysis was then performed on the seven modules (see Table S2 for gene allocation). Firstly, it confirmed the endothelial and hematopoietic identities of the populations corresponding to the EC and HC modules, respectively. Secondly, it indicated that the DT and VT modules had very distinctive features. Most genes of the DT module designated ‘brown‘ were implicated in neurogenesis; many transcription factors essential for neuronal development and differentiation such as Sox2, Pou3f3 and Npas3 were among the most connected genes. Most genes of the DT module designated ‘turquoise’ were implicated in chondrogenesis; the transcription factor Pax1, which is expressed in pre-chondrocyte condensations, and the ECM cartilage collagen components Col2a1 and Col9a1 were among the most connected genes. The VT module designated ‘yellow’ was poorly defined because it included many unexplored and predicted genes (olfactory receptors, Riken, Gm genes). In line with the pairwise DT to VT comparison (as shown on Fig. 1B), the VT module designated ‘blue’ was of most interest because it included genes implicated in either adhesion-related processes (‘cell adhesion’, ‘basement membrane’, ‘ECM’, ‘collagen’, ‘homophilic cell adhesion’, ‘focal adhesion’) or neuronal-related ones (‘synapse’, ‘ion channel’, ‘GABA-A receptor activity’, ‘neuroactive ligand-receptor’, ‘neurotrophin signaling’) (Fig. 2D). Among the most connected genes were the ECM components Dcn and Tnc, the neurotrophin-3, glutamate and netrin receptors Ntrk3, Grid1 and Unc5c, and molecules implicated in synapse functioning and/or neuronal development such as Clstn2, Amph, Syt1 and Cntn4 (Fig. 2E,F). Genes encoding for soluble factors such as Bmp6, Gdf6, Csf1, Fgf14, Igfbp1 and Igfbp5, Ntf5 and Postn also characterized the ‘blue’ module (Fig. 2E).

These data suggested the presence of three cell populations within the dorsal and ventral tissues, two dorsally-located (above the aorta ceiling), either expressing neuronal or chondrocytic genes, and one ventrally-located (beneath the aorta floor) expressing both adhesion- and neuronal-related genes. To find out how those cell populations would integrate within the AGM region we performed single cell transcriptomic analysis on the entire E11.5 AGM.

The entire AGM was dissected as shown on Fig. 1A (area delineated by a green line). After quality control (Fig. S3A) we retained 15,154 out of 15,862 cells, with 2698 genes detected on average. After normalization and dimension reduction, the Louvain algorithm allowed us to define clusters projected on Uniform Manifold Approximation and Projection (UMAP) embedding (Fig. 3A). Thirteen of the 22 clusters were packed together forming a large cluster, clearly distinct from the nine others that presented as single entities.

The large cluster consisted of mesenchymal cells as indicated by the widespread expression of the interstitial collagens Col1a1, Col1a2 and Col3a1, the lysyl oxidase homolog Loxl1, the microfibril component Fbn2, the ECM-degrading metalloproteases Mmp16 and Bmp1, the proteoglycan Vcan, several molecules implicated in matrix assembly (Ccdc80), matrix organization (Fras1) and cell-collagen interactions (Tgfbi), the receptors for the growth factors of mesenchymal cells (Pdgfra and Pdgfrb), the hematopoietic chemokine expressed by stromal cells (Cxcl12), and finally molecules implicated in EMT, the transcription factor Snai2 and the growth factor Ptn (Fig. 3B; Fig. S3B). Of note, we had found in a previous work that some of these genes (Ccdc80, Snai2 and Tgfbi) were regulators of HSPC development in the AGM (Charbord et al., 2014).

To find out the cell identities of the different clusters aggregated in that large mesenchymal population, we first examined the first-ranked genes. Among the eight major clusters (C0 to C7), each containing more than 600 cells, five corresponded to well-defined identities: cells in C0 were chondrocytic, in C1 and C5 urogenital, in C2 cycling and in C4 smooth muscle (the characteristic genes are indicated on the upper table of Fig. 3A, and the distribution of the first-ranked gene is shown on Fig. 3C). For the three remaining clusters (C3, C6 and C7) the first-ranked genes did not provide any clear clue (upper table of Fig. 3A and C).

Our second approach was therefore to investigate in which clusters the genes of the two DT (turquoise and brown) and the one VT (blue) modules indicated by WGCNA were expressed. Genes of the VT module, but of only one of the DT modules (the ‘turquoise’) projected onto clusters of the mesenchymal population. The projection of the genes of the ‘turquoise’ DT module onto the cluster map indicated that they were expressed in the C0 cluster, with a minor contribution in C2 (Fig. S3C); these results confirmed the chondrocytic identity of C0 cells. Similar projection of the genes of the VT module onto the cluster map indicated that they were expressed in the clusters C6 and C7 (Fig. 3D). Hence, these data suggested the existence of a specific mesenchymal cell population, that we called neuro-mesenchyme as it expressed both adhesive (such as Dcn and Tnc) and neuronal (such as Ntrk3) genes, represented in the C6 and C7 clusters (71% of the genes expressed in C6 were also expressed in C7 at about the same intensity level, suggesting that cells in C6 and C7 clusters belonged to the same population).

The remaining large cluster C3 was then examined again by taking into account the distribution of genes also expressed in other clusters. It revealed that many genes expressed in C6/C7 (such as Ntrk3, Tnc, Hlx, Col25a1), in C0 (such as Col2a1, Col9a1, Foxf1) or in C1 (such as Gata2) were also expressed in C3 (Fig. 3D; Fig. S3C,D), indicating that the C3 cluster was a mix of several cell populations.

Finally for the five minor mesenchymal clusters (containing fewer than 600 cells), there were few (C8-C11) or no (C12) significant genes; C8 expressed some genes implicated in ECM organization or synthesis, and C8-C11 expressed transcription factors implicated at different developmental stages (upper table of Fig. 3A). Many genes implicated in the cell cycle were intensely expressed in C12.

The first ranked genes of the non-mesenchymal clusters revealed without ambiguity their identity (lower table of Fig. 3A). There were two neuronal clusters C13 and C14, the latter including many genes characteristic of neural crest cells (NCC); interestingly, the genes of the DT ‘brown’ module found using WGCNA were expressed in the non-NCC C13 cluster. There was also one endothelial cell cluster (C15) and one containing innate immunity cells (C16), onto which the genes of the endothelial (‘cyan’) and hematopoietic (‘magenta’ and ‘salmon’) modules, respectively, were projected (Fig. S4). Finally, there was a very small cluster (C19) of germinal cells.

Taken together, the single cell analysis confirmed both the heterogeneity of the AGM region taken as a whole, and that of the AGM mesenchyme that included 83% of all AGM cells. The ventrally-located AGM mesenchyme comprised two major populations, one of prospectively ‘neuro-mesenchymal’ cells expressing both adhesion and neuronal genes (C6 and C7 clusters), and the other of urogenital cells (C1 and C5). The dorsally-located AGM mesenchyme comprised one major population of chondrocytic cells (C0). Dorsal versus ventral allocation was meaningless for C4 containing smooth muscle cells located around the aorta, and for C2 containing cycling cells that could be present among all mesenchymal cells. Finally, the location was less clear for C3 that contained a mix of cells that could be ventrally- or dorsally-located, and for the minor clusters C8 to C12 in which genes were too few to allow precise spatial allocation.

The previous analyses strongly suggested that clusters C6 and C7 represented a population of mesenchymal cells expressing neuronal genes. To further demonstrate the existence of such population, we extracted from the ‘blue’ VT module list two gene sets representative of the neural and mesenchymal lineage (Table S3). Gene set scoring using the Seurat AddModuleScore function indicated that 63.9% of the cells in C6 and C7 were positive for both gene subsets (Fig. 4A).

Contrary to the distribution of genes representative of the neural subset that was restricted to C6 and C7, that of the mesenchyme subset was more extended, in particular to cells belonging to C1 (Fig. 4A, third panel from the left). By searching for individual genes expressing such a profile we found that it corresponded to several ECM components: Dcn, and collagen family members Col4a5, one of the major components of basement membranes, and Col6a3 and Col6a2, two chains of the cell-binding protein collagen VI (Fig. 4B).

To get more insight on the organization of genes characterizing the neuro-mesenchyme and to identify direct gene-to-gene interactions and potential causal links, we applied the multivariate information-based inductive causation (miic) algorithm (Verny et al., 2017; Sella et al., 2018) on the VT module genes. To yield significant results the miic algorithm requires a large number of observations compared with variables. We therefore restricted the analysis to the central core of the module including 98 genes and we considered all cells belonging to the large mesenchymal cluster (N=12,516). After filtering out low confidence links, 89 nodes defining direct high-confidence paths were retained (Fig. 4C). The resulting network was pruned compared with that provided by WGCNA. It included four major hubs (with degree ≥17) – Unc5c, Col1a1, the neural cell adhesion molecule Ncam1, and the transcription factor Zmiz1 – and a few less prominent ones, such as the dopa decarboxylase Ddc. Of interest, Ddc and all its eight links (including Tnc, Kcnip4, a channel interacting protein, Nkain2, an ion transport ATPase and Adam23, a disintegrin without catalytic activity) were expressed only in C7, indicating that their mutual interaction links corresponded to a specific gene distribution given by the single cell analysis (Fig. 4C, right panel). Importantly, all of the ECM components were directly linked to neuronal molecules (e.g. Col1a1-Ntrk3-Tnc, Dcn-Ncam1-Tnc, Ntrk3-Col8a1-Ddc, Ntrk3-Tnc-Ddc, Ncam1-Tnc-Ddc), and nine edges were confirmed by the analysis of protein-protein interactions using STRING (Szklarczyk et al., 2019) (Fig. 4D). Finally, the miic algorithm uncovered a number of ‘v-structures’ (e.g. the nuclear receptor Nr1h5 → Tnc ← Ncam1), suggesting cause-and-effect interactions to be investigated in future works.

Finally, to investigate whether this cell population was present only during E11.5 HSPC production, we used the WGCNA ‘blue’ module as a gene set of reference to interrogate the global transcriptomes of ventral micro-dissected tissues at E10.5 and E12.5, as well as the transcriptomes of aortic mesenchyme and urogenital ridge stromal lines. GSEA indicated that the neuro-mesenchymal signature was enriched in E11.5 ventral tissues compared with E10.5 and E12.5 ventral tissues, suggesting that the neuro-mesenchyme is mostly present during HSPC production in the AGM (Fig. S5, left and middle panels). In addition, we observed that the ‘blue’ module is statistically enriched in aorta mesenchymal compared with urogenital ridge stromal lines, confirming the aortic mesenchymal origin of the micro-dissected tissues (Fig. S5, right panel).

Taken together, the data pointed to the existence of a, yet undescribed, E11.5 neuro-mesenchymal population, in which the majority of individual cells express both neuronal and mesenchymal genes that appear to be directly linked to one another in a network with a few highly connected hubs.

To validate the in-silico data, we analyzed the distribution of three characteristic proteins on AGM sections. We selected three candidates, tenascin (Tnc), decorin (Dcn) and the NT3 growth factor receptor (Ntrk3), which we considered as characteristic of the neuro-mesenchyme as they preeminently appeared as highly connected genes in WGCNA, as expressed in the C6 and C7 clusters in scRNA-seq, and as nodes of ECM-to-neuronal links in miic. In addition, Tnc, Dcn and Ntrk3 expression profiles were mapped exclusively onto the ‘stroma’ clusters of the human week 4 AGM single cell transcriptome atlas (Calvanese et al., 2022); interestingly, Dcn was expressed in the three stroma clusters, whereas Tnc and Ntrk3 expression was confined to cells of the larger cluster, many of which expressed both genes (Fig. S6).

Dcn expression was detected by immunofluorescence around the notochord and the aorta, and in the subaortic tissue at E11.5, but not in the aortic endothelial and vascular smooth muscle cells, or macrophages, as shown by the absence of overlapping staining between Dcn and CD31, α-smooth muscle actin (α-SMA) or allograft inflammatory factor 1 (Iba1; Aif1) (Fig. 5A,C; Fig. S7A). Tnc was also highly expressed in cells and associated ECM in large areas of the subaortic tissue (Fig. 5A; Movie 2). Importantly, in agreement with previous report (Marshall et al., 1999), a similar distribution of Dcn and Tnc was observed in human embryos at day 30 of gestation (Fig. 5B). Ntrk3 was present, as expected in the dorsal root ganglion cells, and in cells in extensive areas around and underneath the dorsal aorta (Fig. 5C; Movie 3). Importantly, the distribution of Dcn, Tnc and Ntrk3 staining was quite similar in large areas of the subaortic tissue. The expression of these proteins was clearly distinct from that of cells expressing tyrosine hydroxylase (Th), an enzyme involved in catecholamine biosynthesis, and of axons expressing the tubulin β3 chain (Tubb3), a neuronal specific tubulin expressed in neurons of the peripheral and central nervous systems (recognized by the antibody Tuj1) (Fig. 5C; Fig. S7B-D); immunostaining for Th revealed discrete foci in the dorsal and ventral domains of the aorta, corresponding to primary sympathetic ganglia and adrenal primordia (Fig. 5C); immunostaining for Tuj1 revealed a dense array of axons on the lateral sides, and very few foci underneath the aorta (Fig. 5C; Fig. S7B-D). Moreover, sagittal analyses of E11.5 3D embryo imaging revealed a strong expression of Dcn and Ntrk3 in the middle part of the AGM known to be enriched for HSC activity (Fig. S8) (Easterbrook et al., 2019; Mascarenhas et al., 2009). Finally, transmission electron microscopy revealed the presence of transverse and longitudinal ECM fibers in the intercellular space between subaortic cells at E11.5 (Fig. 5D).

These data confirmed the presence of a subaortic mesenchyme in which the ECM components Dcn and Tnc had distribution similar to that of the cells expressing NTRK3. Moreover, these cells were clearly distinct from neurons expressing Tuji or Th. These results therefore confirmed the presence of a yet unidentified neuro-mesenchymal population.

An obvious issue was whether this neuro-mesenchymal cell population was derived from NCCs. The computational analysis was not in favor of such an assumption, as genes characteristic of NCCs were expressed by cells belonging to the non-mesenchymal C13 and C14 clusters. Indeed, the major NCC marker Sox10, the early (Phox2b) and late markers (Phox2a and Th) of the sympathoadrenal compartment were expressed in C14 and C13 (Fig. 6A). To definitely exclude this hypothesis we took advantage of a Human tissue Plasminogen Activator (HtPA)-Cre mouse labeling neural crest cells and their derivatives in vivo (Pietri et al., 2003). After crossing those mice with GT Rosa Tomato mice, E11.5 embryos were collected and analyzed for the expression of Tomato, Dcn, Tnc and Ntrk3 (Fig. 6B-D). Tomato⁺ cells were detected as expected on the dorsal, lateral and ventral sides of the aorta (Fig. 6B-D). There was no overlap between Tomato expression and that of Dcn, Tnc and Ntrk3 (Fig. 6B-D). These data therefore allowed us to exclude the hypothesis that the neuro-mesenchymal population was NCC-derived.

To functionally validate our computational analysis, we performed in vivo knockdown experiments in zebrafish embryos using antisense morpholinos (MOs). The zebrafish is indeed an excellent system to test candidate gene function in a robust and rapid manner (McKinney-Freeman et al., 2012; Charbord et al., 2014; Yvernogeau et al., 2020; Fadlullah et al., 2022). Dcn was selected as a potential HSPC regulator as it is an extracellular matrix component expressed in the subaortic mesenchyme both in murine and human embryos, and, to our knowledge, not reported yet as playing a role in developmental hematopoiesis.

We first investigated dcn expression in zebrafish embryos. In agreement with a previous study (Zoeller et al., 2009), we detected dcn expression in the developing head and trunk (particularly the somites) using whole mount in situ hybridization (WISH) (Fig. S9A). In line with these observations, it has been recently reported that dcn was preferentially expressed in somite-derived endothelial cells acting as a vascular niche for hemogenic endothelial cells (Sahai-Hernandez et al., 2023). To refine our analysis, using RT-qPCR we evaluated dcn expression on sorted cell populations: flk1⁻cmyb⁻ (non-hematopoietic non-endothelial mesenchymal cells), flk1⁻cmyb⁺ (myeloid progenitors), flk1⁺cmyb⁻ (vascular endothelial cells) and flk1⁺cmyb⁺ (HSPCs) from the trunk of 48 h post-fertilization (hpf) embryos using the flk1:mcherry;cmyb:gfp zebrafish line (Fig. S9B). Our results revealed that dcn was significantly more expressed in flk1⁻cmyb⁻ and flk1⁺cmyb⁻ cells in the zebrafish embryo (Fig. S9C).

The formation of HSPCs in the dorsal aorta was monitored by WISH for the expression of the early (runx1) and late (cmyb) HSPC markers (Fig. 7A,B). The morphants dcn injected with MO^START exhibited strong hematopoietic defects in the dorsal aorta, as shown by significant reduced expression of runx1 and cmyb in morphants compared with respective controls at 30 hpf and 36 hpf, respectively (Fig. 7A,B). Similar phenotypes were observed with the injection of a MO^SPLICE (Fig. S10A, S10B).

To confirm the knockdown experiments and to exclude any ‘off-target’ effect of MOs, we examined zebrafish mutants. As adult dcn^−/− zebrafish (dcn^sa20247 line), harboring a point mutation that induces a stop codon in the dcn coding sequence, were viable and fertile, we crossed dcn^−/− females with dcn^−/− males. Resulting embryos exhibited a significant decrease in the numbers of runx1- and cmyb-positive cells compared with dcn^+/+ embryos from wild-type sibling crosses (Fig. 7C,D), thus confirming the phenotypes of the dcn morphants.

Confocal microscopy using kdlr:mCherry;gata2b:GAL4;UAS:eGFP transgenic animals showed a strong reduction in the number of HSPCs in the dorsal aorta in dcn morphants (Fig. S10C). Importantly, the vasculature of flt1:dtomato;flt4:yfp embryos was intact, indicating that the reduction in HSPCs in dcn morphants was not related to vascular abnormalities (Fig. S10D).

As Dcn was known to bind to Tgfβ1 and to serve as a negative-feedback regulator of Tgfβ signaling (Yamaguchi et al., 1990), we investigated by RT-qPCR the expression of Tgfβ ligands and receptors in dcn morphants and control embryos. We did not detect any significant differences in the expression of tgfb2, tgfb3, tgfr1a and tgfr2 (Fig. S9D). In contrast, tgfb1a expression was significantly reduced in dcn morphants compared with controls (Fig. S9D). These data suggest that the effect of dcn on AGM HSPCs may be related to a deregulation of Tgfβ1 expression, which is required for HSPC emergence (Monteiro et al., 2016).

Finally, as the ‘blue’ module contains several genes encoding for soluble factors, we also selected the Growth Differentiation Factor 6, Gdf6, for functional validation. Interestingly, gdf6a (but not gdf6b) was reported to be expressed in the trunk and around the dorsal aorta in the zebrafish embryo (Krispin et al., 2018). We observed that gdf6a morphants exhibited a significant decrease in the number of runx1-positive cells in the dorsal aorta compared with uninjected controls (Fig. S11). The effect of gdf6a on aortic HSPCs could be indirect, as gdf6a is required for vascular integrity by restraining Vegf signaling (Krispin et al., 2018; Hall et al., 2002).

In this work, we have uncovered at the molecular level the cell identities of the AGM mesenchyme at the time of HSPC formation in the mouse embryo. Our computational strategy consisted of several steps: first, using bulk transcriptomics and additional datasets, subtracting the genetic contribution of endothelial and hematopoietic cells of the sub-dissected aortic tissues; second, using a correlation-based algorithm, unraveling the gene networks corresponding to the ventral and dorsal tissues; third, using single cell transcriptomics, projecting the gene modules onto the Louvain map to characterize the clusters that could not be defined by highly-ranked gene allocation; and fourth, using an information-based algorithm, pruning the gene network corresponding to the cell population of most interest to better discern crucial gene-to-gene links. This computational analysis allowed us to reveal the presence, within a large and heterogeneous mesenchymal cell population, of a ventrally-located neuro-mesenchymal subpopulation, which was confirmed by immunofluorescence studies using antibodies against Dcn, Tnc and Ntrk3, three proteins corresponding to genes Dcn, Tnc, and Ntrk3 that appeared preeminently throughout the in-silico analysis. Moreover, lineage tracing studies indicated that the neuro-mesenchymal population was not derived from the neural crest. Finally, loss-of-function experiments in 30- and 36-hpf zebrafish confirmed the essential role of the characteristic ECM component Dcn for HSPC development.

The functional association between the hematopoietic and neuronal systems has already been documented, both in the adult and embryonic contexts. For example, in the adult bone marrow, stromal cells expressing the β3-adrenergic receptor were shown to be sensitive to noradrenaline rhythmically released by the sympathetic nervous system; this intercellular connection pathway results in the circadian oscillations of Cxcl12, a major extrinsic factor involved in HSPC maintenance and homing (Méndez-Ferrer et al., 2008). In the embryo, it has been reported that Gata3 controls the number of AGM HSCs through the production of catecholamines, showing a strong influence of the sympathetic nervous system on AGM HSCs (Fitch et al., 2012). Catecholamines may exert their effect directly on AGM HSPCs, as the β2-adrenergic receptor expressed by HSPCs is required for HSC production (Fitch et al., 2012; Kapeni et al., 2022). In the zebrafish, trunk NCCs were shown to interact with precursors of aortic HSPCs and to contribute to their supportive microenvironment (Damm and Clements, 2017). In mice, subaortic mesenchymal cells may have a dual origin (mesodermal- and neural crest-derived) (Kapeni et al., 2022; Chandrakanthan et al., 2022) with Pdgfra⁺ mesodermal stromal cells being replaced by NCC-derived cells at E13.5, a time coinciding with the loss of EHT (Chandrakanthan et al., 2022). Remarkably, in our present work we identified at E11.5 a mesenchymal cell subpopulation, not derived from NCCs, containing both neuronal genes, implicated in synapse formation, membrane channel activity and neurotrophin signaling, and mesenchymal genes, in particular some distinctive ECM components. This molecular signature was enriched in the subaortic tissue at E11.5, compared with E10.5 and E12.5. The neuro-mesenchymal cell population may represent a transitory state between a mesoderm-derived mesenchymal state and a neuroectoderm-derived one (Chandrakanthan et al., 2022; Richard et al., 2013). Moreover, its establishment underneath the dorsal aorta may be connected to the concomitant development of the nervous system and the projection of neuronal fibers on the lateral and ventral sides of the dorsal aorta at mid-gestation. Collectively, our findings prompt us to hypothesize that the sympathetic nervous system and the subaortic neuro-mesenchyme may cooperate to regulate HSPC development in the dorsal aorta. The dual molecular identity of the subaortic mesenchyme suggests that these cells may eventually be responsive to signals delivered by the sympathetic nervous system.

So far, very little is known about the composition and functions of the ECM in the AGM. Tnc is known to be expressed underneath the dorsal aorta both in the mouse and human embryo (this paper and Marshall et al., 1999) and to promote the differentiation of hematopoietic and endothelial precursors from human pluripotent stem cells (Uenishi et al., 2014). Matrix metalloproteinases, such as Mmp2 and Mmp9, have been shown to regulate EHT and the migration of aortic HSPCs towards secondary hematopoietic tissues (Travnickova et al., 2015; Theodore et al., 2017). Interestingly, in bone marrow HSPC-supportive stromal cells Dcn is not upregulated but is downregulated compared with less or non-supportive stromal cells (Desterke et al., 2020), which is in agreement with its reported role as an HSPC negative regulator in adult mice (Ichii et al., 2012). These data underscore the context-dependency of some of the mediators: activators in the embryonic niche, but inhibitors in the adult one. Dcn could exert its effects on AGM HSPCs as a structural, basement membrane ECM component interacting with other ECM proteins (such as laminins, nidogens, fibronectins or collagens), and/or as a substrate for cytokines such as TGFβ1 or other soluble mediators, such as those widely known for the bone marrow HSPC microenvironment. Our investigations in zebrafish embryos suggest that the reduction in HSPCs in dcn morphants was not related to vascular abnormalities. However, as dcn was significantly more expressed in flk1⁻cmyb⁻ and flk1⁺cmyb⁻ trunk cells, we cannot rule out that the loss-of-function phenotype for dcn could also be initiated in endothelial cells.

In conclusion, we identified a unique non-NCC derived mesenchymal cell population harboring both an adhesive and neuronal molecular identity within the subaortic tissue. Our data also highlight the role of a pre-eminent ECM component as part of the HSPC-supportive niche in the embryo.

C57BL/6 mice were bred at Janvier Labs and maintained according to institutional guidelines in the animal facility of the Laboratory of Developmental Biology. Mouse models for lineage tracing were the Tg(PLAT-cre)116Sdu, referred to as (HtPA)-Cre (Pietri et al., 2003), and Cre reporter mice B6.Cg-Gt(ROSA)26Sortm9(CAG-tdTomato)Hze/J (Madisen et al., 2010). Pregnant females were killed by cervical dislocation at appropriate time of gestation, with the day of vaginal plugging taken as E0.5. Uteri were taken and placed in 10% fetal bovine serum/1% Penicillin/Streptomycin (Gibco)/phosphate-buffered saline (PBS) (Invitrogen). The embryonic stages were confirmed through morphological analysis.

Human embryonic tissues were obtained following voluntary abortions. Informed consent for the use of the embryo in research was obtained from the mother, and tissue collection and use were performed according to the guidelines and with the approval of the French National Ethics Committee.

We used the Arcturus Laser Capture Microdissection System to isolate the ventral and dorsal aortic tissues surrounding the dorsal aorta from mouse embryos (Fig. 1A area delineated in red). The isolated dorsal tissue was localized in between the notochord and the roof of the aorta and the isolated ventral tissue was underneath the aorta. Unfixed embryos were embedded in OCT medium in presence of isopentane cooled at −65°C. Then 20 µm cryostat sections were prepared and stained before laser micro-dissection. Three to eight micro-dissected tissues, corresponding to a total surface of between 84,431 and 412,076 µm², depending on the developmental stages as indicated in Table S1, were collected on appropriate caps and then RNA was extracted using Acturus Picopure RNA isolation kit. RNAs were eluted in a final volume of 13 µl and their concentrations and integrity quantified on a bioanalyzer (Agilent). For each developmental stage, three to four RNA samples corresponding to DT and VT were prepared from independent embryos and prepared for bulk transcriptomics. Between 200 and 1000 pg of total RNA was reverse transcribed following the Ovation Pico V2 system (Nugen) and the resulting double strand cDNA was used for amplification based on SPIA technology. After purification according to the Nugen protocol, sense target DNA was fragmented and biotin-labeled using Encore Biotin Module kit (Nugen). After control of fragmentation using Bioanalyzer 2100, cDNA was hybridized to GeneChip^® MouseGene 2.0 ST (Affymetrix). Normalization and functional analyses of the datasets were performed using Partek Genomics Suite software.

DT versus VT comparisons were made using a two-way ANOVA with a low stringency P-value (P<0.06) with ‘Polarity’ and ‘Embryo’ as grouping variables for E11.5 and E12.5 samples, or three-way ANOVA with ‘Polarity’, ‘Embryo’ and ‘Series’ (corresponding to two independent transcriptomic datasets) for E10.5 samples. After incorporation of additional transcriptomes (AGM stromal cell lines or AGM sorted cell populations) DT versus VT comparisons were made using two-way ANOVA with ‘Cell/Tissue identity’ and ‘Series’ (the different batches) as grouping variables.

Gene ontology (GO) categories were analyzed using DAVID (https://david.ncifcrf.gov) (Huang et al., 2007).

GSEA aims at identifying the significant genes within a gene expression set contrasting two phenotypes (DT versus VT). The expression set is compared with pre-established datasets corresponding to a number of collection profiles, in this work the ‘hallmark’ database. Datasets are retained if the distribution of their genes is shifted toward either phenotype and the statistical relevance of the shift is established by random permutations (here 1000) of the phenotypes (Subramanian et al., 2005).

To reconstruct gene network and identify specific modules for the different cell populations, we used WGCNA (Langfelder and Horvath, 2008). The chosen values for essential input parameters were: (1) the power β obtained from the scale-free fit model curve (r²>0.8) – the β power of 18 was selected such that the connectivity distribution fitted a scale-free model (Fig. S1); (2) the Pearson's correlation; (3) the signed adjacency; (4) the average linkage hierarchical clustering; (5) the minimum module size of 30. Cytoscape (https://cytoscape.org) was used to visualize the connectivity plot of the most connected genes. Force-oriented layout was used for data representation.

The miic algorithm combines constraint-based and information-theoretic frameworks to reconstruct causal and non-causal networks from large-scale datasets (Verny et al., 2017; Sella et al., 2018). Starting from a fully connected graph, miic iteratively removes dispensable edges, by uncovering significant information contributions from indirect paths. Miic also provides an edge-specific confidence assessment of retained edges, which are oriented based on the signature of causality in observational data.

Seven E11.5 AGMs from three litters were dissected (Fig. 1A, area delineated in green), digested with collagenase and submitted to 10x Genomics Chromium next GEM single cell 3′ kit (version 3.1) (two technical replicates). Only cells expressing 200-7000 genes with a percentage of mitochondrial RNA under 15% were kept. In addition, genes expressed in less than three cells were excluded from the analysis. As no batch effect was observed between the replicates, the two datasets were merged without correction. ScTransform normalization (Hafemeister and Satija, 2019) was performed using the Seurat package v3 (Stuart et al., 2019) with regression of the cell cycle. To regress the cell cycle genes, the default human cell cycle gene list given by Seurat was converted into mouse genes using the biomaRt package (Durinck et al., 2009) and was used following the default recommendation for cell cycle regression with scTransform as described in the Seurat vignette. Clustering was performed with the Louvain algorithm with a resolution of 0.7 on the first 50 PCA components and the cells were projected on a UMAP embedding based on the first 42 PCA components. We then identified the DEGs, selecting genes with a fold-change>0.25 in each cluster compared with the whole dataset using the Wilcoxon sum rank test. Only genes expressed in more than 25% of the cluster cells and with a q-value<0.05 were kept. These genes were then ranked and used to determine each cluster identity. For the gene scoring, we defined two lists of mesenchymal and neural genes belonging to the WGCNA blue VT module and used the Seurat AddModuleScore function. This function returns positive values if the genes from the module display a higher expression than control genes selected at random among the ones with a similar average expression over the dataset. The miic algorithm (Verny et al., 2017; Sella et al., 2018) was used to analyze the genes of the blue VT module, taking into account all 12,516 cells belonging to the large central aggregate defined by the Louvain clustering (C0-C12) and the raw counts for each gene of the VT ‘blue’ module core (89 genes).

E11.5 AGMs were dissociated in 0.125% collagenase type I (Sigma-Aldrich) in PBS supplemented with 10% fetal calf serum (FCS) for 20 min at 37°C. Cells were dispersed by gentle pipetting, washed twice, filtered through a 70 µm nylon cell strainer (BD Falcon) and resuspended in Dulbecco's modified medium (DMEM) (Invitrogen) supplemented with 10% FCS. Cells were labeled with APC anti-mouse CD34 (clone RAM34, eBiosciences, 1:100), PE anti-mouse CD117 (clone 2B8, BioLegend, 1:100) and FITC anti-mouse CD45 (clone 30F-11, BioLegend, 1:100) or with PE anti-mouse CD144 (clone BV13, BioLegend, 1:100) and APC anti-mouse CD45 (clone 30F-11, BioLegend, 1:100). Viability was determined using 7-AAD (Beckman Coulter). Cell sorting was performed on a FACSVantage or a FACSAria (BD Biosciences) cell sorter. Data were reanalyzed with FlowJo X (Treestar) software. AGM CD144⁺CD45⁻ endothelial cells, CD34^lowKit^lowCD45^pos differentiated hematopoietic cells and CD34^highKit^high HSPCs were isolated by flow cytometry from E11.5 C57Bl/6 mouse embryos. Three and five independent experiments were performed for the hematopoietic fractions (75-82 embryos per experience) and endothelial fraction (29-80 embryos per experience), respectively.

Embryos were fixed overnight in 4% paraformaldehyde (PFA; Sigma-Aldrich) at 4°C, transferred in a solution containing PBS/sucrose 15% for 48 h, impregnated in PBS/15% sucrose/7.5% gelatin for 2 h at 37°C, embedded at room temperature, frozen in isopentane at −65°C using liquid nitrogen and stored at −80°C until used for cryosectioning. Frozen sections (20 µm) of E11.5 embryos were cut on a cryostat (Leica). Antibodies used for staining the sections were: monoclonal mouse anti-mouse/human Tnc (1/50, Thermo Fisher Scientific, MA5-16086), polyclonal goat anti-mouse Dcn (1/200, R&D Systems, AF1060), monoclonal mouse anti-human Dcn (1/200, R&D Systems, MAB143), polyclonal goat anti-mouse Ntrk3 (1/200, R&D Systems, AF1404), polyclonal rabbit anti-mouse CD31 (1/200, Abcam, Ab28364), monoclonal mouse anti-mouse α-smooth muscle actin (1/500, Sigma-Aldrich, A5228), polyclonal rabbit anti-mouse Iba1 (1/1000, Wako, SAB5701363), monoclonal mouse anti-mouse Tuj1 antibody (1/1000, BioLegend, MMS-435P) and polyclonal goat anti-human CD105 (1/100, R&D Systems, AF1097). These primary antibodies were visualized with bovine anti-goat IgG-Alexa647 (1/400, Jackson ImmunoResearch, 805-605-180), goat anti-mouse IgG2a-Alexa488 (1/200, Invitrogen, A-21131), donkey anti-mouse IgG-Alexa488 (1/400, Jackson ImmunoResearch, 715-545-150) and donkey anti-rabbit Cy3 (1/400, Jackson ImmunoResearch, 711-165-152). Fluorescence imaging was performed using a fluorescence Leica microscope and confocal Leica TCS SP5 at the imaging facility, Institut de Biologie Paris-Seine, Sorbonne Université.

Embryos were incubated at room temperature (RT) for 3 h in a solution (PBSGT) containing PBS 1×0.2% gelatin, 0.5% Triton X-100 (Sigma-Aldrich) and 0.01% thimerosal (Sigma-Aldrich). Embryos were then incubated with the first antibodies as described above in PBSGT at RT for 2-3 days. After six washes, samples were incubated with secondary antibodies in PBSGT overnight, washed and processed for 3DISCO tissue clearing as described by Belle et al. (Belle et al., 2014). Images were acquired using an ultramicroscope (LaVision BioTec). 3D imaging, stack image conversion and movies were generated using Vision 4D 3.0 software (Arivis).

The mouse embryo sections were immersed in the fixative solution 1.5% glutaraldehyde and 1% PFA in 0.1 M cacodylate (pH 7.4) for 1 h at room temperature and overnight at 4°C. Samples were washed with cacodylate buffer 0.1 M, incubated for 1 h in 0.2% OTE in 0.1 M cacodylate buffer and washed again with cacodylate buffer 0.1 M. To keep them flat, the sections were embedded between two slides with 4% low melting point agarose in 0.1 M cacodylate. Afterwards, the samples were incubated 1 h in 1% osmium tetroxide in 0.1 M cacodylate buffer and washed with deionized water. The mouse embryo sections were dehydrated through graded concentration of ethanol (50, 70, 95, 100%). Samples were pre-embedded with graded concentration of anhydrous acetone and EPON resin mix (3:1, 1:1, 1:3) and embedded with 100% EPON resin (polymerized at 60°C for 24 h) between two glass slides coated with silicone. Finally, the flat sections were mounted on EPON resin blocks (preliminarily polymerized for 72 h) and cut with the Ultracut ultramicrotome (UCT, Leica) at 70 nm thickness. Ultra-fine sections deposited on 75-mesh copper grids were contrasted with 2.5% uranyl acetate and 2% lead citrate and coated with a 2 nm carbon film (ACE600, Leica). The grids were observed at 20 kV, with a STEM detector in a high-resolution Field-Emission SEM (GeminiSEM500, Zeiss) with a 20 µm aperture diameter. Images were acquired with a 1024×768 definition, with a pixel dwell time of 1.6 μs and a line averaging of 15.

The following zebrafish lines were used: Tg(kdrl:mCherry) (Chi et al., 2008), Tg(Gata2b:Gal4/UAS:LA-GFP) (Butko et al., 2015), Tg(flt1:dtomato) (Bussmann et al., 2010), Tg(flt4:YFP) (Hogan et al., 2009), Tg(cmyb:GFP) (North et al., 2007) and dcn^sa20247 (from the European Zebrafish Resource Center, 23113).

Zebrafish were maintained according to the University of California, San Diego, International Animal Care and User Committee and aquatic animal guidelines of Sorbonne Université. Wild-type or transgenic zygotes were injected with an MO solution and incubated at 28.5°C until they reached the stage of interest. Embryos were collected, staged, fixed and processed for in situ hybridization as previously described (Thisse et al., 1993). Antisense dcn MOs (Gene Tools) were used to target/block the 5′-untranslated region/translation start (MO^START) or a splice junction (MO^SPLICE). Dcn MO sequences were: MO^START GACAGGCCGATTTCATGTTGCTGAC and MO^SPLICE GCAGACCTGGGCATTTTGACACAGA (Zoeller et al., 2009). The Gdf6a MO sequence (Gene Tools) was: GCAATACAAACCTTTTCCCTTGTCC (Sidi et al., 2003). To determine the phenotype of dcn morphants and mutants on hematopoietic production, embryos were in embedded in 80% glycerol and photographed using a Leica IC80 HD camera mounted on a Leica M60 stereoscope. For quantification of runx1- and cmyb-expressing cells in the dorsal aorta, all morphant, mutant and control embryos were considered and quantified individually.

For flow cytometry, fluorescent embryos (∼300 embryos) were tail cut at 48 hpf, and the tissue was treated with Liberase (Roche) and subsequently manually dissociated. The cell suspension was filtered through a 30 µm mesh filter (CellTrics) and resuspended in 1× PBS+1 mM EDTA with SYTOX Red. Cells were sorted using FACSAria II 6 lasers. Data was processed using FlowJo software and statistics were calculated using GraphPad Prism 10. RT-qPCR for sorted cells was performed using TaqMan gapdh (#Dr03436842_m1) and dcn (#Dr03093424_m1) primers. For RT-qPCR to investigate Tgfβ ligand and receptor expression, RNA (∼25 embryos per condition) was isolated using the Qiagen RNeasy Micro Kit from zebrafish trunks at 30 hpf, and 200 ng of RNA was used to generate cDNA using the Invitrogen SuperScript III. RT-qPCR was performed with SYBR^TM reagents. Samples were run in triplicate, with four biological replicates using the primers listed in Table S4. Statistics were calculated using GraphPad Prism 10.

Statistical analysis, sample sizes and P-values are described in figure legends and in the appropriate Materials and Methods sections for computational analysis. For the loss-of-function experiments, statistical analysis was performed using GraphPad Prism software and two-tailed unpaired t-tests and Mann–Whitney nonparametric tests were used to evaluate the statistical significance.

We thank Claire Francastel (CNRS UMR7622), Laurent Yvernogeau (CNRS UMR7622) and Pascal Legendre for critical reading of the manuscript, and all laboratory members for helpful discussion. We thank Srinath Ramkumar and Didier Stainier (Max Planck Institute for Heart and Lung Research, Germany) for providing us with the zebrafish dcn construct. We also thank Monara Angelim and Pascal Legendre (CNRS UMR8246) for providing us with the anti-Tuj1 and -Iba1 antibodies, members of the imaging (Institut de la Vision and Institut de Biologie Paris-Seine) facilities for their expertise and advice, Michael Trichet and Alexis Canette (Service de microscopie électronique, Institut de Biologie Paris Seine, Sorbonne Université, CNRS) for their support and advice on the transmission electron microscopy experiments, and Sophie Gournet (CNRS UMR7622) for drawing assistance.