Functional and evolutionary insights from the Ciona notochord transcriptome

The notochord is the eponymous feature of the chordates and is an essential organ in chordate development. Despite its importance in the innovation of the chordate lineage, the evolutionary origins of the notochord are unclear. The details of notochord morphology and morphogenesis vary somewhat between the cephalochordates, tunicates and vertebrates, but in all chordates the notochord forms from dorsal midline mesoderm cells that express the transcription factors Brachyury and FoxA, and undergo striking convergence and extension movements (Flood et al., 1969; Corbo et al., 1997a; Holland et al., 2004; Friedman and Kaestner, 2006; Satoh et al., 2012). These movements transform a relatively isodiametric primordium into a characteristic long thin rod that can have both structural and/or signaling roles depending on the species (Stemple, 2005; Corallo et al., 2015). This rod becomes encapsulated in a perinotochordal basement membrane (PBM) as it forms, and the notochord cells typically form vacuoles or extracellular lumen pockets that help to inflate the PBM sheath to act as a hydrostatic skeleton. This structural role is particularly evident in species with swimming larval stages. Birds and mammals have smaller less-prominent notochords, but the notochord serves as a source of important midline patterning cues that control the dorsoventral patterning of the neural tube and the mediolateral patterning of the somite in all vertebrate embryos.

Notochord morphogenesis has been most studied in Xenopus, zebrafish and the ascidian Ciona. Non-canonical Wnt/PCP signaling has been implicated in mediolateral intercalation or related cell behaviors associated with convergent extension in all three species (Heisenberg et al., 2000; Wallingford et al., 2000; Jiang et al., 2005; Skoglund and Keller, 2010). Laminin, fibronectin and other extracellular matrix (ECM) components are known to be required for proper notochord morphogenesis in diverse chordates (Skoglund and Keller, 2007; Gansner et al., 2008; Veeman et al., 2008; Segade et al., 2016).

Ciona provides a particularly compelling model for studying notochord morphogenesis (Satoh et al., 2003; Passamaneck and Di Gregorio, 2005). The Ciona notochord consists of only 40 cells that develop through fixed and relatively well-understood cell lineages. The Ciona intestinalis genome is extremely small at ∼160 Mb and it diverged before the major genome duplications in the vertebrates (Dehal et al., 2002; Satoh, 2013). Modern molecular phylogenies show that it is the tunicates, including Ciona, and not the cephalochordates that are the closest outgroup of the vertebrates (Delsuc et al., 2006; Vienne and Pontarotti, 2006). Ciona was the second animal genome to be sequenced, and there are excellent community resources for genomics (Tassy et al., 2010; Brozovic et al., 2016).

Tissue-specific gene expression in the Ciona notochord has been of longstanding interest in the context of understanding cis-regulatory transcriptional control (Corbo et al., 1997b; Hudson and Yasuo, 2006; Hotta et al., 2008; Oda-Ishii et al., 2010; Katikala et al., 2013; Di Gregorio, 2017), the mechanisms of morphogenesis (Munro and Odell, 2002; Jiang et al., 2005; Hotta et al., 2007b; Veeman et al., 2008; Deng et al., 2013) and the enigmatic evolutionary origins of this defining innovation of the chordates (Satoh et al., 2012). Several previous studies have led to the systematic identification of Ciona notochord genes. Many of these came from a groundbreaking effort that used subtractive hybridization to identify transcripts upregulated in response to misexpression of the key notochord transcription factor Brachyury (Hotta et al., 1999; Takahashi et al., 1999). Others came from high-throughput in situ hybridization projects and microarray experiments (Satou et al., 2001b; Imai et al., 2004; Kugler et al., 2008; José-Edwards et al., 2011, 2013). Together, these projects identified 106 genes (summarized in Table S1) with either notochord-specific or notochord-enriched expression during intercalation and elongation. These genes have been a valuable resource for diverse functional studies, but we hypothesized that this was a far from complete list.

FACS sorting of fluorescently labeled cells followed by microarray analysis was originally used to identify genes expressed in the border cells of Drosophila ovaries (Borghese et al., 2006) and has since become a popular strategy for tissue-specific transcriptional profiling. FACS/microarray approaches have been used in Ciona to identify genes expressed in heart precursors (Christiaen et al., 2008) and in pigment cells (Racioppi et al., 2014), and also to identify seven notochord-expressed transcription factors (José-Edwards et al., 2011). Here, we take advantage of the greater sensitivity and dynamic range of RNAseq (Malone and Oliver, 2011), combined with flow-sorted notochord-enriched and notochord-depleted cell populations from dissociated embryos, to gain a comprehensive and unbiased view of transcript-level notochord gene expression across the entire genome. Ciona notochord cells were fluorescently labeled by electroporation of fertilized eggs with the reporter construct Bra>GFP, which strongly expresses GFP in the notochord, although there can also be weak expression in the mesenchyme. Three timepoints were used: early in intercalation (late neurula, Hotta stage 16) (Hotta et al., 2007a); late in intercalation (early tailbud, Hotta Stage 19.5); and well after the end of intercalation (late tailbud, Hotta stage 23) (Fig. 1A-C). At each time point, embryos were dissociated to single cells in Ca²⁺/Mg²⁺-free artificial seawater with trypsin, and then flow sorted to isolate populations of notochord-enriched GFP⁺ cells and notochord-depleted GFP⁻ cells (with the caveat that the notochord depletion is likely to be modest because of mosaic transgene expression). Triple biological replicates of purified RNA were analyzed by RNAseq, and both normalized read counts and differential expression were calculated using the Cufflinks software suite (Trapnell et al., 2012). This provides quantitative, genome-wide information on normalized transcript abundance in sorted notochord cells and on the fold enrichment in notochord compared with the rest of the embryo. Our analysis here focuses on the genes enriched in the notochord, but the data can also be mined to identify genes that are expressed in notochord without necessarily being enriched in notochord. The stage 23 timepoint later proved to be partially contaminated with muscle-specific transcripts (described below), so we have removed it from most analyses.

