A branching gene regulatory network dictating different aspects of a neuronal cell identity

The nervous system contains many different neuronal subtypes and understanding cell fate specification remains a major challenge. Mature neurons differ from each other with respect to a number of different properties, including neurotransmitter identity, electrophysiological properties and axon/dendrite morphology. These different properties are not always linked. For example, neurons producing the same neurotransmitter may project axons to different targets, and neurons projecting to the same target may produce different neurotransmitters. In addition, related neurons may originate from different progenitors and different locations with the CNS. This begs the question of how the developmental regulatory machinery can determine which features to combine in order to establish a particular cell fate. Although a number of studies have focused on the genetic programs determining one specific aspect of a neuronal cell fate, few studies have addressed multiple aspects simultaneously. The exception is studies in C. elegans, which have involved a comprehensive analysis of the regulatory programs dictating most, if not all, aspects of terminal neuron cell fate for several specific neurons (reviewed by Hobert, 2016; Hobert et al., 2010). This has unravelled several intriguing principles underlying neuronal subtype specification. However, it remains unclear how regulatory cascades are integrated to determine all aspects of neuronal subtype identity in animals with thousands of neurons and hundreds of cell fates.

The Apterous (Ap) neurons of the Drosophila ventral nerve cord (VNC) constitute a powerful model for addressing these issues, because of the availability of selective markers and genetic tools, as well as the insight available into their specification that stems from previous extensive work (reviewed by Allan and Thor, 2015). Ap neurons constitute a small subgroup of interneurons that express the Apterous (Ap) LIM-homeodomain factor and the Eyes absent (Eya) transcriptional co-factor and nuclear phosphatase (Fig. 1A) (Lundgren et al., 1995; Miguel-Aliaga et al., 2004). A subset of Ap neurons expresses the Nplp1 neuropeptide: the lateral thoracic Tv1 neurons, part of the thoracic Ap cluster of four cells, and the dorsal medial row of dAp neurons (Fig. 1A,B) (Baumgardt et al., 2007; Park et al., 2004). In line with the different position of the Tv1 and dAp neurons, previous studies have revealed that they originate from two different progenitors (neuroblasts, NBs): NB5-6T and NB4-3, respectively (Baumgardt et al., 2009; Gabilondo et al., 2016). In addition to their expression of the same neuropeptide, Tv1 and dAp neurons also project their axons along a common route, ipsilaterally and anteriorly along the medial-most Fasciclin II (FasII)-positive axon bundle (Fig. 1A,B) (Lundgren et al., 1995). In addition to expressing the Nplp1 neuropeptide, Tv1 and dAp are similar with respect to their secretory properties. Specifically, their neuropeptide expression necessitates the co-expression of a number of genes involved in the processing of neuropeptides, as well as the establishment of secretory large dense-core vesicles (Allan et al., 2005; Park et al., 2004, 2008b). These features are collectively referred to as the regulated secretory pathway (Park and Taghert, 2009). Thus, Tv1 and dAp neurons share three defining features: their axon pathfinding; the expression of the Nplp1 neuropeptide; and the expression of the regulated secretory pathway.

A number of studies have addressed the molecular genetic underpinnings of Tv1 and dAp neuron specification, both with respect to axonal projection, the expression of the Nplp1 neuropeptide and establishment of the regulated secretory pathway. These studies have identified two distinct spatio-temporal combinatorial transcription factor codes, one acting in NB5-6T and the other in NB4-3, that dictate Tv1 and dAp fate, respectively (Gabilondo et al., 2016). In NB5-6T, these upstream cues involve the ladybird early (lbe) gene, which activates the collier gene (col; kn – Flybase), encoding a COE/EBF transcription factor (Fig. 1B) (Baumgardt et al., 2009; Gabilondo et al., 2016; Karlsson et al., 2010). Lbe and Col then trigger a feed-forward loop (FFL) consisting of Ap, Eya and the Dimmed (Dimm) bHLH transcription factor, which together activate the Nplp1 gene (Baumgardt et al., 2007; Stratmann and Thor, 2017). In addition, Ap and Eya have been found to play crucial roles in Tv1 and dAp axon pathfinding (Lundgren et al., 1995; Miguel-Aliaga et al., 2004), while Dimm appears to be necessary and sufficient for establishing the regulated secretory pathway (Allan et al., 2005; Hamanaka et al., 2010; Hewes et al., 2006; Park et al., 2011, 2004, 2014, 2008a,b). These genetic studies have shed light upon the cell specification of Tv1 and dAp neurons with respect to the neuropeptide Nplp1 and the regulated secretory pathway. However, they have not addressed the role of each regulator with respect to axon pathfinding, or resolved the full gene expression signature of these neurons regarding regulatory cascades and effector outputs.

To address these issues, we focus here on the Tv1 neurons, which are generated by NB5-6T. We find that both col and lbe are crucial for Tv1 axon pathfinding. In addition, co-misexpression of col/lbe is sufficient to ectopically re-route other neurons into the Tv1 fascicle. The strength of this effect prompted us to conduct transcriptome analysis, attempting to unravel the full Tv1 signature. As anticipated, this identified ap, eya and dimm as being transcriptionally upregulated. We also observed transcriptional upregulation of many genes involved in the regulated secretory pathway, presumably as an effect of dimm upregulation. In addition to ap, eya and dimm, we observed altered expression of two additional transcription factor genes: nerfin-1 and Dbx, as well as of a number of axon pathfinding genes, including uncoordinated-104 (unc-104), golden goal (gogo) and futsch. Mutant and misexpression analysis of the genes identified by transcriptomic analysis revealed different effects upon Tv1 neuron specification. In line with genetic analysis (herein and previously conducted), co-misexpression of ap/eya followed by transcriptome analysis revealed a smaller number of genes that were mis-regulated. Together with previous studies on Tv1 neurons, these findings point to a branching neuronal cell fate specification model, where the different regulators play partly overlapping, but also unique roles with respect to neurotransmitter, axon pathfinding and secretory properties. Such branching regulatory networks are likely to be prevalent in more-complex nervous systems.

