Cis-regulatory properties of medaka synexpression groups

Tissue determination during embryogenesis is governed by dynamic developmental gene expression, and failures of proper gene expression in time and space lead to embryonic defects. Spatiotemporal gene expression is tightly controlled within a developmental gene regulatory network (GRN) (Peter and Davidson, 2009). A major challenge in the genomic era is to decipher this network by identifying its components and the relationships between them.

Subsets of developmental genes called synexpression groups (SGs) are characterised by strikingly similar spatiotemporal gene expression patterns throughout embryo development. They have been proposed to be involved in common biological pathways, highlighting their central role as basic modules of developmental GRNs, promoting the coordinate expression of gene batteries and potentiating rapid evolutionary changes (Niehrs and Pollet, 1999). To identify SGs, a number of large-scale whole-mount gene expression screens have been performed in Drosophila (Tomancak et al., 2002), Xenopus (Gawantka et al., 1998) and mouse (Visel et al., 2007). These studies have highlighted the usefulness of SGs in predicting biological functions for uncharacterised genes, based on the assumption that genes belonging to the same SG are involved in a common biological process (Amaya, 2005).

Niehrs and Pollet (Niehrs and Pollet, 1999) have proposed that genes within an SG display a similar expression pattern because they possess common cis-regulatory elements that are bound by a common regulatory trans-acting factor (t-AF). Thus, genes belonging to the same SG should share common sequence motifs that are reflective of the t-AF DNA binding site. Support for this hypothesis has been provided for several SGs in amniotes (Ernsberger, 2000), Ciona (Brown et al., 2007), Xenopus (Karaulanov et al., 2004), medaka (Ramialison et al., 2008), mouse (Visel et al., 2007) and human (Grade et al., 2009).

Although these studies have been performed in specific SGs, they suggest that their shared cis-regulatory input leads to co-expression. In this pilot analysis we present a systematic investigation of the functional and regulatory properties of selected SGs in the developing medaka fish embryo with a view to uncovering general principles of SG regulation.

First, we identified SGs by collecting similar spatiotemporal gene expression patterns from the Medaka Expression Pattern Database (MEPD) (Henrich et al., 2003; Henrich et al., 2005) and from a new random in situ hybridisation (ISH) screen. From these data, we built a network of co-expressed genes at three developmental time points. We used data reduction to collapse highly similar gene expression patterns into expression ‘mountains’, followed by a higher-order hierarchical clustering of these mountains into ‘mountain ranges’. From this co-expression network, SGs were identified by systematically searching for modules of co-expressed genes that share a common biological process. By looking for over-represented DNA motifs de novo in each SG, we revealed that SGs share common complex putative regulatory input consisting of different cis-regulatory modules (CRMs) that are potentially able to drive synexpression. We performed an in-depth analysis of an SG of genes expressed in proliferative tissues and demonstrated experimentally the ability of various CRMs to recapitulate individual components of the overall expression pattern of the SG. Finally, we systematically examined the chromosomal localisation of genes from each SG and discovered that some SGs contain genes that are clustered at the same genomic locus, a feature of SGs that has not been previously described.

The iCab strain was derived by successive brother-sister inbreeding of wild-type Oryzias latipes from a closed stock. iCab, Cab and Heino (Loosli et al., 2000) strains were maintained as described previously (Koster et al., 1997). Embryos were staged according to Iwamatsu (Iwamatsu, 1994).

The automated ISH screen on whole-mount embryos was performed at stages 18, 24 and 32 of development as described (Quiring et al., 2004). Probes were generated from clones derived from a full-length cDNA library of 50,000 clones (sequenced by the Yunhan Hong group, National University of Singapore).

The data present in the database were translated into a binary array of gene expression values. The substructures of a given tissue were assigned an expression value inherited from the parental tissue following the hierarchical design of the anatomical ontology. For instance, if a gene was annotated as expressed in the ‘eye’, all the eye substructures, i.e. ‘retina’, ‘lens’, etc., will be attributed the value 1. Genes with no expression were excluded from this array. To increase the number of expression patterns taken into account to generate the array, MEPD gene annotations for stages 18 and 19 and for stages 30 and 32 were merged, corresponding to 1 dpf and 4 dpf, respectively, as there were no significant morphological differences between these respective stages.