We identified 891 genes with higher expression levels in the notochord-enriched population compared with the notochord-depleted cell population at stage 16 (Fig. 1D) and 933 upregulated genes at stage 19.5 (Fig. 1E). When data from these two stages were combined, a total of 1364 notochord-enriched genes were identified; 460 of which were enriched at both stages (Fig. 1F and Table S2). Notochord-enriched genes spanned a wide range of expression levels, with normalized read count (fragments per kilobase mapped, FPKM) values ranging from 1 to 4515. They represent 11.8% of the 11557 genes expressed in the notochord-enriched sample at one or both time points above a threshold of 1 FPKM. Our sequencing depth and biological replication were sufficient to detect statistical significance (corrected for multiple comparison) for genes that were only enriched in notochord by 1.47-fold. As seen in Fig. 1F,G, a smooth distribution of fold enrichment values was seen between these modestly enriched genes and the most enriched genes. By comparison, the 20 most enriched genes at stage 16 and at stage 19.5 showed mean enrichments of 49- and 23-fold, respectively. This group includes the key notochord transcription factors Brachyury and FoxAa, as well as previously identified Ciona notochord genes Noto4, Noto8, FGL2, TGFbeta2/3, SLC26Aα and THBS3/4 (thrombospondin) (Hotta et al., 1999, 2000; Imai et al., 2004; Katikala et al., 2013). The top hits also show considerable overlap, with half of the ten most enriched genes at stage 16 also in the top ten genes for stage 19.5 (Table 1).

As an initial validation, we examined the RNAseq expression data for the 106 genes previously identified as having notochord-specific or notochord-enriched expression (highlighted in green in Fig. 1D,E). Our genome-wide assay for notochord-enriched expression identified 81 of these 106 known notochord genes. Of the 23 that we missed, previously published in situ results showed that 16 were qualitatively expressed more strongly in another tissue compared with notochord (an expression pattern that we would expect to miss in our differential expression analysis) and five have conflicting published in situ hybridization data. This provided an initial confirmation that our RNAseq approach was likely to have few false negatives for genes with highly specific notochord expression.

Compared with previous screens, our RNAseq approach identifies a much greater range of expression levels, with less bias for highly expressed genes. At stage 19.5, for example, the mean notochord expression of the genes identified as notochord enriched by RNAseq was 130 FPKM, whereas it was 504 FPKM for the 106 previously characterized notochord genes and 697 FPKM for the 40 of those identified in the original subtractive hybridization screen (Takahashi et al., 1999). Many of the newly identified putative notochord genes were nevertheless highly expressed, with 14 of the 20 most highly expressed genes at stage 19.5 not found by the previous screens.

To further validate our RNAseq data, we selected 111 genes enriched in notochord at stage 16 and/or stage 19.5 and tested them by in situ hybridization at stages 16, 19.5 and 23. To better understand the data, we selected genes that spanned a broad range of notochord expression levels (from 3.5 to 4515 FPKM) and enrichment values (from 1.5- to 90.9-fold) (Fig. S1). We included genes enriched at both single stages and at multiple stages. Predicted gene functions were not a criteria in selecting genes for in situ hybridization, with the exception of the 16 transcription factors discussed later. Of the 98 that were detectable by in situ hybridization, 88 (90%, Agresti-Coull 95% confidence interval: 82-95%) were expressed in notochord. These were sometimes expressed in other tissues as well, but notochord staining was distinct and usually predominant at the predicted stages. A representative sample of these genes is shown in Fig. 2A, and all results are summarized in Table S2. We were able to visualize notochord-specific expression of some genes, such as ARF1/3/6 (5 FPKM and 8.7-fold enriched at stage 19.5, Fig. 2A), that were expressed at very low levels but highly enriched in the notochord. Genes with more modest notochord enrichment scores (enriched twofold or less) needed to have higher notochord expression levels (generally greater than 100 FPKM) to be detectable as enriched. PTK2 is shown as an example in Fig. 2A. Both the timing of expression and the relative enrichment by in situ hybridization were well matched to the RNAseq data. Eight of the remaining ten genes were expressed strongly in mesenchyme cells, and two in muscle cells (examples in Fig. 2B,C). This modest contamination with mesenchymally expressed genes likely reflects the low levels of GFP expression induced in this tissue by the Bra>GFP construct. The stage 23 RNAseq data largely matched the in situ results for genes that first became enriched in notochord at stage 16 or stage 19.5.

We also tested 41 (36 detectable) genes that did not become enriched in notochord by RNAseq until stage 23, and found that these were largely expressed in mesenchyme and/or muscle (Fig. 2D,E). Only 11% (four genes) were notochord expressed. It is not clear whether this reflects failure to fully disassociate older embryos, muscle cell autofluorescence or suboptimal gating choices while FACS sorting. We therefore excluded the stage 23 RNAseq data from our subsequent analyses.