The activation of the Nplp1 neuropeptide in Tv1 and dAp neurons is controlled by a shared coherent FFL, consisting of col, ap, eya and dimm (Baumgardt et al., 2007; Gabilondo et al., 2016). Tv1 and dAp neurons are generated by two different NBs, NB5-6T and NB4-3, respectively, and, hence, the FFL is triggered by two different upstream spatio-temporal combinatorial codes (Baumgardt et al., 2009; Gabilondo et al., 2016; Karlsson et al., 2010). In NB5-6T these upstream cues include the homeobox gene lbe (Fig. 1C).

The Tv1 and dAp neurons both project axons ipsilaterally and anteriorly, into a common fascicle that overlaps with the medial-most bundle identified by the axon fascicle marker FasII (Fig. 1A,B) (Lundgren et al., 1995). ap and eya have both been found to control proper axon pathfinding of Tv1/dAp neurons (Lundgren et al., 1995; Miguel-Aliaga et al., 2004), whereas dimm is not involved in axon pathfinding (Allan et al., 2005). The roles of lbe and col with respect to axon pathfinding have not been addressed previously. To unravel the roles of the identified regulators in axon pathfinding, we focus on the Tv1 neurons.

We have previously found that combinatorial misexpression of col and lbe was sufficient to trigger widespread specification of ectopic Tv1 neurons, evident by ectopic Ap, Eya, Dimm and Nplp1 expression (Gabilondo et al., 2016). We repeated this experiment, using the pan-neural elav-Gal4 driver, and confirmed ectopic expression of Ap, Eya, Dimm and Nplp1 (Fig. 1D,E). As anticipated from previous studies (Baumgardt et al., 2009, 2007; Gabilondo et al., 2016; Stratmann et al., 2016), co-misexpression of col and lbe did not activate expression of the FMRFa neuropeptide, expressed by the Tv4 and SE2 neurons, but rather downregulated its expression (Fig. S1A-H). In addition, based upon the Nplp1 expression, ectopic neurons project axons into a common medial fascicle, similar in location to the normal Tv1/dAp fascicle (Fig. 1D,E). This axon projection effect was also evident when staining for FasII, revealing axon re-routing into the medial-most FasII-positive bundle, and a reduction of the number of FasII-positive bundles (Fig. 1F,G, Fig. S2E-K).

Although both ap and eya are necessary for Nplp1 expression and proper axon pathfinding of Nplp1-positive neurons, combinatorial misexpression of ap/eya was not sufficient to trigger ectopic Nplp1 expression (Fig. S2A,B) (Baumgardt et al., 2007). In addition, we did not observe any apparent changes to the FasII-positive axon scaffold (Fig. S2C,D). The striking difference in axon pathfinding effects when comparing col/lbe co-misexpression with ap/eya co-misexpression indicates that col and lbe have roles in axon pathfinding independent of their activation of ap and eya.

The disorganized FasII-positive axon scaffold observed in col/lbe co-misexpression from the pan-neural elav-Gal4 driver raised the possibility that col/lbe may be chiefly activating Nplp1 in neurons already projecting in a medial pathway, rather than simultaneously triggering Nplp1 expression and axon projection changes. To address this, we misexpressed col/lbe using the more-restricted wg-Gal4 driver, which expresses in only one of the seven rows of NBs: row 5 (Chu-LaGraff and Doe, 1993). To reveal axonal projections, we co-expressed the membrane-targeted UAS-EGFP^F transgene (Allan et al., 2003). In line with previous mapping studies of NB lineages in the embryonic VNC (Schmid et al., 1999; Schmidt et al., 1997), in control, neurons generated from row 5 NBs project axons along several different pathways, including several of the longitudinal FasII-positive bundles and motor nerves (Fig. 2A,C). Strikingly, in wg>col,lbe co-misexpressing embryos, we observed extensive ectopic projections into the medial-most FasII-positive bundle, as well as an apparent complete loss of motor nerve projections (Fig. 2B,D, Figs S3A-N, S4A-C). This effect was evident already at stage (St) 16 but became more pronounced by air-filled trachea (AFT) stage (Fig. 2A-D, Fig. S3A-N).

To determine the individual effects of col and lbe, we misexpressed them alone and in combination, and quantified the extent of axon re-routing by comparing the integrated GFP-intensity in the medial-most FasII-positive bundle in the region between the commissures (Fig. S5K). We found that col single misexpression had a small but significant effect (Fig. S5A,B,E). In contrast, single lbe misexpression resulted in a significant reduction of axons in the medial-most FasII-positive bundle (Fig. S5C,E). The reduction of the medial projection in lbe single misexpression may be due to the effect of lbe acting without col. Specifically, early-born NB5-6T neurons, which express lbe and not col (Baumgardt et al., 2009, 2007; Gabilondo et al., 2016), project to the lateral-most FasII-positive bundle, whereas late-born neurons (Ap neurons), which express col and lbe, project to the medial-most bundle (Fig. 3A-E). col/lbe co-misexpression resulted in striking combinatorial effects, with more than a fourfold increase in GFP signal in the medial-most FasII-positive bundle, when compared with control (Fig. S5D,E).