TERRAIN (TRN) data reduction (Saeed et al., 2003; Stuart et al., 2003; Saeed et al., 2006) was performed using the binary matrix of expression patterns as input, with the following parameters: Euclidean distance; TRN initialisation on genes; number of closest neighbours, ten. The TRN ‘mountains’ were extracted by individually selecting hubs of linked genes. These hubs consist of an independent group of at least four linked genes at a default distance weight threshold ≥0.8. The randomisation of each matrix was performed by shuffling the values within each row (gene) and by keeping the column names identical (anatomy ontology). The mean profile of each TRN mountain was calculated by averaging the expression values of the genes contained in the cluster. We then performed a hierarchical clustering of the matrix created using the EBI expression profiler (Kapushesky et al., 2004) to generate the hierarchical tree that was used to browse the co-expression network. Parameters: Euclidean distance, average linkage clustering, gene tree-based clustering.

We used ErmineJ software (Lee et al., 2005; Gillis et al., 2010) to detect shared over-represented biological terms from Gene Ontology (GO) (Ashburner et al., 2000) [P-value threshold t=0.05, multiple-test correction by Benjamini-Hochberg (default)]. We compiled GO annotations of medaka genes by retrieving the GO annotations of the orthologous vertebrate genes.

Medaka EST sequences from MEPD were blasted against medaka genes in Ensembl (BLASTX, e-value<10^–7) version 46 (Hubbard et al., 2007). For each medaka Ensembl gene, the orthology assignment by comparison with other vertebrate species was performed using the Ensembl Perl API for the Compara database, restricting the orthology search to the terms: ortholog_one2one, ortholog_one2many, inspecies_paralog and apparent_ortholog_one2one. For all the medaka genes, the non-coding regions were retrieved using the Perl API for the Ensembl Core database.

We used the Trawler program (Ettwiller et al., 2007; Haudry et al., 2010) to identify conserved over-represented motifs in each node of the hierarchical clustering tree. Trawler was run with the default options and the number of occurrences (K) was calculated according to the number of genes composing each node (n) as: K=n/2. Orthologous sequences for the conservation calculation were retrieved from the Ensembl database.

Perl scripts using the API Core of Ensembl v48 were written to retrieve the number of genes on each chromosome and to obtain their genomic coordinates.

We have built a repository of medaka synexpression groups, which is available at http://zooserv1.zoo.uni-heidelberg.de/mirana. Included are SGs where there are at least four genes that share the exact same occurrence of the most over-represented motifs.

Regulatory regions were cloned using the following primers (5′-3′, forward and reverse): cdon, TGTACGCTGCAGTTAGCGCC and AAGTGCAGGCCAGACACGCA; otx3, TTTGCATCCTTACACTCTGGTATTCC and ACTCTGTCCAAATAAACCAA; hmgb2, TGTCCTGTGTTTCTTTGCATTTGA and TGGAGCTTGACGCTTAACGG. The corresponding PCR products were subcloned into pCRII-TOPO (Invitrogen) and subsequently inserted into a pBlueScript-based transgenesis vector containing two recognition sites for the meganuclease I-SceI flanking a multiple cloning site and a 3′ cassette containing a minimal promoter, enhanced GFP and an SV40 polyadenylation signal. We used the p35S minimal promoter for cdon and the hsp70 minimal promoter for hmgb2. No minimal promoter was used for otx3 as the regulatory region was cloned including its endogenous minimal promoter.