For stages 16 and 19.5, however, our 90% in situ validation rate indicates that the RNAseq data at these stages has a very low fraction of false positives. By comparison, the previous effort to systematically identify notochord genes by ectopic expression of Brachyury coupled with subtractive hybridization had an in situ validation rate of less than 10% and only identified 40 genes in total (Hotta et al., 1999; Takahashi et al., 1999). Microarray screening for targets of the notochord-expressed transcription factor TBX2/3 had a moderately better validation rate of 29%, with 20/70 detectable genes confirmed as notochord expressed (José-Edwards et al., 2013). Our extensive validation predicts that about 90% of the 1364 notochord genes we have identified by RNAseq are likely to be bona fide. Our results support the tissue-specific transgene/FACS/RNAseq (FACSseq) approach for the systematic and relatively unbiased quantification of tissue-specific gene expression.

We used topGo (Alexa et al., 2006) to identify Gene Ontology (GO) terms enriched in our notochord-expressed gene sets. Gene model to GO code mappings were downloaded from the ANISEED ascidian community database (Brozovic et al., 2016) and are all inferred from sequence conservation and/or orthology with vertebrates. We ran the analyses for the group of genes that are notochord enriched at both stage 16 and stage 19.5 (460 genes), as well as for the 891 and 933 genes that are enriched at stage 16 or at stage 19.5 (Table 2 and Table S3).

For the genes enriched at both stages, the GO code analysis indicated a massive enrichment of genes associated with the ECM and cell-ECM interactions. This includes both commonplace ECM components, such as laminins, fibronectins, nidogen, perlecan and collagens, as well as many more exotic ECM components, such as fibulin, matrilin, hemicentin, FREM1, tenascin, fibrillin, papilin and thrombospondin. Genes involved in ECM modification were also common, including prolyl hydroxylase, lysyl hydroxylase, peroxidasins, heparanase and serpins. Components of the focal adhesion complex linking the cytoskeleton and ECM were also enriched, including integrinA3/6/7, tensin, actinin, talin, parvin and fermitin. These are all consistent with the crucial role of the perinotochordal and intranotochordal ECM in notochord morphogenesis, both during intercalation and after (Skoglund and Keller, 2007; Gansner et al., 2008; Veeman et al., 2008; Segade et al., 2016).

Many of the other GO terms enriched with respect to the genes enriched at stage 16 and 19.5 were initially perplexing (neuromuscular junction, perception of sound, axon guidance, PDGF binding) but closer inspection indicated that these all reflected the enrichment of ECM-related genes. Although significant expression of ECM components was predicted given the known importance of the perinotochordal basement membrane, our data show that the most defining signature of the Ciona notochord transcriptome is the expression of a diverse set of ECM structural components, modifying enzymes and ECM-cytoskeleton linkers. We have identified many ECM components not previously known to be notochord enriched that will be useful in future studies of notochord morphogenesis. This is particularly important given the potential for functional redundancy in the ECM.

The GO codes enriched with respect to all genes enriched in the stage 16 notochord also include several associated with transcription, RNA processing, translation and the COPII ER to Golgi transport complex. These potentially reflect a general upregulation of the basic machinery for producing and secreting proteins, likely in preparation for the secretion of the perinotochordal ECM. Genes whose products are associated with GTPase activator activity (including both regulators of heterotrimeric G proteins and activators of small GTPases) were also over-represented in the notochord at this stage, and provide exciting candidates for regulatory roles in early notochord morphogenesis.

In the set of genes enriched at notochord at stage 19.5, we again found over-represented GO terms related to the ECM and cell-ECM binding. We also found over-represented GO terms related to both cell-cell adhesion and the actin cytoskeleton. In particular, we note that virtually the entire ARP2/3 complex (seven out of nine components) is enriched in the Ciona notochord at stage 19.5 (Table 2). This complex plays a key role in nucleating the growth of branched actin filaments (Swaney and Li, 2016), and has not been previously implicated in notochord morphogenesis. Stage 19.5 is quite late in the process of mediolateral intercalation, so we speculate that the ARP2/3 complex is involved in the subsequent morphogenetic event where the notochord becomes longer via the cells slowly changing from disk-shaped to drum-shaped, which is known to be driven by actomyosin contractility in the cell cortex (Dong et al., 2011).

In addition to our unbiased assessment of enriched GO terms, we also searched the list of 1364 notochord-enriched genes for predicted functions in cell signaling, channel/transporter activity and transcriptional regulation. Signaling genes included numerous components of the noncanonical Wnt/PCP pathway, which has been previously implicated in Ciona notochord morphogenesis (Jiang et al., 2005). Enriched noncanonical Wnt/PCP components included Prickle, Strabismus, Daam1/2 and ROR1/2. There were also several Wnt ligands and receptors that could potentially be involved in canonical or noncanonical Wnt signaling, including Wnt5b, Wnt9, Frizzled5/8 and Frizzled4/9/10. Surprisingly, APC and Casein kinase Iε were also enriched, both of which are thought to be specifically involved in canonical Wnt/β-catenin signaling (Hart et al., 1998; Peters et al., 1999; Angers and Moon, 2009).