In spite of the lack of obvious effects upon the FasII-positive scaffold in response to ap/eya co-misexpression, when driven from elav-Gal4, the higher resolution of the wg-Gal4, UAS-EGFP^F driver set-up prompted us to re-address the effects of ap and eya. Indeed, using this driver-marker combination, we observed weak, but significant, re-routing to the medial-most FasII-positive bundle in both ap and eya single misexpression (Fig. S5F-H,J). These effects were enhanced in the ap/eya co-misexpression embryos (Fig. S5I,J). However, the ap/eya combinatorial effect was weaker than that of col/lbe co-misexpression, again supporting the notion that col and lbe have axon pathfinding roles outside of their activation of ap and eya.

Next, we addressed the mutant effects of col and lbe upon Ap neurons and Tv1 axon pathfinding. To achieve this, we could not use ap, eya, dimm or Nplp1 axonal reporters, because these genes are all downregulated in col and lbe mutants (Baumgardt et al., 2007; Gabilondo et al., 2016). Hence, we made use of the lbe(K)-EGFP transgenic reporter, which labels the entire NB5-6T lineage (Baumgardt et al., 2009). However, because lbe(K)-GFP labels the entire lineage, we first needed to map the subset of NB5-6T projections that correspond to the Ap neurons (including Tv1). Hence, we combined lbe(K)-GFP with ap-tau-lacZ (Lundgren et al., 1995), and co-stained for GFP, β-gal and FasII. This analysis revealed the development of the entire NB5-6T lineage projections, the Ap neuron projections and their relation to the FasII-positive scaffold (Fig. 3A-E). We found that the earlier-born NB5-6T neurons projected into the lateral-most FasII-positive bundle, as well as across the anterior commissure, whereas Ap axons can be observed at St16 joining the medial-most bundle (Fig. 3A-E). With this information at hand, we introduced the lbe(K)-GFP transgene into the col and lbe mutants, and assayed the effects upon the Ap axons, including Tv1. In both col and lbe mutants, we observed a failure of the Ap axons to project into the medial-most FasII-positive bundle, apparent both at St16, St17 and stage AFT (Fig. 3F-H, Fig. S6A-E). In line with the finding that early-born NB5-6T neurons, which express lbe but not col, project laterally or across the anterior commissure (Fig. 3A-E), we observed defects in these projections in lbe mutants, whereas col only affected the medial-most projections (Fig. 3F,G).

The mutant and misexpression analyses demonstrate that col and lbe are necessary, and in part sufficient, for Tv1 axon pathfinding. Similarly, ap and eya are also necessary and, to a lesser extent, sufficient for Tv1 axon pathfinding. This raises the issue of whether or not lbe and col act only via ap and eya or whether they have roles outside this function. To address this, we misexpressed col in the NB5-6 lineage, using a lbe(K)-Gal4, UAS-EGFP^F combination, in control as well as in ap and eya mutants. Because of the effect of single misexpression of lbe from wg-Gal4, i.e. routing of axons into the lateral-most FasII-positive bundle, we overexpressed lbe in ap or eya mutants. We analysed abdominal segments, because the lack of Ap neurons in the NB5-6A lineages (Karlsson et al., 2010) results in the absence of axons from this lineage in the medial-most FasII-positive bundle (Fig. S7A). In control, similar to the effects of elav-Gal4 and wg-Gal4, misexpression of col from lbe(K)-Gal4 triggers ectopic axon projections into the medial-most FasII-positive bundle (Fig. S7B). Strikingly, in both ap and eya mutants, col is still able to trigger these ectopic axon projections (Fig. S7C-F).

To summarize, col, lbe, ap and eya are all crucial for proper axon pathfinding of Ap neurons, including the Tv1 neuron. Although col/lbe co-misexpression robustly activates Ap and Eya, the weaker effects of ap/eya co-misexpression and the misexpression effects of col in ap or eya mutant backgrounds indicate that col and lbe have axon pathfinding roles in addition to their activation of ap and eya.

Co-misexpression of col/lbe results in robust reprogramming of many neurons to Tv1 cell fate, as apparent from ectopic expression of the Tv1 regulatory cascade (Ap, Eya and Dimm), the Nplp1 neurotransmitter and Tv1-like axon pathfinding. These results suggested that col/lbe co-misexpression followed by transcriptome analysis of embryonic RNA could allow for the identification of the entire Tv1 cell signature, as well as identify novel regulators and effectors involved in Tv1 fate. To this end, we co-misexpressed lbe/col and analysed the transcriptome, using whole-embryo RNA-seq, both at St16 and stage AFT.

As anticipated, we found that the entire known Tv1 regulatory cascade was activated, including ap, eya and dimm, as well as Nplp1, both at St16 and stage AFT (Fig. S8A,B). As noted above, owing to the temporal progression of the Ap neuron regulatory cascade (Baumgardt et al., 2009, 2007; Gabilondo et al., 2016; Stratmann et al., 2016), col/lbe co-misexpression does not trigger ectopic expression of the FMRFa neuropeptide (Fig. S1A-N) (Gabilondo et al., 2016), which is expressed by the Tv4 Ap neurons. In line with these previous findings, we did not observe upregulation of the FMRFa gene (Fig. S8A,B).

