Cell lineage specification and signalling pathway use during development of the lateral plate mesoderm and forelimb mesenchyme

The lateral plate mesoderm (LPM) is a transient tissue that forms from the mesodermal germ layer during vertebrate embryogenesis. The LPM produces a remarkable diversity of cell and organ types during development, including the smooth muscle and connective tissues of the cardiovascular, respiratory and digestive systems, and the endoskeleton of the limbs (Nishimoto and Logan, 2016; Prummel et al., 2019, 2020; Selleck and Stern, 1991). LPM specification and limb morphogenesis are complex processes underpinned by coordinated signalling events and gene expression dynamics across development (Loh et al., 2016). The primitive mesoderm elongates along the anteroposterior (A-P) body axis, and undergoes mediolateral specification to form axial (notochord), paraxial (somite), intermediate and lateral plate mesoderm domains, through specific combinations of BMP, FGF and WNT signals (Loh et al., 2016; Tonegawa et al., 1997). The LPM forms as bilateral sheets, comprising anterior (aLPM) and posterior (pLPM) domains, which are spatially defined by FGF8 and retinoic acid signalling boundaries (Cunningham et al., 2013; Zhao et al., 2009). LPM identity is further established and maintained by localized BMP4 signalling and WNT antagonism (Tonegawa and Takahashi, 1998; Tonegawa et al., 1997; Yoshino et al., 2016).

Cells of the undifferentiated LPM undergo further dorsoventral subdivision to form the somatic and splanchnic LPM layers, separated by the embryonic coelom. The somatic LPM gives rise to tissues that form the body wall, amnion and limbs, while the splanchnic LPM forms the viscera and connective tissue lining the digestive system. Differentiation of these diverse tissues is accompanied by activation of key transcription factors during fate specification (Agarwal et al., 2003; Firulli et al., 1998; Harvey et al., 2002; Mahlapuu et al., 2001; Martin et al., 1995; Rallis et al., 2003). However, despite this understanding of LPM differentiation and tissue diversification, little is known about the mechanisms that drive LPM subdivision and specification towards a somatic and splanchnic LPM fate (reviewed by Prummel et al., 2020).

Initial specification of the LPM is accompanied by activation of transcription factors FOXF1, HAND1, OSR1 and PRRX1 (Kuratani et al., 1994; Loh et al., 2016; Martin et al., 1995; Ormestad et al., 2004; Peterson et al., 1997). LPM subdivision is thought to occur through BMP signalling from the overlying ectoderm, through activation of somatic LPM genes (Funayama et al., 1999). In the presence of ectodermal BMPs, PRRX1 becomes restricted to and IRX3 is activated in the somatic LPM (Funayama et al., 1999; Mahlapuu et al., 2001; Ocaña et al., 2012). OSR1 becomes restricted to the intermediate mesoderm, while HAND1 and FOXF1 are maintained in the splanchnic LPM (Funayama et al., 1999; Mahlapuu et al., 2001). However, the mechanisms that drive somatic LPM development remain unresolved. Prrx1, Irx3 and Osr1 are dispensable for somatic LPM formation, as mouse null mutants for each of these genes are viable and do not possess aberrant LPM phenotypes (Li et al., 2014; Martin et al., 1995; Wang et al., 2005). Conversely, Foxf1 null mutants are embryonic lethal with partial to incomplete subdivision of the LPM, misexpression of somatic LPM genes in splanchnic LPM and gut patterning defects (Firulli et al., 1998; Mahlapuu et al., 2001). This suggests that FOXF1 plays a role during LPM subdivision, acting to antagonize expression of somatic LPM genes that are activated by ectodermal BMP signals. However, the mechanisms responsible for establishing somatic LPM identify remain undefined and are essential to our understanding of LPM biology and early embryogenesis.

Once formed, the somatic LPM is specified into limb-forming domains through retinoic acid (RA) signalling and HOX gene expression boundaries (Nishimoto et al., 2015; Tanaka, 2016). RA, HOX and β-catenin/TCF/LEF are thought to initiate limb development through activation of the T-box transcription factor TBX5 (Agarwal et al., 2003; Logan et al., 1998; Minguillon et al., 2012; Nishimoto et al., 2015; Rallis et al., 2003). Limb outgrowth begins with TBX5-induced activation of FGF10, which together are proposed to initiate a localized epithelial-to-mesenchymal transition (EMT) of the somatic LPM (Gros and Tabin, 2014). FGF10 activates FGF8 in the overlying ectoderm, establishes a positive-feedback loop, and drives proliferation of the limb bud mesenchyme and limb outgrowth (Moon and Capecchi, 2000; Nishimoto et al., 2015; Ohuchi et al., 1997). Limb patterning sees further activation of networks of transcriptional regulators and morphogens to promote outgrowth maintenance and patterning, which are well characterized (for comprehensive reviews, see Tickle, 2015; Zuniga, 2015). Yet, despite these known events, those preceding the TBX5-dependant limb regulatory pathway in the somatic LPM are less well understood.

In this study, we resolve the ambiguity underlying early LPM differentiation, subdivision and limb initiation in the developing chicken forelimb field, using single-cell transcriptomics. We construct a high-resolution gene expression atlas across the developing limb field, corroborating known and identifying previously uncharacterized tissue-specific marker genes, signalling pathways and ligand-receptor pairs that are active during LPM differentiation and early limb initiation. We address a long-standing gap in the literature, defining the transcriptional networks underlying subdivision and specification of the somatic and splanchnic LPM. Furthermore, we identify the sequential activation of transcription factors during limb initiation, including early activation of TWIST1 and other genes, upstream of the TBX5 and the FGF10-FGF8 feedback loop. Interrogation of TWIST1 expression and localization in the early somatic LPM reveals a putative role during EMT of the somatic LPM. Finally, we reveal that BMP signalling is not sufficient for LPM subdivision, as previously suggested (Funayama et al., 1999), but is necessary for determining somatic LPM fate and limb initiation through activation of TBX5 and FGF10. Together, these findings provide a robust overview of the developmental landscape underlying formation and specification of the LPM, representing one of the earliest events of vertebrate organogenesis.