Channels and transporters are of great interest with respect to the inflation of the notochord lumen and other aspects of notochord cell behavior. Even though these did not prove to be enriched GO codes, we found that many channels and transporters were enriched in the notochord (see Table S2). These included several predicted potassium channels (KCNB1/2, KCNH3/4/8, KCNJ4/12/8, KCNN1/2/3 and KCNQ2/3/4/5), 25 different members of the Solute Carrier (SLC) superfamily with diverse predicted specificities (Na⁺, K⁺, Cl⁻, Zn⁺², sulphate, anions and cations) and three TMEM family channels. The anion transporter SLC26Aα was previously identified as having notochord-enriched expression as well as playing an important role in notochord lumen expansion (Hotta et al., 1999; Deng et al., 2013). Deng et al. predicted that additional solute transporters and channels would also function during notochord tubulogenesis, and our data provide a list of candidates for future studies.

The notochord-enriched gene set also contained over 100 genes with potential roles in transcriptional regulation, including transcriptional activators and transcriptional repressors, as well as general transcriptional regulators such as chromatin remodeling factors, histone modifying enzymes, Mediator complex components and elongation factors. The transcriptional activators included several previously identified as expressed in the Ciona notochord, such as Bra, FoxAa, TBX2/3, Sall-a, STAT5/6, NFAT5 and Klf15 (Corbo et al., 1997,b; Imai et al., 2004; José-Edwards et al., 2011), as well as many potential new notochord TFs. Some of these can be excluded based on published mesenchyme expression patterns, but we have performed in situ hybridization on 16 transcription factors with ambiguous or no published expression patterns at the relevant stages and found that 11 have subtle but bona fide notochord expression (the rest were expressed in the mesenchyme) (Fig. 3 and Table S2). Additional in situ validation will be needed to define the full set of notochord-enriched candidate transcriptional regulators.

In addition to expression enrichment, the Cufflinks suite can test for differences in promoter use and differential splicing (Trapnell et al., 2013). At stage 16, 124 genes showed evidence of differential transcriptional start site usage in notochord-enriched versus -depleted cell pools (3.6% of the 3492 genes with more than one primary transcript annotated and expressed). Two-hundred and thirty-nine primary transcripts (3.15% of 7589 primary transcripts with sufficient alignments for testing) showed evidence of differential splicing between the two cell populations at this stage. Full results for both stages are shown in Tables S4 and S5. This is the first study to systematically estimate tissue-specific promoter use and alternative splicing in Ciona, so it is not clear at present if these rates are distinctive. One group previously estimated that about 15-20% of all genes in Ciona can be alternatively spliced, but did not look at individual tissues (Kim et al., 2007).

Transcription factor overexpression followed by transcriptional profiling provides an alternate strategy for identifying genes with tissue-specific expression. We identified many more notochord-enriched genes by FACSseq than were previously identified by Brachyury overexpression (Hotta et al., 1999; Takahashi et al., 1999), but it was unclear to what extent this reflected the more quantitative and genome-wide nature of RNAseq compared with subtractive hybridization. Alternatively, it might reflect a more fundamental difference between the genes induced by Brachyury misexpression and the normal notochord transcriptional profile. To test this, we broadly misexpressed Brachyury under the control of the FoxAa enhancer in a similar fashion to Takahashi et al. (1999) (Fig. 4A,B), but then used RNAseq to identify upregulated genes as compared to GFP-expressing control embryos. Given that Brachyury is thought to be at the top of the notochord transcriptional cascade, we hypothesized that most of the 1364 notochord genes identified by FACSseq would be upregulated to some extent by ectopic Brachyury expression.

We identified 925 genes with increased expression in FoxAa>Bra-expressing embryos when compared with Bra>GFP-expressing controls (Table S6). There was unexpectedly little overlap, however, between these genes and the notochord FACSseq gene set, with only 264 found in both (Fig. 4E). Genome-wide there is a positive but weak correlation between notochord specificity (FACSseq fold-change) and upregulation in response to Brachyury misexpression (FoxAa>Bra fold-change) (Fig. 4C), but the majority of notochord-enriched genes show no statistically significant response to ectopic Brachyury expression (Fig. 4D). The sequencing depth and triple biological replicates used showed statistically significant differential gene expression for genes with a log2 fold-change as low as ∼0.4, so this was a very sensitive test.

The FoxAa cis-regulatory region used drives expression in neural, endodermal and mesenchymal lineages in addition to notochord. It has previously been suggested that FoxAa>Bra expression transforms these other cells to notochord fate (Takahashi et al., 1999). Our RNAseq data, however, shows that only a subset of notochord-enriched genes are upregulated. To better explore the transcriptional response, we defined panels of tissue-specific markers by screening the ANISEED gene expression database for genes with distinct tissue-specific expression in the main Ciona tissues at the relevant stages. We found that markers of muscle and epidermal fate were unaffected by Brachyury misexpression. Many markers of neural, endodermal and mesenchymal fate were downregulated, but others were unaffected and a few were paradoxically upregulated (Fig. 4F and Table S7). This again supports the hypothesis that Brachyury overexpression only partially reprograms other cells to notochord fate.