Previous studies have demonstrated that dimm controls the regulated secretory pathway, being both necessary and sufficient to directly activate a number of genes in this pathway (Allan et al., 2005; Hamanaka et al., 2010; Hewes et al., 2006; Park et al., 2011, 2004, 2014, 2008a,b). Because col/lbe activates dimm expression, we analysed the expression of 21 of the previously identified dimm target genes in the regulated secretory pathway. This analysis revealed that col/lbe co-misexpression significantly activated expression of 14 of these 21 genes, with the remaining six showing a trend upwards, albeit non-significant (Fig. S8C).

Next, we turned to the axon pathfinding aspect of Tv1 cell fate. Because the aforementioned results showed that not only col/lbe but also ap/eya are involved in the axon pathfinding of Tv1 neurons (this study; Lundgren et al., 1995; Miguel-Aliaga et al., 2004), we also analysed the transcriptome for combinatorial and single misexpression of all four genes. Focusing first on co-misexpression of col/lbe versus ap/eya, at stage AFT, this analysis revealed twofold changes (2FC) in expression of 1075 and 650 genes, respectively (Fig. 4A,B). A smaller number of genes showed altered expression at St16 (Fig. S9A,E). There was considerable overlap in the gene expression profiles, with 458 genes affected by both col/lbe and ap/eya (Figs S9E, S10B). Gene ontology (GO) filtering of the overlapping set of genes revealed a diverse set of molecular functions, biological processes and cellular components (Fig. S10A-C).

GO filtering for the 416 Drosophila genes identified as ‘neuron projection development’ revealed that col/lbe and ap/eya co-misexpression altered expression of 13 and seven of these genes, respectively (Fig. 4A,B). As anticipated, col itself (knot; kn) was identified in the col/lbe co-misexpression samples, as well as ap and eya in the ap/eya co-misexpression samples (Fig. 4A,B). Surveying the overlap in gene expression between col/lbe and ap/eya within the GO category ‘neuron projection development’, we identified five genes whose expression was changed twofold in both genotypes, indicating that those factors are regulated by col/lbe via ap/eya. In addition to ap and eya themselves, these were: jelly belly, which was up-regulated; and ftz-f1 and Neurotrophin 1 (NT1), which were downregulated (Fig. 4C).

Next, we focused on the eight ‘neuron projection development’ genes regulated by col/lbe and not by ap/eya. To address the combinatorial action of col/lbe, we chose to include the RNA-seq results from the col and lbe single misexpression in the heatmap, regardless of their fold-change value. We find a strong combinatorial effect in the RNA-seq data of col/lbe on nerfin-1 (zinc finger, C2H2-type; vertebrate: NSM1/2), which is upregulated by lbe/col co-misexpression, but not in the col and lbe single misexpression experiments (Fig. 4D). This effect was confirmed by in situ hybridization, which revealed an apparent upregulation of nerfin-1 in the elav>col,lbe embryonic VNC (Fig. 5I-J). Previous studies of nerfin-1 have revealed that its mRNA is expressed by most, if not all, VNC NBs, with Nerfin-1 protein being restricted to a subset of NBs and their daughter cells (Kuzin et al., 2005; Stivers et al., 2000). Both nerfin-1 mutants and misexpression resulted in axon-pathfinding defects in the developing VNC (Kuzin et al., 2005). To test the roles of nerfin-1 in Tv1 axon pathfinding, we analysed mutants and misexpression at stage AFT, and monitored potential pathfinding by staining for FasII and Nplp1. nerfin-1 mutants showed strongly disorganized FasII-positive axon bundles and, in addition, disorganized or absent Nplp1-positive axons (Fig. 4E,F). Misexpression of nerfin-1 from elav-Gal4 showed a milder phenotype with respect to the FasII-positive bundles, but disorganized Nplp1 axon projections, which fail to project towards the medial-most FasII-positive bundle (Fig. 4G). Strikingly, however, Nplp1 expression was apparently unaffected in nerfin-1 mutants or misexpression (Fig. 4J). On analysing Nerfin-1 expression, we found that Nerfin-1 is expressed in the NB5-6T lineage and Ap cluster neurons, evident by overlap with lbe(K)-GFP and Eya, at St14 and late St15 (StL15). At St16, Nerfin-1 expression is absent from the Ap cluster (Fig. S11A-D). This prompted us to test whether Nerfin-1 can be activated or maintained in the NB5-6T lineage at St16, when col/lbe are co-misexpressed from the NB5-6 lineage-specific driver lbe(K)-Gal4. Indeed, we find that col/lbe co-misexpression results in a minor, yet significant, increase of Nerfin-1-positive cells in the NB5-6T lineage (Fig. 5A,B,G). Next, we quantified Nerfin-1-positive cells in col and lbe mutants in the NB5-6T lineage. Based upon the dynamic Nerfin-1 expression, we decided to assess Nerfin-1 cell numbers at an earlier stage than St16, i.e. StL15. However, quantification of Nerfin-1-positive cells in both col and lbe mutants did not show any significant differences in cell number (Fig. 5C,D,G). We further analysed the expression of Nerfin-1 in ap and eya mutants, and, in line with the lack of upregulation of nerfin-1 by elav>ap,eya in the RNA-seq data, mutant analysis did not reveal any apparent effects on Nerfin-1 expression (Fig. S12A-H).