We performed single-cell RNA sequencing of embryonic chicken forelimb fields and reconstructed the underlying cell populations. Briefly, tissues corresponding to the presumptive forelimb field, i.e. lateral to somites 20-25, were dissected from embryonic day (E) 1.5, 2.5 and 3.5 embryos, corresponding to approximate Hamburger-Hamilton stages (HH) 10, 14 and 18, respectively (Hamburger and Hamilton, 1951) (Fig. 1A). These stages cover the emergence of the early LPM, subdivision and specification of the LPM into the somatic and splanchnic layers, and limb initiation and early outgrowth (Newton and Smith, 2020). Dissected tissues were dissociated into single cells, FACS sorted to remove dead and dying cells, then processed through the 10x Chromium system. After quality filtering, a total of 15,355 cells, corresponding to 5273 cells from E1.5, 6856 from E2.5 and 3226 from E3.5 embryos, were recovered, with expression profiles for 16,779 genes. To remove low read count or low diversity cells, we applied an additional strict filtering threshold of 2000 UMI counts per cell, yielding a total of 3262 cells, with 1210 from E1.5, 1313 from E2.5 and 739 from E3.5 embryos. Cell transcriptomic relationships were visualized with global t-distributed Stochastic Neighbour Embedding (tSNE) dimension reduction, which showed a distinct separation of cell types according to germ layer, cell cycle phase and embryonic stage (Fig. 1B,C).

Unsupervised clustering of the chicken E1.5, E2.5 and E3.5 presumptive forelimb cell populations revealed 13 transcriptionally distinct cell clusters (c) that represented embryonic vasculature and tissues derived from the ectoderm, mesoderm and endoderm germ layers, broadly covering known cell types within the presumptive forelimb field (Fig. 1C,D). Tissue and cell type identities were assigned to each cluster based on their differential gene expression profiles (observed as differences in log-fold change between clusters; Fig. 1E,F) and corresponding spatiotemporal expression patterns observed throughout the developing chicken embryo on the GEISHA chicken gene expression database (Bell et al., 2004; Darnell et al., 2007). Cluster-specific gene expression is shown in Fig. 1E,F and Table S1. The embryonic ectoderm comprised two clusters (c2 and c11) defined by expression of FABP3 and WNT6. However, c11 showed expression of FGF8, revealing these cells as progenitors that contribute to formation of the apical ectodermal ridge (AER). The embryonic endoderm (c7) was defined by unique expression of GUCA2B, TTR and SHH (not observed in limb bud mesenchyme due to the early stages sampled). Embryonic vasculature (c8) showed unique expression of CDH5, and red blood cells (c13) expressed haemoglobin subunit HBBA.

The embryonic mesoderm was found to contribute the largest overall number of cells and was largely defined by expression of the established marker PRRX1. The mesoderm was composed of six clusters (c1, c3, c4, c5, c10 and c12), which represented known mesodermal tissue types during development, separated by embryonic stage (Fig. 1B,C). The earliest identified cell type was that of the primitive mesoderm/LPM (c4), which was largely made up of E1.5 cells with expression of primitive streak and early mesodermal markers MSGN1, EVX1 and CDX4 (Alev et al., 2010). Cells of the E2.5 mesoderm largely represented more distinct tissues, including the somatic LPM (c3) – high PRRX1 and low TBX5; splanchnic LPM (c5) – COLEC10; and paraxial mesoderm (c1) – TCF15 (Fig. 1D,E). Finally, E3.5 cell clusters represented more differentiated tissues, such as the extra-embryonic LPM/amnion (c12) – AQP1; and limb bud mesenchyme (c10) – high TBX5 and FGF10. We also detected two clusters (c6 and c9) that possessed ubiquitous expression of ectodermal and mesodermal markers, but also high ribosomal and mitochondrial counts. These were defined as low diversity cells that were excluded from subsequent analyses. Together, these clusters represent all of the known major tissues types within the developing forelimb field that, for the first time, now possess detailed transcriptional profiles (Table S1).

Specification of the mesoderm during embryogenesis is influenced by dynamic intrinsic and extrinsic signalling between the surrounding germ layers and developing tissues. In particular, members of the BMP, FGF, HH and WNT signalling pathways are known to communicate between the germ layers to induce differentiation and development (Loh et al., 2016). However, the precise ligands and receptors that facilitate different aspects of mesoderm development are unclear. We therefore examined global signalling patterns and ligand-receptor crosstalk between the global cell types during specification of the forelimb field from the LPM using CellChat (Jin et al., 2021). This analysis revealed extensive signalling pathway use between different cell types and tissues (Fig. 2A), which dynamically changed during lineage allocation of key tissue types. The E1.5 primitive MES/LPM (c4) showed the highest level of signalling pathway activity among the tissues studied, with active signalling through more than half of the predicted pathways. This included signalling through the important non-canonical WNT, FGF, HH and BMP pathways, as well as presumptive signalling through midkine (MK) and pleiotrophin (PTN), ephrin B, semaphorin (SEMA3 and SEMA6), chemokine ligand (CXCL) and adhesion factors laminin, JAM and nectin.