Of the 661 genes transcriptionally upregulated by Brachyury misexpression but not notochord enriched by FACSseq, many were either not expressed or only minimally expressed in control embryos (311 at FPKM<2, 428 at FPKM<5). We searched the ANISEED gene expression database and found that 18% of the 661 had previously reported in situ hybridization results at comparable stages. Of these 120 genes, only 14% had expression in notochord (typically in combination with some other tissue), 13% in mesenchyme, 8% in muscle, 13% in endoderm, 35% in neural cells and 26% in epidermis. Forty percent showed either no expression or diffuse/nonspecific expression. We performed in situ hybridization for 18 genes that had not been previously characterized and found that half of the genes had widespread expression (Fig. S2). Only three genes showed specific expression in the notochord, again in combination with other tissues. Most of the 18 genes had their strongest expression in muscle and/or mesenchyme. Together, these results show that only a small fraction of the genes that are upregulated by Brachyury overexpression but not notochord enriched by FACSseq are normally expressed in the notochord at all.

One possible explanation for why Brachyury misexpression only induces ∼20% of notochord-enriched genes is that morphogenesis is severely perturbed in FoxAa>Bra embryos and the notochord GRN might perhaps be dependent on the normal progression of morphogenesis. To test this, we blocked morphogenesis entirely after notochord induction by treatment with cytochalasin B and then used RNAseq to look for differential gene expression as compared with control embryos. Genome-wide we found only 43 genes, showing significant downregulation at stage 19.5, and only four of these were notochord enriched by FACSseq (Table S8). Feedback from morphogenesis thus does not appear to be a major factor in the notochord GRN, at least not over the timeframe examined.

We also compared the notochord FACSseq and FoxAa>Bra RNAseq data with the Brachyury ChIP-chip data of Kubo et al., who used chromatin immunoprecipitation of a tagged Brachyury construct coupled with a sparse but genome-wide microarray to determine Brachyury occupancy at the early gastrula stage (soon after notochord induction) (Kubo et al., 2010). They identified 2252 genes with evidence of Brachyury occupancy, 464 of which are notochord enriched according to our FACSseq data. Two-hundred and eighty-six of the genes that are FoxAa>Bra induced but not notochord enriched are predicted to be Brachyury bound, but 375 are not (Fig. 4E and Table S9). This supports a complex model of notochord gene regulation in which many notochord genes are not direct Brachyury targets and most Brachyury-bound regulatory regions do not drive notochord-specific expression.

There are both profound similarities and also distinct differences in notochord anatomy and morphogenesis between the cephalochordates, tunicates and vertebrates. To define a core set of genes with shared notochord expression, we compared our set of Ciona notochord-enriched genes to the gene expression patterns cataloged in the mouse and zebrafish gene expression databases (Howe et al., 2013; Finger et al., 2017). We chose these species because they have the most comprehensive gene expression databases, but the mouse and zebrafish notochords are also quite different in their anatomy and development (Myers et al., 2002; Stemple, 2005; Yamanaka et al., 2007), and thus serve as diverse vertebrate comparison groups. For both species, we queried the databases for all genes expressed in the notochord during roughly comparable stages of notochord morphogenesis and then identified the Ciona orthologs (Brozovic et al., 2016).

We found many genes that are expressed in the notochord in both Ciona and in fish and/or mouse embryos. Of the 132 mouse notochord genes with orthologs in Ciona, 31 were present in our dataset (Table S10). This was significantly more overlap than the 14 genes that would be expected by chance (P=2.1×10⁻⁵ by the hypergeometric distribution test). Similarly, the orthologs of 99/359 fish notochord genes were enriched in our sorted notochord cells (Table S11). This is again greater that the 39 genes expected by chance (P=2.6×10⁻¹⁹). Overall, there were 12 genes identified as notochord expressed in all three species (Table 3).

Although we have identified more than twice as many genes with notochord expression shared between Ciona and vertebrates than would be predicted by chance alone, it is likely that the true overlap is considerably greater given the varying techniques used (RNAseq versus microarrays versus in situ hybridization; genome-wide studies versus gene-by-gene studies, etc.), the varying coverage in high-throughput in situ screens and the precise developmental stages examined. For example, many of the zebrafish notochord genes were identified in high-throughput in situ hybridization screens (Thisse et al., 2004), whereas the most common source of the mouse notochord genes was an in situ validated microarray experiment using flow-sorted Noto>eGFP-expressing cells cross-matched to FoxA2 ChIP data from adult mouse liver cells (Tamplin et al., 2011).

When compared with our full set of 1364 Ciona notochord-enriched genes, the 112 that are also notochord expressed in mice and/or fish have a similar distribution of notochord enrichment values, but higher absolute expression levels. The mean FPKM at stage 19.5 of genes with conserved notochord expression was 2.37-fold higher than for all genes enriched at that stage (D=0.388, P=1.421×10⁻¹¹ by the Kolmogorov–Smirnov test). This largely reflected the absence of genes with low notochord FPKMs in the conserved gene set, which is potentially of evolutionary significance but more likely reflects the more-sensitive RNAseq assay in which the Ciona notochord genes were here identified.