In contrast to nerfin-1, the RNA-seq revealed that Dbx (homeobox factor; vertebrate: DBX1/2, BARX1) was downregulated in lbe/col co-misexpression, and also strongly affected by a single misexpression of lbe (Fig. 4D). Similarly, RNA-seq at St16 revealed a strong combinatorial effect of lbe/col on Dbx (Fig. S9D). In line with these findings, col/lbe co-misexpression, from elav-Gal4, resulted in significantly decreased numbers of Dbx cells, both in thoracic and abdominal segments (Fig. 6A,B,G). However, col or lbe mutants did not show significant differences in Dbx-positive cells, in either thoracic or abdominal segments (Fig. 6C-F,H,I). Previous studies revealed that Dbx mutants did not show an apparent axon pathfinding phenotype, whereas Dbx misexpression led to a decrease in the ability of many neurons to extend axons (Lacin et al., 2009). In agreement with this, we observed apparently intact FasII-positive bundles in Dbx mutants. In contrast, there were severe effects on the axon projections of the Nplp1 axons, characterized by a failure to project along the medial-most FasII-positive bundle (Fig. 4H). Dbx misexpression, under control of elav-Gal4, showed no apparent effect on the FasII-positive bundles, but projection defects of Nplp1 axons along the medial-most FasII-positive bundle (Fig. 4I). Quantification of Nplp1-positive cells in the Dbx mutant and misexpression embryos surprisingly revealed a partial loss of Nplp1-positive cells in both cases (Fig. 4J). These findings prompted us to further address the effect of Dbx on Tv1 cells, by analysing additional Tv1/Ap cluster markers, at stage AFT. Looking first at Col, Ap, Eya and Nplp1, we noted that whereas Ap clusters were present in each hemi-segment of both Dbx and misexpression mutants, there was an apparent mis-specification of the Ap cluster cells, with a frequent loss of Nplp1 expression (Fig. S13A-G). However, we also noted some cases of an extra Nplp1-expressing cell in the Ap cluster (Fig. S13A-G). This notion was confirmed by analysis of Dimm, together with Ap, Eya and Nplp1, which also revealed that both Dbx and misexpression mutants displayed a mixed phenotype with either loss of, or extra, Dimm cells in the Ap cluster (Fig. S14A-E). Previous studies on Dbx revealed it to be expressed in neurons stemming from five NBs, none of which corresponds to NB5-6 or NB4-3 (Lacin et al., 2009) (the NBs that generate Ap neurons); indeed, we observed no apparent overlap between Dbx and Ap in postmitotic cells (Fig. 6A). To determine whether Dbx is expressed earlier in the NB5-6T lineage, perhaps in the NB itself, we analysed Dbx expression in NB5-6T during St12-St14. However, we did not, detect any apparent expression in the lineage at these stages (Fig. S15A-F).

The scarcity of axon pathfinding genes in the 2FC group in the RNA-seq data prompted us to analyse genes with a fold-change of ≥1.5. Not surprisingly, this expanded the number of genes with significantly altered expression to 3174 for elav>col/lbe and 2194 for elav>ap/eya (Fig. S16A,B). This also expanded the subset of genes belonging to the ‘neuron projection development’ GO class: 65 for col/lbe and 52 for ap/eya (Fig. S16A,B). Comparing these two subsets to each other revealed extensive overlap (Fig. S16C). Given the potent role of col/lbe in controlling Tv1-type axon projection, in part independent of ap/eya, we focused further upon the ‘neuron projection development’ genes altered by col/lbe but not ap/eya, and in particular those most effectively regulated by col/lbe co-misexpression as opposed to single misexpression. In this group, several of the upregulated genes stood out based upon their previously described roles in ‘neuron projection development’, including golden goal (gogo), uncoordinated-104 (unc-104) and futsch (Fig. S16D) (Berger et al., 2008; Hummel et al., 2000; Medina et al., 2006; Roos et al., 2000; Tomasi et al., 2008). To address the possible role of these three genes during Tv1 axon pathfinding, we analysed mutants with the Nplp1 and FasII markers. None of the mutants displayed loss of Nplp1 expression in cell bodies (Tv1 or dAp). However, gogo and futsch mutants both displayed an apparent disorganization of Tv1 axons (Fig. S17A,B,D,E). While unc-104 mutants displayed Nplp1 expression in cell bodies, there was a complete loss of Nplp1-positive axonal projections (Fig. S17C). The unc-104 phenotype prompted us to introduce the Nplp1-CRM-myrRFP transgene into the unc-104 background, which revealed that Tv1 axons did indeed project normally in this mutant (Fig. S17F,G). This finding is in line with previous studies, demonstrating that unc-104, which encodes a Kinesin-like molecule, is crucial for dense core vesicle axon transport and, hence, neuropeptide transport (Barkus et al., 2008; Lim et al., 2017).

Previous ChIP-seq analysis of Col, from whole mid-embryo chromatin, identified ∼400 genes (de Taffin et al., 2015). We compared this gene list with our RNA-seq lists for elav>col, both the 1.5-fold change and 2.0-fold change, at St16 and AFT, and found a number of overlapping genes (Fig. S18A-D). These included nerfin-1, Dbx and ap. In addition, the Col ChIP-seq study also identified eya as a direct target of Col (de Taffin et al., 2015). Previous gene-focused studies have also provided support for direct binding of Lbe on col, Col on ap, eya and dimm, Ap on dimm and Nplp1, and Dimm on Nplp1 (Stratmann and Thor, 2017). In addition, genome-wide analysis of Dimm DNA binding (Hadžić et al., 2015) revealed direct interaction with many of the regulated secretory pathway genes (genes in red in Fig. S8), and, furthermore, the direct regulation of Phm has been observed in directed studies (Park et al., 2008a). Hence, a substantial part of the regulatory cascade for Tv1 cell fate apparently involves direct gene regulation.