Pathway use was observed to significantly change during subdivision of the LPM into somatic and splanchnic LPM by E2.5, including decreases in ncWNT, FGF, CXCL and semaphorin signalling. Splanchnic LPM (c5) formation involved HH, AGRN, FN1, CD99 and ephrin signalling, and decreased BMP signalling. In contrast, formation of the somatic LPM (c3) was characterised by decreased HH, maintained BMP, and activated collagen and nectin signalling. Finally, development of the limb bud mesenchyme (c10) by E3.5 was characterised by specific activation of WNT and FGF signalling, confirming known interactions, as well as strong activation of ephrin A and ANGPTL. Although these observations indicate that differentiation of multipotent LPM progenitors involves diverse signalling crosstalk (Fig. 2A), we focused on TGFβ, WNT, BMP, FGF and HH, given their known roles in influencing LPM specification, subdivision and limb development (Loh et al., 2016).

To begin to define these complex signalling relationships, we interrogated TGFβ, WNT, BMP, FGF and HH signalling pathways active within sender and receiver cell types, allowing construction of dynamic signalling networks (Fig. 2B). TGFβ signalling was responsible for development of the embryonic vasculature, with ligands only being received by blood vessels (c8). Early LPM development was observed to occur through BMP and HH signalling. BMPs were sent by the LPM (c4) and ectoderm, and received by the early LPM (c4) and somatic LPM and/or limb bud (c10), recapitulating known local and ectodermal signalling (Funayama et al., 1999; Tonegawa et al., 1997). Interestingly, the endoderm was identified as a significant sender of HH signals, being received by the early LPM, splanchnic LPM and paraxial mesoderm, affirming its known role during gut formation from splanchnic LPM (Roberts et al., 1995; Yoshino et al., 2016). During later limb initiation, the somatic LPM (c3) and limb mesenchyme (c10) were strong senders of WNT signals, received by the E3.5 limb bud mesenchyme (c10) and amnion (c12) (Takeuchi et al., 2003). The limb mesenchyme (c10) was further identified as a sender of FGF signalling, received by the AER ectoderm (c11), confirming the known role of secreted FGF10 signalling during limb initiation (Ohuchi et al., 1997). However, we did not identify the AER ectoderm (c11) as a sender of FGF signals, nor FGF8 expression in the c11 AER cluster (Table S1), suggesting our sampling of AER was prior to establishment of the FGF10-FGF8 feedback loop.

To further define the influence of these signalling pathways, we next examined the underlying BMP, FGF, HH and WNT ligand-receptor pairs facilitating these dynamic events. This revealed disparate patterns of ligand and receptor use among the different cell populations, which are detailed in Fig. 2C. During LPM subdivision, endoderm-mesoderm HH signalling occurred exclusively through SHH activation of PTCH1. Ectoderm-mesoderm BMP signalling is achieved through BMP2 and BMP7, as previously described (Funayama et al., 1999), in combination with ACVR1 and BMPR2, ACVR2A or ACVR2B receptors as heterodimeric complexes (Tajer et al., 2021). Interestingly, we also observed a strong source of local BMP5 signalling from the somatic LPM, although its role in LPM development is unknown. During later limb development, there were diverse patterns of WNT ligand signalling in tissues throughout the limb field, including ectodermal WNT3A, WNT4, WNT6, WNT7A and WNT10A, but restricted expression of FRZD and LRP6 receptors only within limb bud mesenchyme and amnion (extra-embryonic LPM). Together, these data reveal that signalling throughout the limb field is achieved through combinations of both ubiquitous and tissue-specific ligand-receptor expression patterns, which are summarized in Fig. 2D. Overall, the molecular signalling events underlying LPM development and subdivision appeared to be dependent on BMP and HH signalling, while limb development is largely driven through FGFs and WNTs. However, it is important to note that limb development is also dependent on RA signalling (Nishimoto et al., 2015; Stratford et al., 1996), although we were unable to predict these RA ligand-receptor interactions due to the complex biosynthesis pathway.

Next, we sought to determine the mechanisms underlying fate choice of the early LPM during differentiation, subdivision and limb initiation. This was achieved by first sub-setting the dataset to cells of mesodermal origin. Mesodermal cells were reprocessed with UMAP dimension reduction and Leiden clustering, allocating the six previous mesodermal clusters (c1, c3, c4, c5, c10 and c12; Fig. 1A) into 12 new sub-clusters (mc1-12; Fig. 3A). These captured the previously determined cell types, as well additional intermediate, transitional and terminal cell types during LPM development (Fig. 3B). Namely, these clusters labelled cells from the somitic (mc9) and intermediate mesoderm (mc10), primitive mesoderm (mc2), early (undifferentiated) LPM (mc1 and mc7), somatic (mc3) and splanchnic (mc4) LPM, visceral and/or gut precursors (mc11), limb bud mesenchyme (mc5), non-limb flank LPM (mc12), extra-embryonic LPM (mc8) and amnion (mc6; Fig. 3A). Diagnostic cluster gene markers and log-fold changes are listed in Table S2.

With defined cell and tissue types, differentiation pathways throughout the LPM were first explored using estimates of RNA velocity. Cell velocity estimates revealed specific, directional transcriptional trajectories between cells as they transitioned from undifferentiated E1.5 (∼HH10) precursors towards their distinct tissue fates, connecting three out of the four E3.5 (∼HH18) cell populations (Fig. 3B,C). Interestingly, we identified a heterogeneous population of early undifferentiated LPM cells with low directional velocity compared with other neighbouring clusters, despite existing in various stages of the cell cycle (Fig. 3C). This suggests that the multipotent cells of the early LPM exist in a transiently uncommitted state, before rapidly committing towards defined cell fates by E2.5 (∼HH14), likely as a response to secreted ectodermal/endodermal signals (Fig. 2D). Furthermore, the somitic and intermediate cell clusters of the paraxial mesoderm, despite forming from the primitive mesoderm, did not show a continuum of directional RNA velocities from these precursors. As we did not intentionally sample paraxial tissues, it is unclear whether there were not enough cells isolated to represent the complete differentiation trajectory, or whether these cells possess an earlier embryonic origin to the primitive mesoderm cells captured in our data. As such, we chose to not include this lineage in subsequent analyses.