The 12 core genes shared by all three species fell into four distinct classes. Brachyury (Satoh et al., 2012) and Sall-b (also called SAL1/3/4) (Ott et al., 1996; Thisse et al., 2001) (Fig. 3) are both notochord transcription factors. Another key notochord transcription factor, FoxAa, is also expressed in notochord in all three species, but was absent from our high-throughput analysis because the mouse and Ciona genes are somewhat divergent and were not matched in the genome-wide ortholog database in ANISEED, despite their orthology having been previously confirmed (Yagi et al., 2003). Collagen2A1, fibronectin1, prolyl hydroxylase (also called LEPREL1/2) and PAPS Synthetase all encode either structural components of the ECM or proteins involved in ECM post-translational modifications (Faiyaz ul Haque et al., 1998; Hynes and Naba, 2016). Prickle, Daam1/2 and ERM (the single Ciona ortholog of ezrin/radixin/moesin) all encode proteins involved in controlling the actin cytoskeleton, either directly or via roles in noncanonical Wnt/PCP signaling (Tsukita and Yonemura, 1999; Skoglund and Keller, 2010). Prickle and Daam1/2 have both been implicated in mediolateral intercalation (Habas et al., 2001; Takeuchi et al., 2003; Veeman et al., 2003), whereas a Ciona ERM morpholino had a later defect in the disk-to-drum notochord cell shape change that follows intercalation (Hotta et al., 2007b). The fourth class is an ad hoc grouping of LRIG, CRIM1 and NAV3. All three have previously been implicated in aspects of cell signaling and/or the cytoskeleton, but have not been extensively characterized (Wilkinson et al., 2003; Kinna et al., 2006; van Haren et al., 2014; Simion et al., 2014). Their conserved expression across broad phylogenetic distances makes them attractive candidates for future functional studies of their roles in notochord morphogenesis.

The evolutionary origins of the notochord have long been debated. It is unique to the chordates and there are no compelling transitional fossils. The hemichordate stomochord is superficially reminiscent, but lacks Brachyury expression and a notochord-like basement membrane, and is now thought to be derived from pharyngeal endoderm (Peterson et al., 1999; Satoh et al., 2014). It has recently been proposed that a structure homologous to the notochord can be found in the ventral midline mesoderm of protostomes (Lauri et al., 2014; Brunet et al., 2015). This ‘axochord’ was characterized in the annelid Platynereis by virtue of its conserved expression of a set of transcription factors and other genes that are expressed in the notochord of many chordates. Not all were specific to the axochord, but the overlap of their complex expression patterns was relatively specific to this structure.

Given the considerable support for the dorsoventral inversion hypothesis at the level of gene expression (Arendt and Nübler-Jung, 1997), it is in many respects not surprising that ventral midline mesoderm in a protostome might have a conserved transcriptional signature with respect to dorsal midline mesoderm in chordates. It has been suggested, however, that the axochord actually engages in notochord-like behaviors such as mediolateral intercalation (Lauri et al., 2014), and that the notochord might in some respects thus have surprisingly deep origins before the divergence of protostomes and deuterostomes. We would argue, however, that many tissues engage in intercalatory behaviors independently of the noncanonical Wnt/PCP-dependent tissue rearrangements that result in the massive convergence and extension of the chordate notochord (Zallen and Wieschaus, 2004; King et al., 2009; Walck-Shannon and Hardin, 2014). It is thus of great importance to know whether the axochord also expresses a broad range of notochord effector genes.

Ciona lacks notochord-specific expression of several genes, including Noggin, Twist, SoxD, SoxE, Hedgehog (Hh) and Slit, that were used as key markers of notochord identity by Lauri et al. (2014) (Table S12). Many of these are expressed in the notochord in Amphioxus as well as vertebrates, likely reflecting lineage-specific loss of notochord expression in the tunicates (Lauri et al., 2014). These genes were thus likely expressed in notochord in the last common ancestor of the chordates, but it is important to note that they do not all appear to be key effector genes for notochord morphogenesis per se. Noggin and chordin, for example, are involved in the function of the axial mesoderm as the Spemann/Mangold dorsal organizer (Smith and Harland, 1992; Smith et al., 1993; Sasai et al., 1994). Netrin, Slit and Hh are involved in later aspects of neural patterning and axon pathfinding (Echelard et al., 1993; Kennedy et al., 1994; Roelink et al., 1994; Yuan et al., 1999). Tunicates have lost the dorsal organizer as a key mediator of DV patterning (Lemaire et al., 2008) and appear to make less use of the notochord as a source of midline patterning cues (Yamada et al., 2009), but they nevertheless make a stereotypically chordate notochord that undergoes all of the expected morphogenetic behaviors. To better understand the evolutionary origins of the notochord, it will thus be important to determine the protostome expression patterns for a broad range of notochord-enriched effector genes conserved across the chordates, with a particular emphasis on genes likely to play essential conserved roles in notochord morphogenesis. Key questions suggested by our data are whether the protostome axochord is similarly enriched for Prickle or other PCP components, and whether or not it has a similar signature of ECM components beyond collagen.

Ciona notochord morphogenesis involves a broad range of morphogenetic behaviors, including asymmetric division (Veeman and Smith, 2013), mediolateral intercalation (Munro and Odell, 2002; Jiang et al., 2005), boundary capture (Veeman et al., 2008), controlled changes in aspect ratio (Dong et al., 2011) and tubulogenesis (Deng et al., 2013; Denker et al., 2013). This richness of cell behavior can be studied in toto in the 40 Ciona notochord cells, making it an attractive model for integrative studies of how tissue-specific gene expression controls tissue-specific morphogenetic behaviors. Particularly compelling candidate genes for loss-of-function studies include the large set of ECM structural components and modifying enzymes, the ARP2/3 complex components, and the various channels and transporters. LRIG, CRIM1 and NAV3 are also of particular interest given their conserved notochord expression in mice and fish.