A previous study addressing Lbe target genes in muscle cells, by using muscle-specific RNAi and muscle-specific misexpression, resulting in the identification of 117 genes with a high fold change (Junion et al., 2007). We compared these results with our CNS-targeted RNA-seq data for elav>lbe, both the 1.5-fold change and 2.0-fold change, at St16 and AFT. In all cases, we noted minimal overlap between the Lbe CNS and muscle target genes (Fig. S19A-D). Not surprisingly, this indicates that the context (CNS versus muscles) greatly influences target gene selection by Lbe.

The identity of the Tv1 neuron is specified by a genetic cascade that initially commences with a core pathway FFL: lbe→col→ap/eya (Baumgardt et al., 2007; Gabilondo et al., 2016). It then forks into three separate branches, with the first branch being a continued FFL of col→ap/eya→dimm, which activates Nplp1 neuropeptide expression (Baumgardt et al., 2007; Gabilondo et al., 2016). Recent studies indicate that the entire lbe→col→ap/eya→dimm→Nplp1 FFL is executed via sequential DNA binding by the TFs onto each enhancer in the cascade (Stratmann and Thor, 2017).

A second branch consists of dimm alone, which is necessary and sufficient to activate the regulated secretory pathway, including the neuropeptide-processing machinery and the large dense-core vesicle machinery (Allan et al., 2005; Hamanaka et al., 2010; Hewes et al., 2006; Park et al., 2011). Dimm directly activates Phm, and possibly syta and sytb (Park et al., 2014, 2008a), and acts directly on many other genes in the regulated secretory pathway (Hadžić et al., 2015). Evidence points to Dimm also acting directly on Nplp1 (Stratmann and Thor, 2017), as well as on other neuropeptide genes (Gauthier and Hewes, 2006), and it appears to play a dual role in both activating the regulated secretory pathway and acting combinatorially to activate neuropeptide gene expression (Baumgardt et al., 2007). We find that col/lbe co-misexpression, which ectopically activates dimm in the CNS, results in the activation of most of the previously identified Dimm targets in the regulated secretory pathway.

In contrast to the central role of dimm with respect to the regulated secretory pathway, dimm does not affect axon pathfinding of Tv1 neurons (Allan et al., 2005). Instead, axon pathfinding is controlled by the third branch, where the col→ap/eya and lbe→col→nerfin-1 FFLs are important for Tv1 axon pathfinding, ensuring that Tv1 axons project towards the midline and join the medial-most FasII-positive fascicle. In line with this model, nerfin-1 does not affect Nplp1 expression, only pathfinding. The transcriptome and genetic analysis further suggests that the axon pathfinding FFL is itself divided into sub-branches, where col/lbe activate the effectors unc-104, gogo and futsch, while ap/eya activate (for example) jeb and dnc. Because col/lbe activate ap/eya, we cannot rule out a model whereby genes activated by col/lbe but not by ap/eya still require the combinatorial action of col/lbe/ap/eya. However, the lack of ectopic and mutant effects on (for example) nerfin-1 argues against such a model, at least regarding this downstream target. Thus, the specification of the Tv1 neuron commences with a core initial terminal selector FFL pathway, which branches out into different FLLs, each of which executes different sub-routines that are crucial for final Tv1 identity (Fig. 7). The results from both genome-wide and directed DNA-binding studies (de Taffin et al., 2015; Hadžić et al., 2015; Park et al., 2008a; Stratmann and Thor, 2017) support the hypothesis that these FFLs play out via sequential direct DNA-binding of TFs onto the downstream genes (genes in red in Fig. 7 are predicted to be direct targets for one or several of the upstream TFs).

The pleiotropic effects of Dbx on Tv1 neuron specification, with both mutants and misexpression affecting several Tv1 properties, including Nplp1 expression, Ap cluster composition and axon pathfinding, makes it difficult at this stage to precisely place Dbx within the Tv1 cascade. However, this pleiotropic action is reminiscent of other genes acting in the NB5-6T, such as Kruppel, squeeze, nab and seven up, all of which are dynamically expressed during the generation of Ap cluster neurons (Allan et al., 2005; Allan et al., 2003; Baumgardt et al., 2009; Benito-Sipos et al., 2011; Stratmann et al., 2016; Terriente Felix et al., 2007). The surprising apparent lack of Dbx expression in Tv1 neurons or in the NB5-6T lineage may be due to the fact that Dbx is expressed even earlier than was assayed here, or that certain Dbx protein variants are not detected with the antibody used. Alternatively, given the lack of apparent Dbx expression in the NB5-6T lineage, and as the misexpression of Dbx was conducted using elav-Gal4, a pan-neural driver, it is possible that both the mutant and misexpression phenotypes of Dbx reflect non cell-autonomous effects.

One, perhaps unexpected, finding was that although the transcriptome analysis of col/lbe co-misexpression did reveal twofold upregulation of ap, eya, dimm and nerfin-1, it did not reveal upregulation of any additional TFs. The absence of such TFs may suggest that the identified regulators for axon pathfinding (lbe→co→ap/eya/nerfin-1), for Nplp1 (col→ap/eya→dimm) and for regulated secretory pathway (dimm) may indeed constitute the main regulators of these subroutines in Tv1 neurons. This notion is also supported by a recent forward genetic screen, scoring for expression of FMRFa-eGFP in the Tv4 Ap cluster neurons, which, although it re-identified lbe, col and eya, did not identify any other late-acting TFs (Bivik et al., 2015). Another unexpected result was that we did not observe twofold upregulation of any obvious axon pathfinding genes, but did identify such genes in the lower (1.5-fold change) group. This could stem from their generally broader expression, which may reduce the fold-change effects in our whole embryo RNA analyses. However, it should also be noted that filtering solely on the GO term ‘neuron projection development’ may fail to identify some bona fide axon pathfinding genes, and also include genes that are not directly involved in axon pathfinding. Hence, future analysis of the 1075 genes with twofold change and the 3174 genes with 1.5-fold change, affected by col/lbe co-misexpression, may help identify additional axon pathfinding genes, as well as further refine the Tv1 neuron genetic cascade.