To further define the pathways of differentiation that arise throughout the LPM, we explored these through LPM lineage generation. The E1.5 primitive mesoderm (mc2) represented the earliest identified cell type and therefore was determined as the root node of mesoderm differentiation. A principal neighbour graph was fit with monocle3 (Trapnell et al., 2014), revealing four major lineages (L1-L4; Fig. 3D), which supported RNA velocity estimates (Fig. 3C). These lineages describe the transition from E1.5 (∼HH10) primitive mesoderm cells to E3.5 (∼HH18) visceral and/or gut precursors (mc11; L1), non-limb somatic LPM (mc6 and mc12; L2), limb bud mesenchyme (mc5; L3) and somitic/intermediate mesoderm (mc9-10; L4). Importantly, a distinct bifurcation point was observed within the early undifferentiated LPM (mc1/mc7), marking LPM subdivision and lineage specification into the somatic (mc3) and splanchnic (mc4) LPM tissue layers (Fig. 3D).

The emerging somatic and splanchnic LPM layers are accompanied by specific expression of IRX3 and FOXF1, respectively (Funayama et al., 1999; Mahlapuu et al., 2001). To confirm the accuracy of this lineage bifurcation point, IRX3 and FOXF1 expression were examined in mesodermal cells. Indeed, cells displayed mutually exclusive expression of IRX3 or FOXF1 in complementary domains following the LPM bifurcation point (Fig. 3D) with FOXF1+ cells observed in the splanchnic mesoderm and viscera precursors (mc4/7/11), and IRX3+ cells in the somatopleure, flank LPM and limb bud mesenchyme (mc1/3/5) (Fig. 3D). These data confirm the accuracy of somatic and splanchnic LPM bifurcation. However, also included within this lineage trajectory was an irregular branch linking the splanchnic LPM to the extra-embryonic LPM and amnion (Fig. 3D,E, dashed line), which does not accurately represent its origins in vivo. Although we were able to resolve this branch through additional k-mean principal graphs, these each came at the expense of the somatic-splanchnic LPM bifurcation point, with neither correct topology present in a single graph without losing cluster resolution (Fig. S1). Therefore, we focused solely on the somatic-splanchnic LPM branching, as confirmed by FOXF1 and IRX3 expression, for subsequent lineage trajectory analysis.

Using this defined lineage of LPM differentiation, key gene expression dynamics underlying subdivision of the somatic and splanchnic LPM, and specification towards a limb or viscera fate were examined. This was achieved through calculations of pseudotime along the viscera and limb lineages (Fig. 3E) (Trapnell et al., 2014), where genes with significant expression changes along each branch were identified through Moran's I spatial autocorrelation tests (Tables S3 and S4). The top differentially regulated genes along a given lineage were visualized through expression heatmaps, producing gene expression modules that were repressed or activated across pseudotime. This revealed a pseudo-temporal hierarchy of dynamically expressed genes during LPM subdivision, including marker genes that were specifically activated immediately after LPM bifurcation (Fig. 3F,G; Figs S2,S3). Each lineage shared co-expression of primitive mesoderm and/or early LPM gene modules, which were repressed before subdivision. Immediately afterwards, there was activation of shared or lineage-specific genes and gene modules (Fig. 3F,G). During splanchnic LPM differentiation and mesodermal viscera development, nine modules of genes were co-activated, including, but not limited to, activation of transcription factors GATA6, NKX2-3, TCF21 and HAND1, and repression of PRRX1 (Fig. 3F; Fig. S3). Somatic LPM differentiation featured co-activation of six transcriptionally distinct modules, including, but not limited to, increasing PRRX1 expression, proximal activation of OLFML3, NR2F2 and the principal retinoic acid synthesis gene ALDH1A2 (Raldh2), followed by transcription factors IRX3, IRX6 and TWIST1 (Fig. 3G; Fig. S3). Immediately afterwards, genes modules involved in limb initiation were activated, beginning with expression of TBX5 and TBX2, and followed by limb bud mesenchyme genes LMX1B, PRRX2, WIF1 and FGF10. Interestingly, although TWIST1 is known to be expressed in the somatic LPM and limbs (Tavares et al., 2001), its rapid activation prior to the limb initiation module indicates a previously uncharacterised role during early limb development.

Limb outgrowth is first initiated by a localized EMT of the somatic LPM, proposed to be driven by TBX5 and FGF10 (Gros and Tabin, 2014). However, before TBX5 and FGF10 expression in the somatic LPM, we observed strong expression of PRRX1 and activation of TWIST1 (Fig. 3G), which are known EMT inducers during development and cancer (Fazilaty et al., 2019; Ocaña et al., 2012). We therefore further investigated PRRX1, TWIST1 and TBX5 expression in pseudotime and by in situ hybridization, confirming their spatiotemporal expression profiles during limb development. PRRX1 was expressed in the forelimb field and LPM as early as E1.5/HH10 before subdivision, before becoming regionalized and strongly expressed in the somatic LPM and limb bud (Fig. 4A). In comparison, TWIST1 was specifically activated within in the somatic LPM at ∼E2.0/HH12-13, immediately before TBX5, which was activated in the somatic LPM by ∼E2.5/HH14 (Fig. 4A).