Injections into one-cell stage medaka embryos were performed as previously described (Rembold et al., 2006). Transgenic embryos were fixed overnight at 4°C in 4% paraformaldehyde, washed three times in phosphate-buffered saline containing 0.1% Tween 20 and dechorionated. Embryos older than stage 24 were grown in 0.003% N-phenylthiourea (Sigma-Aldrich) to inhibit pigmentation. Samples were mounted using 1% SeaPlaque GTG Agarose (Lonza) and glass-bottom culture dishes (MaTek). Images were obtained using a Leica TCS SPE confocal microscope, 488 nm excitation line and ACS APO 10×/0.30 objective lens. Image stacks of up to 35 confocal slices were transformed into maximum intensity projections using ImageJ software (Version 1.41o; http://rsbweb.nih.gov/ij/). Colour balance was minimally corrected and the filter Unsharp Mask [radius, 1.0 pixels; mask weight, 0.2) was applied to all images.

In order to systematically identify and analyse SGs, we combined experimental and bioinformatics approaches (summarised in Fig. 1). In brief, gene expression patterns and their respective annotations were deposited in the Medaka Expression Pattern Database (MEPD). Annotation according to a standardised medaka anatomical ontology was used for clustering genes with similar expression patterns. GO annotations and non-coding sequences of medaka genes were retrieved in order to identify common biological functions and shared DNA motifs among similarly expressed genes.

We collected the spatiotemporal expression of randomly picked medaka genes with cDNA probes from a library of full-length sequenced clones in an automated screen on whole-mount preparations using in situ hybridisation (Fig. 1). A total of 782 clones were used to probe gene expression in embryos at three developmental time points: 1 day post-fertilisation (dpf) (late neurula stage); 2 dpf (16 somite stage); and 4 dpf (somite completion stage). Of these, the expression patterns of 407 clones were manually annotated and integrated into MEPD. By combining the data obtained from this screen and the pre-existing MEPD records from previous ISH screens (Nguyên et al., 2001; Quiring et al., 2004), the expression patterns for 750 UniGene clones, corresponding to 560 Ensembl genes, have been recorded with precise descriptions in at least one embryonic stage.

To identify groups of co-expressed genes, we performed pairwise comparisons of expression patterns documented in MEPD, using the medaka ontological annotations. Genes sharing highly similar annotations were regrouped into co-expression ‘mountains’ using the TERRAIN (TRN) algorithm (Saeed et al., 2003; Saeed et al., 2006) (Fig. 2); 18, 32 and 63 TRN mountains were obtained for 1, 2 and 4 dpf, respectively (supplementary material Tables S6-S8). Consistent with the tissue complexity that arises during development, we observed an increasing number of TRN mountains, reflecting emergent domains of co-expression as development proceeds. Interestingly, this increase in co-expression domains does not arise from new specific gene expression but mainly from a restriction of early ubiquitous expression into late specific expression. Indeed, we observed that although a large proportion of genes showed ubiquitous expression at early stages (58% at 1 dpf, 47 % at 2 dpf), they became spatially restricted at later stages, particularly into the developing neural tissue and central nervous system (CNS) (supplementary material Fig. S1). By 4 dpf, only 7% of the genes remain ubiquitously expressed. This specific landscape does not occur by chance, as similar analysis on randomised annotations revealed no changes in landscape (supplementary material Fig. S2; see Materials and methods).

If they shared a common domain of expression, the co-expression mountains were further grouped into ‘mountain ranges’ by hierarchical clustering (see Materials and methods and http://zooserv1.zoo.uni-heidelberg.de/mirana). For instance, at 4 dpf, the TRN mountain of genes expressed in ‘somites’ was linked to the TRN mountain of genes expressed in ‘somites, pectoral fin’ to form a mountain range (Fig. 2). These two TRN mountains might have a common mode of regulation as they share a common expression domain (i.e. somites) and therefore might participate in the same SG. By systematic browsing through each node of the hierarchical tree, at 4 dpf 74 nodes representing TRN mountains or mountain ranges were identified to share common expression domains representing potential SGs (supplementary material Tables S1-S4).

We investigated whether co-expressed genes belong to the same biological process by calculating GO term over-representation in the 74 groups. Of these, 30 (41%) shared statistically significant over-represented biological categories, thereby complying with the definition of SGs by Niehrs and Pollet (Niehrs and Pollet, 1999) and constituting novel medaka SGs.

We surveyed the literature to substantiate the correlation between the observed expression patterns and the over-represented GO categories (Fig. 3). The GO term ‘cytoskeleton’, which is over-represented in the enveloping layer (EVL) SG (Fig. 3A), reflects the active cell migration that EVL cells undergo during epiboly (Zalik et al., 1999). ‘Gas exchange’ (such as oxygen transport) is the primary function of the cardiovascular system (Burggren and Pinder, 1991) and is the over-represented biological process in the ‘cardiovascular system’ SG, which includes well-known genes such as α-globin (Fig. 3B). The ‘contractile fibre’ GO term is over-represented in the ‘somites’ and ‘somites, pectoral fin’ SGs, matching the ontogeny of these muscle cells (Fig. 3C). More strikingly, we identified an SG that contains genes that are expressed in the ciliary marginal zone (CMZ) of the retina or in the tectum proliferative zone (TPZ) of the brain (Fig. 3D). These neural domains are known to be highly proliferative and are suspected to contain stem cell niches (Wittbrodt et al., 2002; Alunni et al., 2010). Accordingly, the ‘DNA-dependent DNA replication’ biological process is over-represented in this SG, suggesting a role in cell proliferation as exemplified by the presence of proliferating cell nuclear antigen (pcna) in the SG.

Finally, biological categories related to ribosomal proteins (‘constituent of ribosome’, ‘ribosomal subunit’, ‘cytosolic ribosome’), were the most over-represented GO categories in 18 SGs. Genes within the medaka ‘ribosome biogenesis’ SG were ubiquitously expressed early on and became restricted in neural derivative tissues. This dynamic spatiotemporal expression pattern corroborates the specific role of ribosomal genes during development, more particularly in pan-neural derivative tissues. Indeed, in C. elegans it has been suggested that ribosome biogenesis plays a role in the regulation of developmental processes (Saijou et al., 2004). In zebrafish, knockdown of 21 ribosomal proteins led to tissue-specific defects (Uechi et al., 2006). Similarly, in mouse, ribosomal protein Rpl38 knockout also led to a tissue-specific phenotype (Kondrashov et al., 2011). Several medaka orthologues of mouse ribosomal proteins showing specific expression in that study were overlapping with the genes in our ‘ribosome biogenesis’ SG (namely rpl13a, rpl36, rps4, rps3a and rps8).

In summary, we have identified groups of co-expressed genes with a shared biological function that correlates with their expression pattern, in agreement with the definition of SGs. Some of the identified SGs have not been described previously (e.g. ‘cardiovascular system’), whereas others were also found in Xenopus, suggesting evolutionary conservation [e.g. ‘cell cycle’ (Baldessari et al., 2005), ‘muscle and epidermis’ (Pollet et al., 2005), ‘ribosome biogenesis’ (Wischnewski et al., 2000)].

Although independent studies in specific SGs have highlighted that these genes are co-regulated by common t-AFs, the generality of these principles to other SGs remains unproven. We hypothesised that: (1) if all the genes within an SG are regulated by the same t-AF, we would expect them to harbour a common DNA sequence (or motif) in their regulatory regions that represents the binding site of the t-AF; and (2) this shared motif should be enriched in the set of genes within the same SG, as compared with its occurrence in the regulatory regions of other genes. Using the Trawler program (Ettwiller et al., 2007; Haudry et al., 2010), we systematically searched for common sequence motifs that represent likely t-AF footprints within the non-coding 5 kb upstream of the transcription start site (TSS) of genes encompassed in the 74 co-expressed TRN mountains or mountain ranges.

We identified significantly over-represented motifs de novo in 68% (50/74) of the co-expressed groups (supplementary material Tables S1, S2). About half (48%, 24/50) of these co-expressed groups fit the definition of SGs as they also share a common over-represented GO category (supplementary material Table S1). Interestingly, 52% of the co-expressed groups sharing an over-represented motif did not display an over-represented GO category (supplementary material Table S2). This latter group offers great potential for identifying new SGs because not all the medaka genes have been annotated with GO terms [only 20% (3824/18076) have been associated with a GO category to date]. We observed that SGs share a large number of over-represented motifs (20 on average), suggesting that they could be regulated by multiple t-AFs (Fig. 4A).

We examined the evolutionary conservation of these motifs by running Trawler with orthologous vertebrate sequences. Previous studies have shown that evolutionarily conserved regulatory sequences are more often present around developmental genes (Woolfe et al., 2005). As expected, we identified conserved putative regulatory motifs: 21% of shared motifs are conserved among teleosts and 1% among vertebrates (Fig. 4B). However, the largest proportion of newly discovered motifs were medaka specific (78%), which suggests that non-conserved DNA motifs can potentially have a regulatory function.

To further investigate the functionality of this shared input, we focused on the proliferative zone SG for experimental validation. We noticed that the positions of the shared motifs along the upstream regions of the genes are not random but are rather clustered into modules (Fig. 5, light-pink squares). To confirm this observation, we investigated whether there is a bias in the distances between the positions of two consecutive motifs, or whether they are randomly distributed along the putative regulatory region. Comparing the distances between motifs in the SG against distances in a dataset of random positions, we observed an enrichment of clustered motifs separated by short distances (less than 100 bp) in the SG (Fig. 4C). This modular organisation of the over-represented motifs in SGs supports their regulatory potential, as transcription factors are well known to interact with other co-factors, likely aided by colocalisation into compact CRMs (Howard and Davidson, 2004; Wasserman and Sandelin, 2004). Moreover, the multiple occurrence of a de novo identified motif within the same CRM is in agreement with previous studies that have shown that local clustering of transcription factor binding sites is a feature of vertebrate enhancers (Gotea et al., 2010).

To gain further insight into the relevance of the newly identified motifs, we investigated the specific signalling pathways involved in the regulation of the ‘proliferative’ SG and attempted to link these pathways with specific instances of the de novo motifs. On the one hand, we performed a signalling pathway enrichment analysis using DAVID (version 6.7) (Huang et al., 2009). The first enriched pathway from the KEGG database was ‘DNA replication’, confirming our previous GO analysis (supplementary material Fig. S3A). Interestingly, the first enriched pathway from the BIOCARTA database is the ‘p53 signalling pathway’, of which three well-described genes, namely e2f, cdk2 and pcna, are present in this SG (supplementary material Fig. S3B). On the other hand, the first over-represented motif identified with Trawler (Fig. 3D, lower panel; supplementary material Fig. S3C, upper panel) was predicted to overlap with the known Foxl1 binding site (MA0033 from the Jaspar database) (Portales-Casamar et al., 2010) (supplementary material Fig. S3C, lower panel), a transcription factor belonging to the class O forkhead box family. Strikingly, several studies have suggested a parallel role of the FOXO and p53 pathways: crosstalk between the p53 and FOXO transcription factors has been described and they have been demonstrated to share several common target genes (You and Mak, 2005). Taken together, our motif discovery and pathway analyses suggest that Foxl1 is a likely candidate to regulate genes belonging to the p53 signalling pathway. This example highlights the value of our pipeline to generate robust hypotheses for further experimental validation.

To experimentally validate that the CRMs encompassing shared DNA motifs can drive specific expression in related tissues, we assayed the activity of six CRMs from genes within the ‘proliferation’ SG (cdon, e2f, hmgb2, nup85, otx3 and rpf20). Three of them drove expression in the proliferative zones (Fig. 6).

Cdon (or Cdo, cell adhesion molecule-related/down-regulated by oncogenes) has been shown to promote the differentiation of several cell lineages in mouse (Oh et al., 2009), where it is initially expressed in the CNS and sensory organs and becomes restricted to the dorsal ventricular part of the brain (Mulieri et al., 2000), which is known to contain proliferative cells. This expression pattern is similar to that of medaka cdon, which is specifically restricted to the proliferative zones of the retina and brain at later stages of development. cdon harbours one conserved predicted motif in a region proximal to the TSS (Fig. 6A, red cross). We tested the regulatory activity of a 511 bp conserved region (Fig. 6A, green rectangle) containing this motif. We found that it recapitulates the endogenous expression pattern in the forebrain proliferative zone (FPZ) and midbrain-hindbrain boundary at 2 dpf, and also in the posterior part of the FPZ at 4 dpf (Fig. 6A). This exemplifies the ability of an evolutionarily conserved CRM encompassing the newly identified over-represented motif to drive specific spatiotemporal expression.

otx3 (or dmbx1) is a transcriptional repressor that plays a crucial role in brain morphogenesis (Zhang et al., 2002; Kimura et al., 2005). In zebrafish, it has been shown to regulate cell cycle exit during neurogenesis (Wong et al., 2010). This function correlates with its specific expression in the TPZ of the brain, which represents an area of proliferating as well as post-mitotic cells (Nguyên et al., 2001). Two conserved regions were identified in the upstream regulatory region of otx3, but only the region distal to the TSS contained an over-represented motif as predicted by Trawler (Fig. 6B, green circle). When we tested the regulatory activity of a 3 kb fragment (Fig. 6B, green rectangle) that contains both conserved regions, it recapitulated the endogenous expression pattern of otx3 in the TPZ. Deletion of the region with the over-represented motif, however, results in the complete loss of reporter expression (data not shown). In this case, a stretch of conserved DNA does not necessarily drive expression, but the conserved fragment that contains the identified over-represented motif is necessary to successfully recapitulate the endogenous expression pattern.

hmgb2 (high mobility group 2) is a cell proliferation promoting factor implicated in several oncogenic processes (Kwon et al., 2010). In accordance with this function, it is strongly and specifically expressed in all proliferative areas of the medaka embryo: the CMZ of the retina, the FPZ and TPZ of the brain, the rhombic lips and the marginal zones of the fins (Fig. 6C). In this case, we tested the activity of a 696 bp non-conserved CRM containing two copies of the same motif (Fig. 6C, purple triangles) as predicted by Trawler. This short fragment recapitulates a subset of the endogenous expression pattern in the anterior part of the FPZ (Fig. 6C) and in the ventral CMZ, the otic vesicles and pectoral fins (data not shown).

For the three other CRMs tested, the e2f CRM recapitulated only transiently the early ubiquitous expression pattern, the nup85 CRM did not drive a specific expression pattern and the rpf2 CRM drove an expression pattern different to that of the rpf2 gene itself. It is possible that these CRMs regulate other genes or act in synergy with other CRMs to perform their regulatory function.

The three positive examples shown above indicate that half of the computationally predicted modules containing de novo over-represented motifs drive specific gene expression.

It has been generally assumed that genes belonging to the same SG can be positioned anywhere in the genome, provided that they share a common regulatory input. We sought to systematically investigate whether this is a general property of SGs by examining the genomic location of genes within each of the 50 SGs that share over-represented DNA motifs. To our surprise, 30% (15/50) of 50 SGs sharing a common regulatory input contained pairs of genes located on the same chromosome and separated by fewer than ten genes (supplementary material Table S5). This spatial segregation is statistically significant as compared with random pairs of genes on the same chromosome (P<0.03 for 100,000 iterations, except for one pair of genes located on ultracontig88, where only 21 genes were annotated).

To further address the significance of this arrangement, we investigated whether the other genes lying between or near co-expressed pairs of genes also exhibit a similar expression pattern (Fig. 7; supplementary material Fig. S4). Six chromosomal locations (supplementary material Table S5, bold) revealed strikingly similar expression patterns in the newly investigated proximal genes. In particular, two distinct locations with clusters of genes expressed in the proliferative zones (Fig. 7A; supplementary material Fig. S4A) revealed four additional genes displaying identical expression patterns. In one of these clusters, the expression pattern of one gene (ENSORLG0000001937) is inverted, i.e. it is specifically repressed in the proliferative zones (supplementary material Fig. S4A), hinting at a common regulation pathway with the potential for both positive and negative effects on gene expression. Two loci include two gene pairs expressed in the ‘somites, pectoral fin’ SG (Fig. 7B; supplementary material Fig. S4B). One contains paralogous Troponin genes but also includes myoD, the role of which in muscle formation is very well characterised (Rudnicki et al., 1993). The expression of myoD in the somites was confirmed in the second round of ISH along with that of rfwd3 (ring finger and WD repeat domain 3), which is part of this cluster but has never been shown to be involved in muscle formation (Fig. 7B). Two additional genes were further confirmed by ISH to be co-expressed in the somites in the second locus (supplementary material Fig. S4B). ISH also highlighted co-expression of two intervening genes in a cluster of four genes previously grouped into the ‘EVL’ SG on chromosome 1 (Fig. 7C). Finally, two genes belonging to the ‘retina, forebrain, midbrain, cerebellum, pectoral fin’ SG (supplementary material Table S5) and located together on chromosome 12 (Fig. 7D) display striking similarity of expression, especially in the otic vesicle structure. These two genes are separated by a single gene, gfm2, which is expressed in the otic vesicle domain as judged by ISH analysis (Fig. 7D, white arrows).

One explanation for this chromosomal clustering is that these genes arose by gene duplication (Brunet et al., 2006). We therefore investigated whether the clustered genes in SGs share a common evolutionary origin and found that 12 SGs contain clustered co-expressed genes without any structural or ancestral relation. Only four SGs contain paralogous genes (supplementary material Table S5, italics); namely, the ‘yolk syncytial layer’ (YSL) SG with tandemly duplicated Apolipoprotein genes, the ‘cardiovascular system’ SG with a Globin gene cluster that arose from duplication of ancestral genes (Hardison and Miller, 1993), the ‘somites’ SG with a group of Troponin genes that are also arranged into clusters of paralogous pairs in human (Cullen et al., 2004) (Fig. 7B), and the ‘EVL’ SG cluster, with four of six uncharacterised genes sharing a conserved protein domain suggesting a paralogous relationship (Fig. 7C).

As the majority of synexpression clusters did not arise by gene duplication, we examined whether this spatial arrangement was selected under an evolutionary constraint by investigating the synteny conservation within the clusters. Four SG clusters display noteworthy synteny conservation between medaka and mouse [two ‘proliferative’ SGs (Fig. 7A; supplementary material Fig. S4A) and the ‘retina, forebrain, midbrain, cerebellum, pectoral fin’ SG (Fig. 7D)]. This conservation extends to human in the ‘somites, pectoral fin’ SG (supplementary material Fig. S4B). This locus is syntenic with human chromosome 19, which is known to harbour a concentration of skeletal muscle-specific transcripts (Bortoluzzi et al., 1998).

To reinforce the relevance of this genomic clustering, we investigated whether the expression pattern of the genes conserved in synteny is likewise conserved by screening the zebrafish and mouse expression pattern repositories [ZFIN (Bradford et al., 2011), GXD (Finger et al., 2011), Eurexpress (Diez-Roux et al., 2011)]. Strikingly, in all of the three clusters that show conserved synteny, we found at least two genes that display similar expression patterns in mouse (Fig. 7A,B,D, genes in green rectangles). For instance, genes expressed in the TPZ in medaka and zebrafish are also expressed in the subventricular zone in mouse, a region that is known to contain neural stem cells. Only two synexpression clusters display poor synteny conservation: the ‘somite’ SG contains only two genes with conserved synteny compared with mouse, which are the known Troponin paralogues (Fig. 7B), and within the ‘EVL’ SG the synteny conservation could not be assessed as no orthologous genes were found, suggesting that this is a cluster of medaka-specific genes (Fig. 7C).

In summary, we have shown that 30% of SGs include genes that are clustered at the same chromosomal locus, possibly representing an evolutionarily conserved feature along with the conservation of gene expression.

Although this represents a pilot study, it serves as a pioneer work to investigate general rules of SG regulation. The whole-genome study of SGs is mainly constrained by technical limitations in the simultaneous recording of precise spatiotemporal expression patterns during development. DNA microarray assays, for instance, have been successfully applied for studying genome-wide transcriptional regulation (Furlong et al., 2001; Baldessari et al., 2005) as they allow quantitative and temporal monitoring of expression during embryonic development, but lack the three-dimensional information concerning tissue connectivity within the embryo. Large-scale ISH on whole-mount embryos is an approach that can provide spatial information at high resolution for gene expression on an entire organism level (Lynch et al., 1995; Gawantka et al., 1998; Kawashima et al., 2000; Neidhardt et al., 2000; Kudoh et al., 2001; Pollet et al., 2005; Lecuyer et al., 2007). In Xenopus, microarray analysis not only confirmed SGs identified by ISH, but also allowed the discovery of novel SGs. Nevertheless, because they provide a low-resolution analysis of gene expression patterns compared with ISH, the risk of including false positives will increase (Baldessari et al., 2005). Studies combining both techniques would be ideal, as demonstrated in invertebrate models (Tomancak et al., 2002). Medaka is well suited as a vertebrate model organism on which to perform such an analysis because of its amenable properties for high-throughput screens (such as rapid and external development, embryo transparency) (Wittbrodt et al., 2002) combined with the completed assembly of its compact genome (Kasahara et al., 2007). The development of virtual embryos (Bryson-Richardson et al., 2007; Keller et al., 2008; Richardson et al., 2010) and subsequent three-dimensional expression pattern recordings offers great promise to accelerate the analysis of SGs.

We discovered several different DNA motifs that are shared between genes of the same SGs, although the combinatorial input differs from one gene to another. This observation is in accordance with several studies performed in other developmental model organisms that propose that similar expression patterns are directed by different regulatory codes. In Ciona, precise dissection of the regulation of muscle-specific genes highlighted a great variability in the number and types of shared regulatory motifs that are able to drive co-expression in this tissue (Brown et al., 2007). In D. melanogaster, combinations of binding events that direct gene co-expression were found to be variable in occupancy and intensity (Zinzen et al., 2009).

Nonetheless, this diversity of regulatory input might confer robustness against perturbations owing to the redundancy of the combinatorial input. Small changes in combinations of regulatory input or mutations in the regulatory repertoire might not be lethal for the embryo, but could trigger subtle expression pattern changes that would translate into phenotypic changes leading to animal diversity.

A key finding in our study is the identification of 12 genomic loci at which structurally non-related genes belonging to the same SG were clustered. Since this work only represents a pilot study, the prospect for identifying further synexpression clusters is high. The analysis of tissue-specific EST collections showed concentrations of multiple genes at the same loci (Bortoluzzi et al., 1998; Ko et al., 1998; Li et al., 2005; Vogel et al., 2005). Genes with common biological properties, such as highly expressed genes (Caron et al., 2001), housekeeping genes (Lercher et al., 2002), genes involved in the same biological pathway (Lee and Sonnhammer, 2003) and cancer genes (Zhou et al., 2003), were also shown to be clustered at the same chromosomal loci. The mechanisms that promote this chromosomal clustering remain largely unknown. Bi-directional promoters could explain why pairs of genes are co-expressed but do not explain why more than two genes share similar expression patterns (Cho et al., 1998; Cohen et al., 2000).

SGs have been considered as the eukaryotic counterparts of bacterial operons (Niehrs and Pollet, 1999). By analogy, co-expressed genes might be maintained along the same chromosome, owing to a single shared regulatory input, in order to lower the amount of regulatory information (Price et al., 2005). On the contrary, our data suggest that the regulatory information is redundant and combinatorial in SGs with genes clustered at the same locus. It is likely that the two mechanisms of regulation are not mutually exclusive, and that synexpression members are clustered at the same locus owing to a common mode of regulation. It is possible that they are co-regulated by global control regions (Spitz et al., 2003) that reside in a distal enhancer. Owing to computational restrictions, we limited our analysis to look for common over-represented motifs in the 5 kb upstream regions of the SG genes. In the future, increasing computational power will overcome this restriction. Given that the synexpression clusters display striking synteny and similar expression patterns, it is tempting to speculate that they might also show a sequential onset of gene expression, similar to the collinearity described for Hox gene clusters (Duboule, 1998).

We have performed the first systematic investigation of the regulatory properties of SGs. We initially identified medaka SGs and thereafter revealed that SGs share a common complex input that is potentially regulatory and embedded within CRMs. We also provide the first evidence that a significant proportion of SGs contain genes that are closely clustered at genomic loci. De novo unbiased approaches such as ours will aid in deciphering the complexity of interactions that occur in the regulation of gene co-expression during development.

We thank the J.W. laboratory for their valuable input on this project; Stuart Archer, Alvaro Mateos-Gil and Steffen Lemke for critical comments on the manuscript; Yunhan Hong for end-sequencing the 50 K medaka library; and Franck Bourrat for his valuable contribution to the annotations.