Beginning with the early effort of Takahashi and colleagues to systematically identify notochord-specific genes by Brachyury overexpression and subtractive hybridization (Hotta et al., 1999; Takahashi et al., 1999), the Ciona notochord has become an important model for studying the mechanisms of tissue-specific gene expression (Passamaneck et al., 2009; José-Edwards et al., 2013, 2015; Katikala et al., 2013; Farley et al., 2016). The 1364 genes with notochord-enriched expression identified here represent a more than 10-fold increase in known notochord genes, and provide a comprehensive, genome-wide view of gene expression in the notochord during key stages of mediolateral intercalation. The Ciona notochord transcriptome provides a foundation for systematically dissecting the notochord effector gene regulatory network and systems-level studies of tissue-specific enhancer function.

Although used inconsistently in the literature, master regulator genes are typically defined as genes at the top of tissue-specific transcriptional cascades that are both necessary for initiating tissue-specific transcriptional programs and also capable of reprogramming other cell types to that fate when ectopically expressed (Chan and Kyba, 2013), as has been demonstrated with well known examples such as MyoD (Tapscott et al., 1988) and Drosophila eyeless (Halder et al., 1995). Brachyury has long been treated as a de facto master regulator gene for the Ciona notochord, in the sense of being both necessary for notochord fate but also sufficient to induce notochord fate when expressed in other cell types. Its necessity is well established, as various loss-of-function assays, including a likely null mutation show that it is required for notochord development (Satou et al., 2001a; Yamada et al., 2003; Chiba et al., 2009). If truly sufficient for notochord fate, however, misexpressing it outside the notochord would recapitulate the normal notochord transcriptional program in these other cell types. This clearly happens on a certain level, as Brachyury misexpression induces expression of many notochord-specific markers (Hotta et al., 1999; Takahashi et al., 1999). These arguments are confounded, however, by the fact that most previously known notochord markers were identified specifically by screening for transcriptional targets of Brachyury misexpression.

Here, we have identified a large and unbiased set of notochord-enriched genes and find that relatively few of them are upregulated in response to ectopic Brachyury (Fig. 4E) FoxAa>Bra-expressing cells instead appear to take on a hybrid character in which some but not all notochord markers are upregulated, some but not all markers of endodermal, mesenchymal and neural fate are downregulated, and many other genes are ectopically upregulated (Fig. 4F). This potentially reflects the combinatorial control of gene expression as influenced by Brachyury in combination with other transcriptional regulators expressed in FoxAa-expressing cells but not in wild-type notochord cells. The idea that Brachyury misexpression only partially transforms other cell types to notochord fate is supported by the observation that most FoxAa>Bra-expressing cells outside the notochord form an overtly undifferentiated mass and never undergo mediolateral intercalation or other characteristic notochord cell behaviors. Notochord-enriched gene expression appears almost entirely normal in cleavage-blocked embryos, so the failure of FoxAa>Bra-expressing cells to morphologically differentiate as notochord likely reflects these transcriptional differences and not simply defects in embryonic architecture. An intriguing question for the future is whether the observation of numerous genes that are notochord enriched but not upregulated by FoxAa>Bra is more reflective of factors expressed outside the notochord that repress the notochord transcriptional regime, or whether other transcription factors might act in parallel to Brachyury at the top of the notochord GRN. In support of the latter hypothesis, we note that, despite being a likely null allele, the Ciona brachyury mutant chobi has a relatively mild phenotype in which a distinct notochord forms but has defects in intercalation and elongation (Chiba et al., 2009). In addition, at least one notochord transcription factor is expressed normally in chobi mutant embryos (José-Edwards et al., 2011).

Ciona intestinalis type A were collected in San Diego (CA, USA) and shipped to Kansas State University by Marine Research and Educational Products. Although long identified as Ciona intestinalis, recent analysis has suggested that type A should potentially be considered a distinct species named C. robusta (Pennati et al., 2015). Adult Ciona were maintained in a recirculating aquarium. Standard fertilization, dechorionation and electroporation protocols were used (Veeman et al., 2011). For sorting experiments, 80 μg of Bra>GFP plasmid was electroporated.

For each of the three biological replicates we collected 2000 electroporated embryos and 500 unelectroporated control embryos at each timepoint. Embryo dissociation was modified from methods described previously (Christiaen et al., 2008). Embryos were washed twice with Ca²⁺/Mg²⁺-free artificial seawater [CMF-ASW; 449 mM NaCl, 9 mM KCl, 33 mM Na₂SO₄, 2.15 mM NaHCO₃, 10 mM Tris-HCl (pH 8.0), 2.5 mM EGTA] in BSA-treated eppendorfs, spinning at 730 g for 1 min to pellet. Embryos were resuspended in 600 μl CMF-ASW and 60 μl of 1% trypsin (0.1% final concentration). Samples were pipetted for 3 min to dissociate embryos, then 600 μl of ice-cold CMF-ASW+0.1% BSA was added to stop the reaction. Cells were spun for 3 min at 730 g, washed once with 1 ml ice-cold CMF-ASW+BSA and resuspended in 300 μl ice-cold CMF-ASW+BSA. Cells were filtered through a BSA-treated 70 μm filter into 200 μl ice-cold CMF-ASW+BSA. Final volume was 500 μl.