Additional immunofluorescent labelling in HH12 chicken forelimb fields revealed strong, specific expression of TWIST1 protein in somatic LPM cells before limb bud EMT (Fig. 4B). At this stage, the somatic LPM existed as a pseudostratified columnar monolayer with meso-epithelial characteristics, expressing the transmembrane protein N-cadherin, but not the typical epithelial marker E-cadherin, which was restricted to the ectoderm and endoderm. TWIST1 was specifically detected in cells of the somatic LPM and somite, with some TWIST1⁺ cells appearing to delaminate from the somatic LPM (Fig. 4B, arrows). By E2.5/HH14 the somatic LPM was no longer a distinct monolayer, indicated by increased numbers of TWIST1⁺ cells and reduced N-cadherin, which became localized to the apical edge of the somatic LPM. By ∼E3.0/HH16, the early forelimb bud was distinct, populated by large numbers of TWIST1⁺ mesenchymal cells (Fig. 4B), which were also present in the sclerotome. We also examined whether other EMT transcription factors were co-expressed with TWIST1 in the somatic LPM during EMT but did not detect enrichment of other EMT transcription factors in somatic LPM or limb bud mesenchyme clusters (Fig. 4C), or along the trajectory (Fig. 3C; Fig. S3). As such, the early co-expression of PRRX1 and TWIST1 in the somatic LPM (Fig. 4A,B) suggests a putative role in somatic LPM EMT, before TBX5-induced limb initiation and FGF10-dependant outgrowth.

The BMP family of signalling factors, enriched in our receptor-ligand dataset (Fig. 2), have previously been implicated in LPM subdivision and early somatic LPM cell fate (Funayama et al., 1999). Finally, we sought to define their influence on LPM specification, subdivision, somatic LPM fate and early limb initiation through targeted BMP inhibition between the early ectoderm and developing LPM. To do this, we performed electroporation of the secreted BMP antagonist noggin (CAG-GFP-T2A-NOG) into the ∼HH9 forelimb field ectoderm to inhibit BMP signalling to the underlying LPM (Fig. 5A), during the window where it is most active (Fig. 2). Control embryos electroporated with GFP (CAG-GFP-T2A) reproducibly developed limb buds, whereas embryos electroporated with noggin failed to form forelimb buds (Fig. 5A) and, in instances where the expression domain expanded posteriorly, failed to form hindlimb buds (Fig. S4). Interestingly, inhibition of ectodermal-LPM BMP signalling was not sufficient to disrupt LPM subdivision, as previously implied (Funayama et al., 1999). Rather, the somatic and splanchnic LPM still formed and were separated by the coelom, observed by robust FOXF1 expression in the splanchnic LPM, albeit with a reduced limb bud (Fig. 5A,B). However, ectodermal BMP signal inhibition greatly stunted somatic LPM formation, observed by noticeable reduction of PRRX1 expression in the limb field and somatic LPM tissue sections, compared with controls (Fig. 5B).

Inhibition of ectodermal BMP signalling via noggin resulted in reduced outgrowth of the limb bud mesenchyme and a failure in formation of the developing forelimb bud (Fig. 5C). Strikingly, this was observed to occur in response to attenuated expression of TBX5 and to failed activation of both FGF10 in somatic LPM and FGF8 in the AER (Fig. 5C; Fig. S4B), suggesting ectodermal BMP signalling regulates the TBX5-FGF10-FGF8 limb initiation mechanism. Conversely, we found that BMP inhibition did not influence expression of TWIST1, which showed similar expression and protein localization in forelimb sections compared with controls, despite reduced overall proportions of limb bud mesenchymal cells (Fig. 5A,C). Thus, TWIST1 activation appears to be dependent on other signalling inputs that could not be determined, which may include RA signalling (Tavares et al., 2001). Together, our data confirm that ectodermal-mesodermal BMP signalling crosstalk is not required for LPM subdivision, as previously suggested (Funayama et al., 1999), but is necessary for somatic LPM specification and early limb initiation through activation of the TBX5, FGF10 and FGF8 feedback loop to regulate forelimb outgrowth.

Single-cell transcriptomics have provided powerful methods with which to define cell lineage trajectories underlying multiple aspects of mesoderm and limb development (Feregrino et al., 2019; Gerber et al., 2018; Han et al., 2020; Johnson et al., 2020; Loh et al., 2016; Mahadevaiah et al., 2020; Pijuan-Sala et al., 2019; Scialdone et al., 2016). However, the initial developmental events underlying specification and diversification of the LPM were previously undefined. We add to this body of literature, revealing key mechanisms underlying specification of the early LPM towards a limb or viscera cell fate. LPM specification originates in primitive mesodermal precursors, which follow four distinct differentiation pathways from E1.5 (∼HH10) to E3.5 (∼HH18) (Fig. 3). Initially, LPM progenitor cells (FOXF1+) displayed a large degree of cellular heterogeneity (Fig. 3A-D), which was accompanied by extensive signalling pathway use (Fig. 2). This revealed that the emerging LPM forms in a transient, multipotent state, where cells are transcriptionally primed to respond to rapid changes in the extrinsic signalling environment to initiate LPM lineage specification. This was observed where the LPM heterogeneity was rapidly resolved by E2.5, with cells committing to defined differentiation lineages, including LPM subdivision, in response to dynamic changes in signalling pathway activation and repression.

Early establishment and differentiation of the LPM was predicted to occur through signalling interactions within the BMP and HH pathways, originating from the ectoderm and endoderm, respectively, supporting previous observations (Funayama et al., 1999; Roberts et al., 1995; Yoshino et al., 2016). Although our analyses were unable to directly predict the transcriptional targets of these signalling interactions during LPM development, it is possible that HH is the inducing signal for early LPM specification. This hypothesis is supported where a SHH-FOXF1-BMP4 signalling axis drives ureter development (Bohnenpoll et al., 2017), and SHH activates BMP4 during gut development (Roberts et al., 1995). Given that we observe SHH signalling from the endoderm to PTCH receptors in the early LPM (Fig. 2) (Yoshino et al., 2016), and that FOXF1 and BMP4 are markers that specify the early LPM (Fig. 3D) (Mahlapuu et al., 2001; Tonegawa et al., 1997; Winnier et al., 1995), it is plausible that the SHH-FOXF1-BMP4 axis participates in initial LPM patterning and in establishing its identity.