The specification of the terminal and unique identity of any given neuron could conceivably be controlled by a number of different mechanisms (Allan and Thor, 2015; Hobert, 2008, 2016; Hobert et al., 2010). Theoretically, one can envision two extreme models: (1) that all aspects of the identity of a neuron could be determined by a single TF; or (2) that all aspects of the fate of a neuron are combinatorially regulated, and that, in addition, there could be considerable branching in the regulatory networks.

Evidence for the first model stems in particular from C. elegans, where the TF CHE-1 determines most, if not all, aspects of the identity of a specific gustatory neuron (Etchberger et al., 2007). These studies, and related findings, led the authors to coin the term terminal selectors (Wenick and Hobert, 2004). Similarly, in mammals, the TF Pet-1 acts as a terminal selector, being a key regulator of serotonin neuron specification, controlling genes needed for serotonin synthesis, reuptake and vesicular transport, as well as several other unique features of serotonin neurons (reviewed by Spencer and Deneris, 2017). This scenario, of one factor dictating the entire fate of a neuron, may perhaps have been favoured during evolution, as it may help reduce regulatory complexity.

Support for the second model emerges from the common anatomical scenario where neurons expressing the same neurotransmitter may project axons to different targets. There are many examples of such divergence between transmitter and morphology, e.g. different groups of dopaminergic and serotonergic neurons in the mammalian CNS (Björklund and Dunnett, 2007; Gaspar and Lillesaar, 2012), as well as neuropeptide-producing neurons in many systems (Hokfelt et al., 2000; Park et al., 2008b). Hence, there is an apparent dissociation between neurotransmitter identity and axonal projections, a biological reality that is not easily accommodated by one single TF dictating the entire cell fate. Indeed, studies in C. elegans have revealed that the LIM-homeobox gene ceh-14 only regulates the neuropeptide expression in BDU neurons, but not other aspects of this cell fate (Gordon and Hobert, 2015). Based upon the results reported here, and on previous studies, we find that Tv1 fate is not controlled by one single terminal selector TF, but rather by several combinatorial TFs codes and one single TF (dimm), each of which launch different subroutines that are crucial for Tv1 cell fate.

Having the regulatory code for neurotransmitter identity diverge from that controlling axonal pathfinding may better accommodate, even drive, the diverging identities of many neuronal sub-classes during evolution. Furthermore, it would allow for association and disassociation between different terminal neuronal cell fate features. These include association/dissociation of neurotransmitter and cellular morphology, simultaneous or separate production of neuropeptides and small neurotransmitters, coupling-uncoupling of neuropeptide expression with the regulated secretory pathway (scaling). It is tempting to speculate that evolution may have progressed from the basic scenario of one-selector-one-fate to a branching landscape, where every sub-routine of a neuron's property is combinatorially controlled, and where several codes act simultaneously to govern different sub-routines.

Driver lines used were wg⁷⁵⁸ (referred to as wg^Gal4; provided by Konrad Basler, University of Zurich, Switzerland), lbe(K)-Gal4,UAS-EGFP^F [double, i.e. on both chromosome (chr) II and chr III], lbe(K)-EGFP (Ulvklo et al., 2012), apC2.1-tau-lacZ (chr III) (Lundgren et al., 1995) and elav^C155=elav-Gal4 (Bloomington Drosophila Stock Center 458).