Cell sorting was performed using a Biorad S3 at the Flow Cytometry Lab in Kansas State University College of Veterinary Medicine. For each timepoint, an unelectroporated cell sample was first analyzed to determine the levels of background autofluorescence, and gates on forward scatter and GFP fluorescence were conservatively selected to enrich for single cells with or without strong GFP signal. The electroporated sample was then sorted directly into RNA lysis buffer (Zymogen Quick-RNA miniprep kit R1054). For GFP-positive notochord-enriched cells, 2000-15,000 cells were purified into 500 μl buffer, while 100,000 GFP-negative (notochord depleted) cells were simultaneously sorted into 2 ml of lysis buffer. Samples were stored at −80°C until RNA purification.

RNA was purified with the Zymogen Quick-RNA miniprep kit (R1054) and eluted in 50 μl elution buffer. The yield was 0.6 to 1.25 μg of RNA from notochord-depleted cells and 40-140 ng from notochord-enriched cells. RNA (10-60 ng per sample) was treated with DNase (Roche) and concentrated to 10 μl with the Zymogen Clean and Concentrator kit (R1015). The final purified RNA concentration was 0.5-4.0 ng/μl.

Amplified cDNA was prepared using the NuGEN Ovation RNA-Seq System V2 (7102). Purified RNA (2.5 μl) was used in half-reactions following their standard protocol. cDNA samples were fragmented with a Covaris S220 Ultrasonicator using the recommended settings to obtain an average fragment length of 300 bp and the NEBNext DNA Library Prep Master Mix set for Illumina from New England BioLabs (NEB) was used for library construction. The resulting libraries were quality checked with an Agilent 2100 Bioanalyzer and quantified by RT-PCR. Paired-end (2×100 bp) sequencing was performed on an Illumina HiSeq 2500 at the Kansas University Genome Sequencing Core with a read depth of 30-80 million reads per sample.

Fertilized dechorionated eggs were electroporated with either 40 μg of FoxAa>Venus-Bra plasmid or 40 μg of control plasmid (Bra<GFP) on three separate days. For each of the replicates, 400 experimental and 400 control embryos were collected at stage 19.5 into RNA lysis buffer and stored at −80°C until RNA purification. RNA was purified with the Zymogen Quick-RNA miniprep kit (R1054), genomic DNA was removed using the in-column DNase treatment step in the kit and purified RNA was eluted in 50 μl water, at a concentration of 52-60 ng/μl. Libraries were constructed with the TruSeq stranded mRNA library kit (Illumina) using standard protocols. Libraries were quality checked with an Agilent Tapestation and quantified by qPCR. Single-end (1×100) sequencing was performed on an Illumina HiSeq 2500 at the Kansas University Genome Sequencing Core with a read depth of 33.5 to 38.5 million reads per sample.

Dechorionated embryos were treated with 10 μM cytochalasinB (1.4 μl of 21 mM cytochalasin-DMSO stock in 3 ml of ASW) or with DMSO alone as a control. Drug was added at 3.75 hpf at 21.5C (64-cell stage). Embryo harvesting at stage 19.5, RNA purification, library preparation and sequencing were as described for FoxAa overexpression. Read depth was 8.5 to 11.6 million reads per sample.

Cufflinks (Trapnell et al., 2012, 2013) was used for expression quantitation. Reads were aligned to the Ciona KH2012 gene models (KH.KHGene.2012.gff3, retrieved from ANISEED, http://www.aniseed.cnrs.fr/) (Brozovic et al., 2016) with Tophat, abundances were calculated with Cuffquant, and differential expression, splicing and promoter usage were calculated with Cuffdiff using the default parameters. Statistical significance was reported as q-values adjusted for multiple comparisons using a false discovery rate of 0.05.

TopGO (Alexa and Rahnenfuhrer, 2016) was used to test for enrichment of GO terms in our notochord-expressed gene lists. GO annotations for the KH2012 gene models were downloaded from ANISEED. Enriched gene lists for each stage were tested against a limited universe of genes with an FPKM of at least 0.5 in the notochord-depleted and/or notochord-enriched cell populations at the appropriate stage. Enriched GO terms were identified using Fisher's exact test and the weight01 algorithm in topGo with a P-value cutoff of 5×10⁻⁴. This algorithm does not require adjustment for multiple comparisons.

Probe synthesis, embryo collection and in situ hybridization were performed as described (Reeves et al., 2014). Clones used as templates for probe synthesis were obtained from the Ciona Gateway whole-ORF gene collection (Satou et al., 2002), or amplified from cDNA and cloned into pBSII SK(−). The gene collection IDs or PCR primers are listed in Table S13.

Genes expressed in the fish ‘notochord’ or ‘axial chorda mesoderm’ from shield to 22- to 25-somite stages were retrieved from the Zebrafish Model Organism Database (ZFIN) (http://zfin.org/) in September 2016. Genes expressed in the mouse notochord from Theiler stage 12-21 were retrieved from the Gene Expression Database (GXD) (http://www.informatics.jax.org) in September 2016. KH2012 orthologs were downloaded from ANISEED and used to identify genes with conserved notochord expression.

The authors thank the University of Kansas Genome Sequencing Core, the KSU CVM Confocal Core, the KSU CVM Flow Cytometry Lab and the KSU Integrated Genomics Facility for technical support.