The mechanism underlying LPM subdivision have been previously suggested to occur through BMP2/7 signals derived from the ectoderm (Funayama et al., 1999). As our global signalling analysis recovered these BMP signalling events during early LPM development (Fig. 2), we investigated whether these are indeed responsible for formation of the somatic and splanchnic LPM. Targeted inhibition of ectodermal-mesodermal BMP signalling, using the secreted BMP antagonist noggin, was not sufficient to drive LPM subdivision, observed by formation of both layers separated by the coelom (Fig. 5). This reveals that other undefined mechanisms are responsible for LPM subdivision, but may be achieved through endodermal HH signalling through activation and maintenance of FOXF1. Foxf1 null mice possess incomplete LPM subdivision, somatic-splanchnic LPM fusion and misexpression of somatic LPM genes (i.e. Irx3) (Firulli et al., 1998; Mahlapuu et al., 2001). Thus, it is possible that a SHH-induced FOXF1-dependent mechanism initiates LPM subdivision, before ectodermal BMPs induce a somatic LPM identity on the newly formed dorsal layer.

Previously, PRRX1 and IRX3 were the only defining markers of the early somatic LPM (Funayama et al., 1999), yet are both dispensable for somatic LPM development (Li et al., 2014; Martin et al., 1995). Through our analyses, we have defined additional suites of genes and regulators that accompany LPM subdivision into somatic and splanchnic LPM fates (Fig. 3F,G; Figs S2,S3). LPM subdivision occurs through initial repression of a shared module of early mesoderm genes, before seeing linage-specific activation of somatic and splanchnic gene modules (Fig. 3F,G). Of note, the first stage of LPM subdivision saw lineage-specific maintenance and activation of basic helix-loop-helix (bHLH) transcription factors HAND1 and TWIST1 within the splanchnic and somatic LPM lineages, respectively. HAND1 and TWIST1 possess known roles in later gut and limb development (Firulli et al., 1998; Krawchuk et al., 2010; Loebel et al., 2012; Wu and Howard, 2002), although our results suggest that they may play unappreciated roles during initial LPM fate specification. bHLH transcription factors form combinations of homo- and heterodimers with unique binding partners to dynamically regulate different biological processes (Fan et al., 2020). We observe expression of both ubiquitous (HAND2) and lineage-specific bHLH binding partners (TCF12 and TCF21) (Figs S2,S3) (Fan et al., 2020) throughout the LPM cell types. Thus lineage-specific activation of TWIST1 and HAND1 may mediate complementary, yet distinct, roles in target gene activation during the initial stages of LPM subdivision and differentiation.

Although PRRX1 and TWIST1 have known expression in the somatic LPM (Gitelman, 1997; Kuratani et al., 1994; Nohno et al., 1993; Tavares et al., 2001), our data reveal a potential, unreported role of these during somatic LPM development and EMT. Limb initiation begins with a localized EMT of the somatic LPM, suggested to depend on TBX5 and FGF10 (Gros and Tabin, 2014). However, PRRX1 and TWIST1 are known regulators of EMT in both development and cancer (Fazilaty et al., 2019; Ocaña et al., 2012; Qin et al., 2012), suggesting they may possess an uncharacterized role during early EMT or initial induction of the forelimb mesenchyme. As noted above, TWIST1 elicits different biological functions and activity thresholds through dimerization with other transcription factors (Fan et al., 2020; Krawchuk et al., 2010; Loebel et al., 2014), which notably includes PRRX1 in neural crest cells (Fan et al., 2021). Indeed, we identify robust, tissue-specific co-expression of TWIST1 (and PRRX1, not shown) in the meso-epithelial somatic LPM monolayer, before and after EMT and proliferation of the limb bud mesenchyme (Fig. 4). This putative, cooperative role is supported where exogenous prrx1a drives aberrant EMT and migration of the LPM in zebrafish, which can be rescued by co-knockdown of twist1b (Ocaña et al., 2012). Furthermore, Twist1 null mouse mutants possess severely atrophied forelimb buds, as a potential consequence of reduced TWIST1-induced EMT of LPM precursors (Chen and Behringer, 1995). Together, these data implicate PRRX1 and TWIST1 as early mediators of somatic LPM specification and EMT, although this requires further validation in vivo. Of note, it is unclear whether PRRX1 and/or TWIST1 may cooperate with TBX5 to facilitate EMT and induction of limb bud outgrowth.

The unexpected finding of this study was the severe inhibitory effect of BMP antagonism on early somatic LPM fate and limb specification and outgrowth. Inhibition of ectodermal BMP signalling, through ectopic expression of noggin, both reduced commitment and proportions of somatic LPM cells through inhibition of PRRX1, and severe atrophy of limb bud outgrowth through attenuation of TBX5, FGF10 and FGF8 (Fig. 5B). Activation of TBX5 in the somatic LPM has been previously suggested to depend on RA, with cooperative activation through HOX genes and β-CAT/TCF/LEF through the canonical WNT signalling pathway (Minguillon et al., 2012; Nishimoto et al., 2015). However, the WNT ligand for this activation has not been reported, and our data show no evidence of active WNT signalling in the early LPM/somatic LPM (Fig. 2). Supporting this, BMP-induced differentiation of iPSCs into LPM-like cells, defined by induction of TBX5, occurs in the absence of both WNT and RA (Loh et al., 2016). Instead, our data reveal a novel role for ectodermal BMP signalling in limb initiation through somatic LPM fate determination and TBX5 activation. Addition of BMP2 beads to chick hindlimb buds can induce ectopic TBX5 (Rodriguez-Esteban et al., 1999), supporting our finding that BMP induces activation of TBX5 in the forelimb. BMP antagonism is known to affect limb outgrowth by reducing AER formation via attenuated activation of FGF8 and MSX2 in the ectoderm (Capdevila and Johnson, 1998; Pizette and Niswander, 1999; Pizette et al., 2001). However, our data reveal that this is likely an indirect effect of attenuated TBX5 activation and reduced maintenance of the FGF10-FGF8 feedback loop. Together, these observations establish a new model of early forelimb development, where ectodermal BMPs regulate somatic LPM formation and limb initiation (Fig. 6) within the predefined HOX and RA-specified forelimb-forming domain.