Mutant lines used were: lbe^12C005 (BL#59385); Df(lbl-lbe)B44 (referred to as lbe^Df; provided by Krzysztof Jagla, Clermont University, Clermont-Ferrand, France); col¹ and col³ (Crozatier and Vincent, 1999) (provided by Alain Vincent, University Paul Sabatier, Toulouse, France); Df(2L)BSC354 (referred to as eya^Df) (BL#24378); eya^cli-IID (BL#3280); ap^P44 (Bourgouin et al., 1992); nerfin-1^Df (BL#7565); nerfin-1^D159 (provided by Alexander Kuzin and Ward Odenwald, NINDS, Bethesda, MD, USA); Dbx^D48 (BL#56824); gogo^H1657 and gogo^D1600 (provided by Takashi Suzuki, Tokyo Institute of Technology, Yokohama, Japan); futsch^N94 (BL#8805); unc-104^P350 (BL#24630) and unc-104^R757 (BL#24631).

Misexpression lines used were: UAS-col-HA and UAS-myc-lbe (Stratmann et al., 2016); UAS-Flag-ap and UAS-eya-HPC4 (this study); UAS-nerfin-3xHA (FlyORF #F000461); and UAS-Dbx-3xHA (FlyORF # F001875).

Nplp1-CRM-myrRFP-HA (this study) flies were maintained over GFP- or YFP-marked balancer chromosomes. Control OregonR was used as wild type. Staging of embryos was performed according to Campos-Ortega and Hartenstein (Campos-Ortega and Hartenstein, 1985).

Primary antibodies used were: chicken anti-proNplp1 (1:1000), rabbit anti-proFMRFa (1:1000) and guinea pig anti-Dimm (1:1000) (Baumgardt et al., 2007); mouse mAb anti-Eya 10H6 (1:250) and mAb anti-FasII ID4 (1:30) (Developmental Studies Hybridoma Bank); rabbit anti-Ap (1:1000) (Bieli et al., 2015) (provided by Dimi Bieli and Markus Affolter); guinea pig anti-Dbx (1:500) (provided by James B. Skeath, Washington University, St Louis, MO, USA); guinea pig anti-Nerfin-1 (1:500) (provided by Ward Odenwald, NINDS, Bethesda, MD, USA); and goat anti-GFP (1:1000) (ab6673, Abcam).

A nerfin-1 full-length cDNA (LD18634) was obtained from the Drosophila Genomics Resource Center. The authenticity of the cDNA was confirmed by sequencing. The cDNA was linearized using XbaI, and a ribo-probe generated by T7 polymerase. In situ hybridization was conducted as previously described (Baumgardt et al., 2007).

UAS-Flag-ap and UAS-eya-HPC4 were generated by the codon optimization for expression in Drosophila (www.kazusa.or.jp/codon/). EcoRI site and consensus start codon (Cavener and Ray, 1991) was added to the 5′ end. Three different stop codons (amb, och, opa) followed by an XbaI site were added to the 3′ end (see Table S1 for DNA sequences). DNAs were generated by gene synthesis (Genscript) and cloned into pUASattB (Bischof et al., 2007) as EcoRI/XbaI fragments. DNAs were injected into landing site strains BL#9723 (28E) for UAS-Flag-ap and BL#9744 (89E) for UAS-eya-HPC4 (BestGene).

To generate the Nplp1-CRM-myrRFP-HA.attB landing site construct, first the pEGFP.attB landing site vector (Bischof et al., 2013) (provided by Konrad Basler and Johannes Bischof, University of Zurich, Switzerland) was cut with NdeI and XmaJI in order to remove the GFP cassette. Cutting with NdeI creates a TA overhang compatible with a VspI restriction site, which also creates TA overhangs. A VspI-myrRFP-HA-XmaJI sequence containing the L21 sequence upstream the ATG start codon (Pfeiffer et al., 2012) was synthesized at Genescript and delivered in a pUC57 vector (sequences listed in Table S2). The VspI-myr::RFP-HA-XmaJI construct was cut from the pUC57 vector with VspI and XmaJI, and ligated into the empty attB landing site vector to create a myrRFP-HA.attB vector. The Nplp1 enhancer (Stratmann and Thor, 2017) was generated via PCR on OregonR DNA (primers listed in Table S2) and cloned into the TOPO pCR2.1-TOPO TA vector according to the manufacturer's protocol (Invitrogen). In order to clone the Nplp1-CRM into the myrRFP.attB landing site vector, the myrRFP.attB vector was digested with BglII and XbaI, and the Nplp1-CRM was cut from the TOPO vector with BamHI and XbaI and subsequently ligated upstream the myrRFP-HA cassette. The Nplp1-CRM-myrRFP-HA.attB landing site construct integrated into the fly genome via site directed phiC31 mediated integration (Bischof et al., 2007) at BestGene. Stocks for injection were BL#8622 and BL#9723.

Per each sample ∼200 embryos of St16 and stage AFT were collected and RNA was isolated from whole embryos (Qiagen RNeasy Mini kit, 74104). For each genotype and stage, three independent biological replicates were prepared. RNA-sequencing library preparation was carried out using the NEBNext Ultra RNA Library Prep Kit for Illumina by following manufacturer's recommendations (NEB). The samples were sequenced on the Illumina Hiseq 2500 platform using a 1×50 bp configuration for all stage AFT samples and a 2×150 bp paired end configuration for all St16 samples, with a depth of ∼30-50 million reads. (GeneWiz). FASTQ files were pseudoaligned to the Drosophila melanogaster coding sequences (CDS, BDGP6, ensemble.org), abundances of transcripts were determined using Kallisto (release 0.43.1) as transcripts per million (TPM) units (Bray et al., 2016). Genes of interest were identified based on TPMs being different by at least 2FCs and overlapping with the GO-term gene list GO0031175 neuron-projection-development for Dmel (amigo.geneontology.org/amigo/term/GO:0031175). Genes are presented in heatmaps and Venn diagrams. The RNA-seq files have been deposited in GEO under accession number GSE113260.

Zeiss LSM700 or LSM800 confocal microscopes were used for fluorescent images; confocal stacks were merged using LSM software or ImageJ FIJI. Statistic calculations were performed in Graphpad prism software (v4.03). Cell counts and reporter (GFP) measurements were carried out using ImageJ FIJI and numbers were transferred to Graphpad prism. To address statistical significance, Student's t-test was used. Images and graphs were compiled in Adobe Illustrator. The numerical data underlying the graphs in the figures are included in Table S3.

We are grateful to Ward Odenwald, Alexander Kuzin, James B. Skeath, Alain Vincent, Krzysztof Jagla, Dimi Bieli, Markus Affolter, Takashi Suzuki, Konrad Basler, Johannes Bischof, the Drosophila Genomics Resource Center, the Developmental Studies Hybridoma Bank at the University of Iowa, the Bloomington Stock Center and the FlyORF stock center for sharing antibodies, fly lines and DNAs. We thank Doug W. Allan, John B. Thomas and Chris Q. Doe for critically reading the manuscript. We thank Erik Hultin for his help on this project during the course of his Bachelor thesis project. Behzad Salmani Yaghmaeian, Carolin Jonsson and Annika Starkenberg provided excellent technical assistance.