A limitation of this study was our inability to predict the precise function of TWIST1 and the role of RA signalling during LPM specification and limb outgrowth (Cunningham et al., 2013; Gibert et al., 2006; Nishimoto et al., 2015; Zhao et al., 2009) due to its complex biosynthesis pathway. However, our somatic LPM lineage analysis revealed activation of the RA synthesis gene ALDH1A2 and RA-responsive transcription factor NR2F2 in the first module of somatic LPM differentiation, immediately followed by activation of TWIST1 (Fig. 3; Fig. S3). As TWIST1 was not responsive to BMP signalling (Fig. 5C), it is possible that TWIST1 is activated through a RA-mediated mechanism. In support of this, TWIST1 is induced in limb buds by ectopic RA-soaked beads (Tavares et al., 2001). Although additional data are required to determine the role of TWIST1 during somatic LPM differentiation and early limb initiation, these data suggest that somatic LPM fate is specified through the combined action of BMP and RA signalling (Fig. 6).

This study fills a long-standing gap in the literature, establishing the first transcriptional atlas of multipotent LPM progenitors and their differentiation through transitional and maturing cell types in the embryonic chicken forelimb field. We have defined the global signalling pathways and gene expression dynamics during early LPM development and shed light on the early cell-fate decisions preceding initiation of the vertebrate forelimb. Furthermore, although we were unable to identify the signal that initiates LPM subdivision, we defined a necessary role for ectodermal BMP signals in somatic LPM identity and limb bud formation. Although we have begun to characterize the gene regulatory networks underlying LPM specification before limb outgrowth (Fig. 6), more work is needed to define the role of key genes and signalling pathways during early specification of the LPM. Nevertheless, in this study we present the detailed overview of full cellular and developmental events from LPM specification to forelimb initiation, and provide exciting new avenues for scrutinizing the molecular events underlying fate specification, development and evolution of the vertebrate limbs.

Chicken eggs were collected at embryonic day (E) 1.5 (stage 10), E2.5 (HH14) or E3.5 (HH18). Embryos were separated from extra-embryonic membranes, rinsed in ice-cold DPBS and the LPM was then dissected. LPM tissues were digested with 0.05% Trypsin and EDTA, and incubated at 37°C for 15 min, with mechanical dissociation every 5 min until no clumps were visible. Enzymatic activity was stopped with addition of 10% FCS. The dissociated cells were spun at 400 g for 5 min, then resuspended in 1× EDTA and Propidium Iodide in DMEM (Gibco). Cells were filtered through a 70 μm Flowmi Cell Strainer (Scienceware), and viable cells were isolated by flow cytometry (Flowcore, Monash University). Samples were submitted to Micromon Genomics (Monash University) for analysis using the 10X Genomics Chromium Single Cell 3′ V3.1 chemistry, as per the manufacturer's instructions. Samples were subjected to 10 cycles of PCR for cDNA amplification and 16 cycles for library amplification. The prepared libraries were sequenced using the MGITech MGISEQ2000RS platform with MGIEase V3 chemistry, on two lanes with 100 bp paired end reads.

Reads were aligned to the chicken GRCg6a reference using CellRanger (v4.0.0, using option: –force-cells 15,000). Owing to the number of reads observed immediately downstream of annotated genes, the gene annotation (from Ensembl release 100, gene biotypes: protein coding, lincRNA and antisense) was edited to include 1000 bp downstream of each gene. Single cell analysis was performed in R using the packages scran (Lun et al., 2016), scater for QC (McCarthy et al., 2017), monocle3 for trajectory analysis (Trapnell et al., 2014) and iSEE for interactive viewing (Rue-Albrecht et al., 2018). Gene names were used for analysis and, where they mapped to multiple Ensembl ids, the Ensembl ID with the highest number of counts was kept. Cells with low total UMI counts (<2000) were excluded. The cell cycle was annotated with cyclone in the scran package (Lun et al., 2016) using the mouse reference from Scialdone et al. (2015) mapped to its one-to-one chicken orthologs.

The top 1000 genes with the highest biological variance were identified with modelGeneVar function of scran (Lun et al., 2016), blocked on the sequencing sample and mitochondrial genes or genes on the Z or W chromosomes were excluded to minimise sex effects. PCA was calculated on these, and the first 15 PCs used to generate a global chicken tSNE layout. Clusters were defined with the walktrap method, on a SNN graph (k=10) (Lun et al., 2016), and cluster identities were determined from gene logFC changes and spatial expression profiles in the Gallus Expression In Situ Hybridization Analysis (GEISHA) database (Bell et al., 2004; Darnell et al., 2007).

Global signalling patterns throughout the chicken forelimb field were examined using the R package CellChat (Jin et al., 2021), where signalling communication networks were constructed based on 1:1 gene orthology with a curated Homo sapiens database and default parameters. Visualizations of pathway and ligand receptor signalling were generated with CellChat and edited with Adobe Illustrator.

Mesodermal cell clusters were subset from the full tSNE for additional focused analyses. Briefly, mesodermal clusters were subset to a new object, and PCA and UMAP dimension reduction were recalculated using the previously determined highly variable genes and corrected for cell cycle effects. Next, the object was imported into Monocle3 (Cao et al., 2019; Trapnell et al., 2014) for clustering (k=4). Cluster labels were confirmed by identifying differentially expressed marker genes through regression analysis implemented in monocle fit_models function, producing distinct 12 sub-clusters covering all known cell types within the developing limb field mesoderm.

Estimations of RNA velocity were calculated where reads were aligned to the reference a second time with STAR solo (v2.7.5) to identify proportions of spliced and unspliced transcripts (Dobin et al., 2013). Velocity analysis was performed with velocyto.R (La Manno et al., 2018), and directional transcriptional velocities between cells were visualized. Lineage trajectories throughout the mesoderm were additionally constructed in Monocle3 using reverse graph embedding (k=4, minimum branch length=15, rann.k=50), which produced four major lineages that originated in E1.5 cells and terminated in E3.5 cells. Lineage bifurcation points were corroborated using known LPM marker gene expression through the plot_cells function, then pseudotime was calculated by selecting the origin of the lineages using the order_cells function. To identify genes that dynamically changed their expression across pseudotime, key lineages throughout the mesoderm were subset using the chose_graph_segments function and graph tests were run to identify lineage-specific, differentially regulated genes and filtered based on Moran's I statistic and q value. Modules of genes that significantly changed across pseudotime were visualized by hierarchical clustering through the R package ComplexHeatmap. Gene expression in individual cells across pseudotime were further visualized using the plot_gene_in_pseuodtime function in Monocle3.

Whole-mount in situ hybridization for spatial mRNA expression was carried out as described previously (Smith et al., 2016) with minor modifications. Briefly, whole HH8-HH22 chicken embryos were fixed overnight in 4% paraformaldehyde, dehydrated in methanol and rehydrated in PBTX (PBS+0.1% Triton X-1000). Tissues were permeabilized in 10 µg/ml proteinase K for up to 1 h, depending upon stage, then re-fixed in glutaraldehyde/4% PFA. Tissues underwent pre-hybridization (50% formamide, 5× SSC, 0.1% Triton X-100, 0.5% CHAPS, 5 mM EDTA, 50 mg/ml heparin, 1 mg/ml yeast RNA and 2% blocking powder) overnight at 65°C. Riboprobe templates were provided as gifts, generated from public sources, or designed and synthesised in-house. Primer sequences and/or source are listed in Table S5. Where applicable, templates were amplified from limb and whole-embryo cDNA using gene-specific primers. Fragments were resolved by 1% agarose electrophoresis, excised and purified using a Nucleospin PCR clean-up kit and subcloned into p-GEM T-easy (Promega). Antisense RNA probes were synthesized using T3, T7 or SP6 RNA polymerases and the DIG-labelling kit (Roche, 11277073910) as per the manufacturer's instructions. Precipitated probes were added to pre-hybridized tissues (∼5 ml/tube) and hybridization was carried out overnight at 65°C. Tissues were then subjected to stringency washes, blocked in BSA, then treated overnight with anti-DIG antibody conjugated with alkaline phosphatase. Tissues were exposed to BCIP/NBT colour reaction at room temperature for up to 3 h (340 mg/ml NBT and 175 mg/ml BCIP in NTMT [100 mM NaCl, 100 mM Tris-HCl (pH 9.5), 50 mM MgCl₂ and 0.1% Tween-20].

Chicken embryos were fixed in 4% PFA/PBS for 15 min at room temperature then cryo-protected in 30% sucrose. Embryos were snap frozen in OCT and 10 mm frozen sections were cut. For immunofluorescent detection of TWIST1, antigen retrieval was performed by the Monash Histology platform. Briefly, slides were heated at 60°C for 30 min followed by retrieval in citrate solution (pH 6). All other sections were left in PBS. For co-detection of TWIST1 and other markers (Table S5) in the LPM, antibody incubations were performed on successive tissue sections. Sections were blocked and permeabilised in 1% Triton X-100, 2% BSA/PBS for 1-2 h at room temperature, then incubated with primary antibody (1:100) in 0.5% Triton X-100, 1% BSA/PBS incubation overnight at 4°C. Secondaries (1:1000) were added the following day and left at room temperature for 1 h. Slides were washed once with DAPI diluted in PBS, followed by two washes in PBS.

Electroporation of chicken ectoderm was performed using custom parameters. Briefly, eggs were incubated for ∼36 h until stage HH8-HH10. Here, a solution containing TOL2 Transposase and Tol2-CAG-GFP-T2A or Tol2-CAG-GFP-T2A-NOG plasmids (available on request) at final concentrations of 1 μg/µl were mixed with 0.1% Fast Green and injected between the vitelline membrane and embryo. Electroporation was performed by placing the negative electrode above the presumptive forelimb field on the right side of the embryo, and positive electrode under the embryo above the yolk, and delivered through three 10 V, 60 ms width and 50 ms spaced pulses (Intracel TSS20 Ovodyne Electroporator). Eggs were sealed with tape and the embryos were incubated for a further 48-72 h. Embryos were then harvested and GFP fluorescence were captured on a fluorescence dissecting microscope (Zeiss SteREO Discovery V12). GFP-positive embryos were then fixed in 4% PFA overnight at 4°C. Electroporations were performed across multiple independent experiments to ensure consistent results.

In situ hybridization was performed on GFP-positive embryos, as described above. For co-detection of in situ gene expression and GFP immunofluorescence, positive embryos were cryosectioned at 10 μm and subjected to antigen retrieval. Sections were processed for immunofluorescence, as above, using an anti-GFP (Table S5) antibody. Tissue sections were imaged on a bright-field or fluorescent microscope (Zeiss Axio Imager A1).

We thank Micromon Genomics for assistance with sequencing design, Monash University Flowcore for flow cytometry and Monash Histology platform for histological processing. We additionally thank Nadia Davidson, Belinda Phipson, Alex Combes and Kieran Short for helpful suggestions for choice of bioinformatics methods, and constructive comments from all members of the Smith Lab.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.200702.reviewer-comments.pdf