BACKGROUND
the specification and differentiation of cell types during animal development requires that genes be expressed in appropriate spatiotemporal patterns.
defining the regulatory mechanisms controlling this patterning is a central goal of developmental biology research.
one powerful tool to infer regulatory networks is to identify genes preferentially expressed in a cell type and screen experimentally or computationally for transcription factors  likely to bind those genes’ regulatory sequences.
this approach is especially powerful in model organisms such as worms and flies, whose smaller genomes reduce the amount of dna to search for regulatory function.
the nematode c. elegans is well suited for such a comprehensive study of developmental regulation because of its stereotyped development from zygote to adult, with each adult hermaphrodite developing through an identical pattern of cell divisions  <cit> .
each animal has the same number and organization of cells of each type, with  <dig> cells present at the end of embryogenesis.
in addition, the signaling pathways controlling cell type specification, including the notch , ras , and wnt  <cit>  signaling pathways, are conserved with humans and other animals.
time-lapse imaging of fluorescent reporters has generated cellular resolution expression information for many genes , and automated image analysis methods make it possible to identify all expressing cells in embryos or larvae  <cit> .
recent studies have defined the in vivo  <cit>  and in vitro  <cit>  binding and binding motifs  for a substantial proportion of c. elegans tfs, and have experimentally measured tf binding at scale in vivo  <cit>  and in vitro  <cit> , providing a basis for regulatory inference.
integrative analysis of coexpression, genetic and protein-protein interactions, and other data sources allow predicting the functions of many genes  <cit> .
imaging of animals using reporter genes  <cit> , rna fish probes  <cit> , or antibodies  <cit>  can detect developmental expression patterns across all cells of the embryo.
however, logistics limit the number of genes whose expression can be measured at high resolution by these methods.
alternatively, individual cell types can be isolated by flow cytometry from dissociated embryos  <cit>  or larvae  <cit> , and assayed for mrna levels genome-wide.
similarly, tissue-specific mrna can be isolated based on its association with an epitope-tagged poly-a binding protein expressed under the control of a tissue-specific promoter  <cit> .
these approaches have been applied to a subset of terminally differentiated cell types  <cit> , but a comprehensive analysis across cell types is limited by the lack of individual markers for most unique cells, and by the labor and cost associated with isolating and analyzing large numbers of cell types individually.
furthermore, even different cells of the same type  can have different expression profiles depending on their lineage history and position within the animal  <cit> .
previous studies of differential expression in the embryo assayed expression in terminally differentiated cell types, mostly as non-overlapping populations.
here, we developed a strategy, “profiling of overlapping populations of cells ” that uses expression measurements from overlapping cell populations to identify genes differentially expressed in arbitrary patterns.
we previously showed that measuring expression in multiple partially overlapping groups of cells can provide information about differential expression across the entire lineage, and is thus more comprehensive than sorting based on “cell type-specific” markers whose expression is minimally overlapping  <cit> .
here, we applied this concept to identify patterned gene expression across all cells of the c. elegans embryo by measuring expression genome-wide in multiple overlapping cell populations isolated by flow cytometry .
we show that these overlapping expression measurements provide broad information about where genes are expressed in the c. elegans embryo and we define  <dig> gene expression clusters, many of which correspond to groups of genes that are coregulated in particular tissues.
we identify  <dig> tfs whose motifs or in vivo binding are enriched near genes in  <dig> clusters; in many cases the putative regulators are coexpressed with their proposed targets.
we validate these findings by identifying novel gene expression and regulation in the pharyngeal glands and ciliated neurons, and by comparing with existing genomic resources.
these results identify general features of embryonic gene expression patterns and their regulation, and provide powerful resource for future studies of embryonic regulation.fig.
<dig> experimental strategy.
a summary: we facs sort embryonic cells, based on expression of markers with known expression patterns, and measure expression in cells expressing  a particular marker using rna-seq.
genes expressed in similar sets of cells are enriched in a similar set of samples.
b expression patterns of cells used for sorting , and in spencer et al.
<cit> , shown in yellow.
cell fates are shown in the colored bar at the top.
c expression pattern of unc- <dig>  in the abpl sublineage, with cell fates colored as in .
d comparison of overlap of groups of cells used for sorting in this paper, with similar overlap for the groups of cells used in spencer et al.
<cit>
RESULTS
selection and characterization of overlapping sort markers
we selected fifteen c. elegans transgenic reporter strains expressing gfp, mcherry, or both fluorescent proteins in specific embryonic cells   <cit> .
we identified all gfp or mcherry-positive cells in each strain through the hypodermal enclosure stage by automated lineage tracing of 4d confocal movies  <cit> .
this provided a cellular resolution atlas of each reporter gene’s expression, and identified new expressing cells and dynamics of expression for many reporters .
in general, the reporters used for sorting were expressed in multiple terminal cell types.
for example, pros-1::gfp, which was previously reported to be expressed and required in the excretory canal cell  <cit> , is also expressed in many sheath type glia cells, coelomocytes, pharyngeal glands and some neurons .
similarly, unc-130::gfp is expressed in progenitors of diverse cell types including a subset of muscle and hypodermal cells, the excretory system, several types of neurons and a few pharyngeal and rectal cells   <cit> .
the average overlap between our cell populations is much higher than in previous genome-wide analyses of cell-specific expression, which largely focused on distinct terminal cell types .
rna-seq from sorted cell populations reproducibly detects differentially expressed genes
we dissociated cells from embryos and used flow cytometry to purify cells based on these strains' fluorescent marker.
we analyzed both fluorescent “positive” cells and matched non-fluorescent “negative” cells from the same sort.
we prepared rna from each sample and quantified expression using strand-specific rna-seq on the solid platform  <cit> .
this resulted in nearly a billion mapped reads .
we detected expression of  <dig>  genes in at least one facs-sorted sample, at a level of at least one rpm , with between  <dig> and  <dig>  genes detected in each individual sample .
we detected more unique transcripts in cell populations with fewer embryonic cell types, and more genes with enriched or depleted expression , as compared with populations containing more embryonic cell identities.
this effect was significant for depletions .
this suggests that measuring transcriptomes in smaller groups of cells increases sensitivity to detect rare, cell-type-specific transcripts.fig.
<dig> data quality of expression measurements of facs-sorted cells.
a enrichment of genes in two replicates of sorting by a pha- <dig> reporter.
known pharyngeal genes defined as early or late embryonic in gaudet et al.
<cit>  are shown in red and blue, respectively.
b enrichment of mrnas corresponding to markers used for sorting.
promoter fusions are shown in red, while protein fusions are shown in green.
c comparison of number of sorted cell identities  in a sorted fraction with the number of genes enriched  or depleted .
d number of genes enriched or depleted in different numbers of sorted fractions.
e enrichment of time-specific genes in cells sorted by cnd- <dig>  the proportion of the total cells expressing the cnd- <dig> reporter is shown in blue.
f same as , except for cells sorted using a pros- <dig> reporter.
g comparison of tissues present in sorted fractions, with enrichment of known tissue-specific genes found by spencer et al.
<cit> .
there is one point for each pairing of a tissue t with a sort fraction s. the x-axis shows the enrichment of cells with known tissue t in sort fraction s. the y-axis shows the average enrichment in sort fraction s of genes annotated as expressed in tissue t by spencer et al.
<cit>
we identified genes whose expression was enriched or depleted in each group of facs-sorted cells by comparing each annotated gene’s expression between the positive sample and the paired negative control sample; we used mock-sorted cells as a control for two samples where the paired negative was not available.
both normalization methods gave similar results, but using the matched negative samples resulted in higher measured enrichment levels and thus an increased sensitivity to detect genes with modest expression enrichment .
these enrichments were reproducible across biological replicates for independent sort markers , indicating high overall reproducibility.
many genes were enriched or depleted in specific sort fractions;  <dig> genes were enriched or depleted 4-fold in at least one sample, and  <dig> were enriched or depleted in two or more samples .
this provides a conservative list of genes likely to be differentially expressed in the embryo .
expression of most marker genes  was enriched in their own positive sort fraction .
this enrichment was strongest for translational reporter markers where gfp is fused to the c-terminus of the protein and the gene is surrounded by its normal genomic context .
in contrast, enrichment was lower for “transcriptional” reporter markers where the marker gene’s promoter was used to drive a stable mcherry-histone fusion .
this may reflect the fact that many of the marker genes are expressed transiently during embryogenesis  <cit> , with the mcherry-histone fusion protein persisting long after the endogenous rna.
consistent with this, protein levels of the translational reporters often show dynamic regulation mirroring that of the corresponding mrna and are often expressed more transiently than promoter fusion reporters for the same gene  <cit> .
facs gating for single cells preferentially enriches for specific cell types
the combined expression of the positive and negative fractions was similar to, but not identical to, expression in bulk embryonic cells.
this could be because the forward-scatter and side-scatter “gates” used during cell sorting to exclude cell clumps also preferentially exclude certain cell types.
to test this, we compared expression between “singlet” cells that had been gated to exclude cell clumps with “ungated” cells that were run through the facs machine but not gated.
we identified  <dig> genes preferentially expressed in the ungated cells.
these genes were enriched for genes expressed in the hypodermis and intestine  <cit> , and in late embryonic cells .
accounting for the effects of singlet gating improves the similarity between ungated cells and the combined positive and negative expression profiles .
gating for single cells during flow cytometry thus provides information about an additional partially overlapping embryonic cell population that likely includes hypodermal and intestinal cells.
we therefore included “singlet enrichment” in the clustering analysis described below.
rna-seq from facs sorted cell populations identifies spatiotemporal gene expression signatures
since our lineage data identifies which cells should be contained within each sort fraction , we asked whether genes known to be expressed in specific cell types were enriched in the expected fraction.
in some cases cell types predicted to be present or absent in a given cell population have been previously characterized for genome-wide expression.
for example, the pha-4::gfp fraction specifically labels pharynx, intestine and rectal cells, and genes identified previously as expressed in the pharynx  <cit>  were preferentially expressed in that fraction .
we tested this more broadly by asking whether genes previously identified as tissue-specific by the modencode project  <cit>  were enriched in sort fractions that preferentially contain cells from that tissue .
we identified a significant relationship , consistent with the different fractions having the expected tissue compositions.
we identified many anatomy, expression, and gene ontology  annotation terms significantly associated with expression in specific sorted fractions .
each sorted fraction except for the singlet cells had at least one anatomy term significantly enriched .
these were generally consistent with the tissue identities of the cells present in that fraction .
similarly, many go terms enriched in particular fractions were consistent with the cell types present in each fraction and in some cases predicted novel gene classes .
for example, the mir- <dig> fraction, which preferentially contains hypodermal cells, was enriched for the anatomy term “hypodermis” and the go terms “structural constituent of cuticle” and “extracellular region,” consistent with the role of hypodermal cells in secreting the cuticular exoskeleton  <cit> .
similarly, the pha- <dig> fraction, which consisted mostly of pharyngeal cells, was enriched for genes associated with metalloendopeptidase activity.
such proteases have been implicated in remodeling of extracellular matrix during postembryonic organ growth  <cit> , and thus may also play a role in the developing pharynx, which undergoes complex morphogenetic changes and extracellular matrix remodeling  <cit> .fig.
<dig> annotation of facs-sorted cells.
enrichment of chip peaks, motifs, go terms, expression clusters, and anatomy terms associated with genes enriched in each sort fraction.
selected pha- <dig>  and mir- <dig>  enrichments mentioned in the text are boxed in red and blue, respectively
our smallest fraction was the ceh-6;hlh- <dig> double-positive cells, which consists of only four cells: the excretory duct and pore cells , and db <dig> and db <dig> motorneurons .
genes preferentially expressed in this fraction were enriched for annotations associated with db neurons .
intriguingly, this population was also enriched for genes such as grl- <dig>  grl- <dig>  ptr- <dig>  wrt- <dig>  grd- <dig>  and ptc- <dig>  which are associated with “hedgehog signaling”.
this pathway is not thought to be active for signaling in c. elegans  <cit> , but many genes with homology to the ligands and receptors exist.
some hedgehog-related genes have been shown to be expressed in epithelial cell types and consistent with this, a grl- <dig> reporter is expressed in the excretory duct and pore cells  <cit> .
genes homologous to the hedgehog receptor patched have been shown to be important for lumen formation in c. elegans glia  <cit> .
in total eleven genes associated with the hedgehog pathway were enriched in the ceh-6;hlh- <dig> double-positive cells.
this suggests the possibility that many hedgehog-related genes may be involved in lumen formation in the tubular cells of the excretory system.
while we chose our sort markers mostly with the goal of maximizing our ability to measure spatial patterns, these reporters may also contain information about the timing of gene expression.
we tested this by first identifying a group of “temporally-specific” genes expressed at different times in an rna-seq time-course from whole embryos  <cit>  .
we then asked whether these temporally-specific genes were enriched or depleted in each sort fraction.
fractions predicted to contain early embryonic cells had higher expression of “early genes,” while we observed higher expression of “late genes” in cells sorted based on reporters expressed later.
for example, cnd- <dig> is expressed in many cells early in embryogenesis, while pros- <dig> is expressed late, in a smaller fraction of cells.
genes expressed early in whole embryos tended to be enriched by cnd- <dig> sorting , while sorting by pros- <dig> depleted for early genes and was enriched for a subset of later genes .
each fraction was significantly enriched for specific temporal stages .
based on this, we conclude that our expression data includes information about both spatial and temporal expression differences between embryonic cells.
motif enrichment predicts regulators acting in each cell population
to identify tfs that may regulate genes in each cell population, we searched for tfs that preferentially bind near genes enriched in that population as measured by chip data from modencode  <cit> , and for tf motifs overrepresented upstream of the genes enriched in each fraction.
we compiled a database of  <dig> chip experiments from c. elegans  <cit>  and  <dig> tf motifs from multiple species  <cit>  including  <dig> motifs for  <dig> c. elegans tfs  <cit> .
this identified motifs and chip signals significantly associated with each of the seventeen facs sorting experiments .
these represent candidate regulators of gene expression within the cells in each population.
for example, upstream intergenic sequences of genes in the pha-4::gfp fraction are enriched for the foxa <dig> motif recognized by pha- <dig>  <cit> , and for binding of pha- <dig> as measured by chip  <cit> , consistent with the known role of pha- <dig> in pharynx cell identity and gene expression  <cit> .
the foxa <dig> motif, but not pha- <dig> chip binding, was also significantly enriched in genes expressed in the pal-1::gfp fraction, which contains a high fraction of rectal cells.
since pha- <dig> mutants have major rectal defects  <cit> , pha- <dig> may directly regulate many rectal-specific genes, similar to its role in the pharynx, but these genes may be less easily identified by chip on whole embryos because the rectum represents a much small fraction of all embryonic cells than the pharynx.
this suggests that the limitations of whole-organism chip in identifying regulators important for expression in small cell populations may be partially overcome by analysis of motif enrichment.
clustering overlapping sort fraction expression data identifies genes coexpressed across diverse embryonic cell types
since our experiments assayed expression in many partially overlapping populations of cells that collectively cover the full embryo , these data contain information about the expression patterns of every cell type  <cit> .
for example, pharyngeal gland cells are enriched in pha- <dig>  and pros- <dig>  fractions .
therefore genes preferentially expressed in pharyngeal gland cells should be enriched in these fractions and depleted in other fractions that do not contain these cells.
more generally, we predict that genes with similar patterns of enrichment and depletion across sort fractions are expressed in similar tissue-specific patterns.
we provide a web-based tool to allow users to find genes with an embryonic expression pattern similar to that of a query gene .
we used hierarchical clustering to identify groups of genes with similar expression patterns across all samples, suggesting they are coexpressed in the embryo .
we tested different correlation cutoffs for cluster inclusion, and selected a cutoff resulting in  <dig> clusters that maximized our ability to detect candidate regulators of clusters by motif and chip analysis .
we did not use the temporal rna-seq data from whole embryos  <cit>  as part of the clustering, but examining the temporal data for genes within a clusters makes it possible to predict the temporal order of expression for genes within a cluster .fig.
<dig> clustering of enrichment.
a average enrichment for genes grouped into  <dig> clusters.
the timeseries data is from li et al.
<cit> .
b msa lineage, showing expression of pha- <dig>  and pros- <dig> ; yellow indicates overlap.
pharyngeal gland cells are shown as red rectangles.
c cluster  <dig>  enriched with genes known to be expressed in pharyngeal gland cells.
d cluster  <dig>  enriched with genes known to be expressed in ciliated neurons.
e mean expression, and mean absolute enrichement, for each cluster.
clusters with known enriched anatomy annotation are shown in red; selected clusters are labeled.
f overlap of expressed and tissue-specific clusters
many of the clusters correspond to specific tissues, based on significant enrichment of previously annotated tissue specific genes as curated by wormbase in that cluster  <cit>  .
<dig> of the  <dig> clusters had at least one significantly enriched anatomy ontology term at an fdr of  <dig>  .
an additional  <dig> clusters were significantly enriched for tissue-specific genomic expression signatures representing  <dig> of the  <dig> embryonic tissues assayed by modencode  <cit>  .
finally,  <dig> clusters were significantly associated with one or more genome-wide expression datasets that did not explicitly assay tissue-specific expression.
these and the other clusters represent groups of genes that may be coregulated in distinct patterns not previously assayed by genome-wide methods, since such experiments can contain implicit information about cell type specificity  <cit> fig.
<dig> annotation of clusters.
a enrichment  of chip signals, tf motifs, go terms, expression clusters, and anatomy terms associated with genes in clusters.
b expression pattern of hlh- <dig> and nhr- <dig> in comma-stage embryos, measured by rna-fish.
c expression pattern of hlh- <dig> and ceh- <dig> in a three fold embryo, measured by rna-fish.
d enrichment of co-clustered genes in wormnet  <cit>  annotations
we tested whether genes in the same cluster are coexpressed across cells by comparing to the epic dataset of cellular resolution expression profiles for  <dig> genes and to existing larval patterns for  <dig> genes  <cit> .
for two genes with high-resolution expression data, genes in the same cluster had much more similar expression patterns than genes in different clusters .
this similarity was stronger for embryonic  than for the larval expression patterns .
this consistency is striking given that the rna-seq data includes information about later embryonic stages not assayed in the imaging data.
thus, known tissue-specific annotation and expression patterns support the idea that genes which occur in a given cluster are expressed in a similar set of embryonic cells.
we further validated the clusters by comparing them with wormnet  <cit> , which combines many c. elegans genomic resources in a network model.
genes in the same cluster were linked by annotations in wormnet 5- <dig> fold more often than random, depending on the annotation .
this enrichment was strongest for genes whose fly and yeast orthologs undergo protein-protein interactions, consistent with genes in a cluster acting together.
we also compared the resolution of our data with the spencer et al.
<cit>  tissue-specific dataset, by seeing if genes which clustered together in one dataset were tightly coexpressed in the other dataset .
we clustered the  <dig> embryonic experiments from spencer et al.
<cit> , and for each cluster, compared the average within-cluster correlation in that dataset, with the same correlation computed in our dataset .
we also did the reverse comparison .
in each case, many clusters of genes that had very similar expression patterns in one dataset often had different expression patterns in the other dataset.
for instance, the cluster of genes shown in additional file 6: figure s8c had a mean within-cluster correlation of  <dig>  in the spencer data, and  <dig>  in our dataset.
these genes were enriched in the a-type motor neuron sample in the spencer data, while in the facs data, several of the genes  appear somewhat different data.
possible explanations include that these genes are expressed in different subsets of a-type motor neurons, or differentially expressed in progenitor cells.
conversely, genes in a cluster which is highly correlated in our dataset  hlh- <dig> fractions; additional file 6: figure s8d) are not highly correlated in the spencer et al.
<cit>  dataset suggesting they correspond to a cell type not assayed in that dataset.
this analysis suggests that while our data and the spencer et al.
<cit>  overlap in their coverage, each dataset can find similarities between genes that the other dataset cannot resolve.
in some cases, a cluster is enriched for genes known to be expressed in a particular cell type but also predicts novel additional genes to be expressed in those cells.
for example, cluster  <dig> is defined primarily by high expression in the pha-4::gfp and pros- <dig> sorted fractions, and the only cells that are included in both of these fractions are the pharyngeal gland cells .
furthermore this cluster contains seven genes  of the sixteen known to be expressed in the pharyngeal gland cells   <cit> .
however, this cluster also contains an additional  <dig> genes; we predict that many of these are novel pharyngeal gland-expressed genes.
these genes are enriched for transcriptional regulators, especially nuclear hormone receptors, suggesting an important role for these factors in the gland cells.
we validated this by using single molecule rna fish  <cit>  to examine the expression of two tfs from this cluster: nhr- <dig> and ceh- <dig>  both of these genes showed expression overlapping with a reporter for the known regulator of gland cell development hlh- <dig> , indicating that they are also expressed in gland cells.
taken together, our results suggest that membership of a gene in a cluster associated with known anatomy terms is predictive that the gene is expressed in that part of the anatomy.
some clusters of the  <dig> contain mostly genes which are expressed at very low levels, suggesting they may represent molecular or technical noise.
other clusters have high expression but little variation between fractions suggesting they contain genes that are more ubiquitously expressed.
the clusters also differ in size, from  <dig> to  <dig> genes.
we predicted the tissue specificity of each cluster using the mean of the absolute value of enrichments across all the sorting experiments.
we observed that 86 % of the clusters that are enriched for known tissue-specific annotations had a mean absolute enrichment > <dig>  and log-expression > <dig>  based on this cutoff, we estimate that the genes in at least  <dig> of the clusters have cell type-specific expression .
only about half of these cell-type specific clusters were enriched for either anatomy ontology terms or previously described tissue specific expression ; this is not surprising, as existing annotations are limited for most cell types.
most c. elegans genes’ expression has not been characterized comprehensively across cells, and only a few cell types have been annotated with genome-wide approaches.
this suggests that although we only sorted for fourteen markers, the dataset contains information about a much larger number of cell types.
enrichment of motifs and tf binding predicts novel regulators of embryonic gene expression
if genes coexpressed in a cluster have common upstream regulators, motifs or binding of these regulators should be enriched in that cluster  <cit> .
each cluster thus provides an opportunity to identify cell-specific regulators based on enrichment of regulatory motifs or experimentally defined tf binding.
we tested each of the  <dig> motifs and  <dig> chip-seq data sets described previously for enrichment within upstream intergenic sequences of genes in each cluster.
we refer here to these upstream regions as “promoters”, but they likely include both promoter and enhancer elements.
we found  <dig> tf binding site motifs and  <dig> tf chip signals enriched in genomic sequence upstream of genes in at least one cluster, .
in many cases these motif enrichments were consistent with known regulators.
for example, cluster  <dig> is highly enriched for genes expressed in ciliated neurons .
promoters of genes in this cluster are significantly enriched for the x-box homeodomain motif  recognized by the c. elegans rfx homolog daf- <dig>  which is known to regulate expression in ciliated neurons  <cit> .
based on this motif enrichment, we predicted that other genes in this cluster are also regulated by daf- <dig>
similarly, genes in the pharyngeal gland cluster  and six other clusters associated with pharyngeal annotations were highly enriched for pha- <dig> chip binding and a forkhead binding motif predicted to be bound by pha- <dig>  consistent with the broad role of pha-4/foxa in regulating pharyngeal expression  <cit> .
the pharyngeal gland cluster was also enriched for an e-box motif predicted to be bound by hlh factors, likely hlh- <dig>  which has highly correlated expression enrichments to this cluster centroid and is known to regulate pharyngeal gland fate  <cit> .
we tested these predictions by using qpcr to examine the expression of three genes from each cluster in mutants for the predicted regulator  .
expression of all six predicted targets was reduced, with 67 %  reaching statistical significance, indicating they are regulated directly or indirectly by the predicted factors.fig.
<dig> predicted regulatory relationships.
a enrichment of rfx <dig> motif upstream of genes in cluster  <dig>  b enrichment of an eor- <dig> motif upstream of genes in cluster  <dig>  c significance of motifs being more or less conserved, or nearer or further from the tss  d expression of known  and predicted targets of daf- <dig> and hlh- <dig>  when either of those tfs is mutated.
e enrichment of tf-cluster pairs in y1h data from reece-hoyes et al.
<cit>
intriguingly, the rfx motif instances in cluster  <dig> were not uniformly distributed; instead they were highly biased towards positions close to the 5′ end of the annotated transcript , and in conserved sequences, as compared to rfx motifs near genes outside this cluster .
the enrichment in conserved sequences is consistent with the known functional importance of daf- <dig> in regulating these genes.
the enrichment near the transcription start site suggests that daf- <dig> primarily acts by binding promoter proximal regulatory elements rather than distal enhancers.
in contrast, we identified other cases where an enriched motif was preferentially located further from the 5′ end, suggesting it may act primarily in distal enhancers .
based on this, we tested different cutoffs for sequence conservation, and gene-motif distance, and found that tfs differ in those characteristics, with some enriched at specific positions or in conserved sequence, and others more uniformly distributed across upstream sequences .
motifs with the highest motif-cluster enrichments tended to be biased for locations further than 1 kb from the 5' end of the annotated transcript, and for higher conservation, although motifs for several other factors in addition to daf- <dig> were enriched for proximal locations as well.
we expect that some of the regulators of clustered genes will be expressed in similar patterns to their targets.
consistent with this, many known tissue identity regulators’ expression was highly correlated  with the centroid of a cluster containing genes expressed in that tissue, and also had its predicted binding motif significantly enriched in the same cluster.
in total we identified  <dig> tfs coexpressed with a cluster above a correlation coefficient threshold of  <dig>  and whose predicted binding motif was enriched at fdr corrected p <  <dig>  with  <dig> clusters, providing many novel candidate regulators for diverse embryonic cells.
for example, genes in the “coelomocyte” cluster  were enriched for the presence of a forkhead binding motif in their promoters , and expression of a forkhead tf predicted to bind that motif, let- <dig>  was highly correlated with that cluster’s centroid .
let- <dig> is known to be important for development of postembryonic-derived coelomocytes  <cit> , and our work suggests that it also regulates embryonic coelomocyte development.
other prominent examples include daf- <dig> in the ciliated neuron cluster, pha- <dig> in the pharyngeal gland cluster  and hlh- <dig> in a presumed body-wall muscle cluster .
direct binding of most c. elegans tfs to the promoters of other tfs has been assessed using yeast 1-hybrid interactions  <cit> .
we found that tfs that can bind to the promoters of one or more genes in a cluster were significantly more likely to have their motif enriched in that cluster, compared to random pairs of genes.
this enrichment was higher for some classes, such as homeodomain and zinc fingers .
this supports the idea that genes in our clusters often share biological functions, and that the motifs we find often correspond to actual regulatory relationships.
identification of cell type-specific patterns of noncoding rna expression
improved array and sequencing technology have revealed many expressed non-coding transcripts  <cit> , including long noncoding rnas  and rnas that are antisense to protein-coding genes   <cit> .
this noncoding transcription is often tissue-specific  <cit> .
our strand-specific rna sequencing data allowed us to differentiate the expression of non-coding rnas and nearby genes, even if they overlap on opposite strands.
consistent with previous studies, we saw higher average expression and facs enrichment for coding transcripts compared with previously annotated lincrna and ancrnas.
also consistent with previous studies  <cit> , we found a positive correlation  between expression of lincrnas and that of the nearest annotated coding gene.
we found a similar positive expression correlation between antisense “ancrnas” and the overlapping gene .
the mean and distribution of coexpression of lincrnas with neighboring genes is similar to that seen between adjacent protein coding genes   <cit> .
we found a similar positive expression correlation between antisense “ancrnas” and the overlapping gene .fig.
<dig> non-coding rnas.
a cluster containing linc- <dig>  linc- <dig> and genes with known neural expression patterns.
b-d correlation of  ancrnas,  lincrnas, and  all pairs of genes with their nearest neighboring gene
we identified many non-coding rnas that cluster with tissue specific genes .
for instance, cluster  <dig> includes linc- <dig> and linc- <dig>  and is also significantly enriched for genes expressed in pharyngeal ganglia .
since these genes are enriched in a similar set of fractions, we expect that they are also expressed in a similar pattern of cells.
in total,  <dig> lincrnas and  <dig> ancrnas were expressed at levels greater than  <dig> rpm in one or more samples within  <dig> clusters for which we have an annotated tissue type, and another  <dig> lincrnas and  <dig> ancrnas were expressed in putative cell type-specific clusters for which the tissue type is unknown .
we conclude that our data identifies many noncoding rnas likely to be differentially expressed in the embryo.
discussion
a path towards profiling gene expression genome-wide at cellular resolution across the entire organism
while previous studies focused on purifying specific cell types, our “profiling of overlapping populations of cells ” approach to identifying tissue specific expression in principle provides information about all embryonic cell types.
using facs sorting and rna-seq, we found groups of genes which are expressed across a panel of partially overlapping cell populations that collectively cover the entire embryo.
each embryonic cell type has a specific pattern of enrichment and depletion across these cell populations, allowing us to identify genes expressed in tissues that have not been profiled by genome-wide approaches, as demonstrated for the pharyngeal glands.
thus we have substantially expanded our knowledge of spatially patterned gene expression across the embryo.
our dataset provides a powerful resource for c. elegans developmental geneticists interested in specific cell types.
by identifying genes coexpressed with genes they already know are expressed in their cells of interest , they can identify potential players in the development of those cells.
similarly, the tfs for which binding or motifs are enriched near those genes provide candidate regulators of those cells’ development.
in our analysis, we identified coexpressed genes by clustering.
in principle, because our data include both the expression of each gene in each cell population and the identity of the cells in each population, it should be possible to predict the expression of each gene in each cell.
previous deconvolution experiments in other systems such as the arabidopsis root  <cit>  or zebrafish embryo  <cit>  estimated expression essentially on a grid.
in contrast, our strategy uses cell populations that share partially overlapping lineage histories to increase our resolution to identify differentially expressed genes.
we previously developed a computational unmixing strategy to perform this deconvolution  <cit> .
applying this approach to simulated data suggested that at least  <dig> fractions are needed for this strategy to yield accurate high-resolution expression patterns, but that smaller numbers of fractions can yield useful information about cell populations.
consistent with this, applying these methods to the data in this study gave statistically meaningful predictions.
for example, the cell-specific expression predictions resulting from deconvolution could be used to predict the expression within a cell fraction not used for deconvolution  .
while the accuracy of these predictions is not yet high enough to be useful at resolution approaching single cells, our results suggest that with additional data, the deconvolution approach may allow high-resolution inference of expression genome-wide across all cells.
implications for regulatory inference
this study identified substantially more motif-cluster enrichments than were observed previously.
for example,  <cit>  identified  <dig> motifs enriched upstream of clusters defined by coexpression across cell types  <cit> , while we identified  <dig> motif-cluster associations.
this could reflect differences in the information content of the underlying data, the clustering approach, or the motif database.
applying our clustering and motif enrichment approach to the data in  <cit>  or  <cit>  identifies a number of motifs comparable to what we found in this study, suggesting that the increased sensitivity reflects differences in the motif association itself.
the past few years have seen dramatic growth in our knowledge of experimentally defined tf binding specificities .
consistent with this, we identified many more enriched motif-cluster pairs from the experimentally determined binding sites than when we used fire  <cit> , a de novo motif finding algorithm .
clustering the spencer data similarly into  <dig> clusters, and running fire on the resulting clusters yielded a similar number of motifs  this is more motifs than spencer et al.
<cit>  found, but still smaller than the number we found using known motif data.
thus, the recent influx of data on tf binding specificity provides a dramatic boost to regulatory inference.
our observation that motif enrichment was often biased towards particular levels of conservation or positions relative to the transcription start site suggests that new algorithms to integrate motif enrichment with these and other types of information  may further improve regulatory predictions.
CONCLUSIONS
in summary, we described here a new resource of genome-wide transcript level measurements for partially overlapping cell populations from c. elegans embryos.
we used this dataset to identify thousands of genes differentially expressed across the developmental lineage and to infer new regulators of embryonic development.
this resource provides useful data for future studies of embryonic gene regulation in c. elegans and the methods may be broadly useful in other systems where partially overlapping cell populations can be isolated.

 the specification and differentiation of cell types during animal development requires that genes be expressed in appropriate spatiotemporal patterns. defining the regulatory mechanisms controlling this patterning is a central goal of developmental biology research. one powerful tool to infer regulatory networks is to identify genes preferentially expressed in a cell type and screen experimentally or computationally for transcription factors likely to bind those genes’ regulatory sequences. this approach is especially powerful in model organisms such as worms and flies, whose smaller genomes reduce the amount of dna to search for regulatory function. the nematode c. elegans is well suited for such a comprehensive study of developmental regulation because of its stereotyped development from zygote to adult, with each adult hermaphrodite developing through an identical pattern of cell divisions. each animal has the same number and organization of cells of each type, with cells present at the end of embryogenesis. in addition, the signaling pathways controlling cell type specification, including the notch, ras, and wnt signaling pathways, are conserved with humans and other animals. time-lapse imaging of fluorescent reporters has generated cellular resolution expression information for many genes, and automated image analysis methods make it possible to identify all expressing cells in embryos or larvae. recent studies have defined the in vivo and in vitro binding and binding motifs for a substantial proportion of c. elegans tfs, and have experimentally measured tf binding at scale in vivo and in vitro, providing a basis for regulatory inference. integrative analysis of coexpression, genetic and protein-protein interactions, and other data sources allow predicting the functions of many genes. imaging of animals using reporter genes, rna fish probes, or antibodies can detect developmental expression patterns across all cells of the embryo. however, logistics limit the number of genes whose expression can be measured at high resolution by these methods. alternatively, individual cell types can be isolated by flow cytometry from dissociated embryos or larvae, and assayed for mrna levels genome-wide. similarly, tissue-specific mrna can be isolated based on its association with an epitope-tagged poly-a binding protein expressed under the control of a tissue-specific promoter. these approaches have been applied to a subset of terminally differentiated cell types, but a comprehensive analysis across cell types is limited by the lack of individual markers for most unique cells, and by the labor and cost associated with isolating and analyzing large numbers of cell types individually. furthermore, even different cells of the same type can have different expression profiles depending on their lineage history and position within the animal. previous studies of differential expression in the embryo assayed expression in terminally differentiated cell types, mostly as non-overlapping populations. here, we developed a strategy, “profiling of overlapping populations of cells ” that uses expression measurements from overlapping cell populations to identify genes differentially expressed in arbitrary patterns. we previously showed that measuring expression in multiple partially overlapping groups of cells can provide information about differential expression across the entire lineage, and is thus more comprehensive than sorting based on “cell type-specific” markers whose expression is minimally overlapping. here, we applied this concept to identify patterned gene expression across all cells of the c. elegans embryo by measuring expression genome-wide in multiple overlapping cell populations isolated by flow cytometry. we show that these overlapping expression measurements provide broad information about where genes are expressed in the c. elegans embryo and we define gene expression clusters, many of which correspond to groups of genes that are coregulated in particular tissues. we identify tfs whose motifs or in vivo binding are enriched near genes in clusters; in many cases the putative regulators are coexpressed with their proposed targets. we validate these findings by identifying novel gene expression and regulation in the pharyngeal glands and ciliated neurons, and by comparing with existing genomic resources. these results identify general features of embryonic gene expression patterns and their regulation, and provide powerful resource for future studies of embryonic regulation.fig. experimental strategy. a summary: we facs sort embryonic cells, based on expression of markers with known expression patterns, and measure expression in cells expressing a particular marker using rna-seq. genes expressed in similar sets of cells are enriched in a similar set of samples. b expression patterns of cells used for sorting, and in spencer et al., shown in yellow. cell fates are shown in the colored bar at the top. c expression pattern of unc- in the abpl sublineage, with cell fates colored as in. d comparison of overlap of groups of cells used for sorting in this paper, with similar overlap for the groups of cells used in spencer et al.  selection and characterization of overlapping sort markers we selected fifteen c. elegans transgenic reporter strains expressing gfp, mcherry, or both fluorescent proteins in specific embryonic cells. we identified all gfp or mcherry-positive cells in each strain through the hypodermal enclosure stage by automated lineage tracing of 4d confocal movies. this provided a cellular resolution atlas of each reporter gene’s expression, and identified new expressing cells and dynamics of expression for many reporters. in general, the reporters used for sorting were expressed in multiple terminal cell types. for example, pros-1::gfp, which was previously reported to be expressed and required in the excretory canal cell, is also expressed in many sheath type glia cells, coelomocytes, pharyngeal glands and some neurons. similarly, unc-130::gfp is expressed in progenitors of diverse cell types including a subset of muscle and hypodermal cells, the excretory system, several types of neurons and a few pharyngeal and rectal cells. the average overlap between our cell populations is much higher than in previous genome-wide analyses of cell-specific expression, which largely focused on distinct terminal cell types. rna-seq from sorted cell populations reproducibly detects differentially expressed genes we dissociated cells from embryos and used flow cytometry to purify cells based on these strains' fluorescent marker. we analyzed both fluorescent “positive” cells and matched non-fluorescent “negative” cells from the same sort. we prepared rna from each sample and quantified expression using strand-specific rna-seq on the solid platform. this resulted in nearly a billion mapped reads. we detected expression of genes in at least one facs-sorted sample, at a level of at least one rpm, with between and genes detected in each individual sample. we detected more unique transcripts in cell populations with fewer embryonic cell types, and more genes with enriched or depleted expression, as compared with populations containing more embryonic cell identities. this effect was significant for depletions. this suggests that measuring transcriptomes in smaller groups of cells increases sensitivity to detect rare, cell-type-specific transcripts.fig. data quality of expression measurements of facs-sorted cells. a enrichment of genes in two replicates of sorting by a pha- reporter. known pharyngeal genes defined as early or late embryonic in gaudet et al. are shown in red and blue, respectively. b enrichment of mrnas corresponding to markers used for sorting. promoter fusions are shown in red, while protein fusions are shown in green. c comparison of number of sorted cell identities in a sorted fraction with the number of genes enriched or depleted. d number of genes enriched or depleted in different numbers of sorted fractions. e enrichment of time-specific genes in cells sorted by cnd- the proportion of the total cells expressing the cnd- reporter is shown in blue. f same as, except for cells sorted using a pros- reporter. g comparison of tissues present in sorted fractions, with enrichment of known tissue-specific genes found by spencer et al.. there is one point for each pairing of a tissue t with a sort fraction s. the x-axis shows the enrichment of cells with known tissue t in sort fraction s. the y-axis shows the average enrichment in sort fraction s of genes annotated as expressed in tissue t by spencer et al. we identified genes whose expression was enriched or depleted in each group of facs-sorted cells by comparing each annotated gene’s expression between the positive sample and the paired negative control sample; we used mock-sorted cells as a control for two samples where the paired negative was not available. both normalization methods gave similar results, but using the matched negative samples resulted in higher measured enrichment levels and thus an increased sensitivity to detect genes with modest expression enrichment. these enrichments were reproducible across biological replicates for independent sort markers, indicating high overall reproducibility. many genes were enriched or depleted in specific sort fractions; genes were enriched or depleted 4-fold in at least one sample, and were enriched or depleted in two or more samples. this provides a conservative list of genes likely to be differentially expressed in the embryo. expression of most marker genes was enriched in their own positive sort fraction. this enrichment was strongest for translational reporter markers where gfp is fused to the c-terminus of the protein and the gene is surrounded by its normal genomic context. in contrast, enrichment was lower for “transcriptional” reporter markers where the marker gene’s promoter was used to drive a stable mcherry-histone fusion. this may reflect the fact that many of the marker genes are expressed transiently during embryogenesis, with the mcherry-histone fusion protein persisting long after the endogenous rna. consistent with this, protein levels of the translational reporters often show dynamic regulation mirroring that of the corresponding mrna and are often expressed more transiently than promoter fusion reporters for the same gene. facs gating for single cells preferentially enriches for specific cell types the combined expression of the positive and negative fractions was similar to, but not identical to, expression in bulk embryonic cells. this could be because the forward-scatter and side-scatter “gates” used during cell sorting to exclude cell clumps also preferentially exclude certain cell types. to test this, we compared expression between “singlet” cells that had been gated to exclude cell clumps with “ungated” cells that were run through the facs machine but not gated. we identified genes preferentially expressed in the ungated cells. these genes were enriched for genes expressed in the hypodermis and intestine, and in late embryonic cells. accounting for the effects of singlet gating improves the similarity between ungated cells and the combined positive and negative expression profiles. gating for single cells during flow cytometry thus provides information about an additional partially overlapping embryonic cell population that likely includes hypodermal and intestinal cells. we therefore included “singlet enrichment” in the clustering analysis described below. rna-seq from facs sorted cell populations identifies spatiotemporal gene expression signatures since our lineage data identifies which cells should be contained within each sort fraction, we asked whether genes known to be expressed in specific cell types were enriched in the expected fraction. in some cases cell types predicted to be present or absent in a given cell population have been previously characterized for genome-wide expression. for example, the pha-4::gfp fraction specifically labels pharynx, intestine and rectal cells, and genes identified previously as expressed in the pharynx were preferentially expressed in that fraction. we tested this more broadly by asking whether genes previously identified as tissue-specific by the modencode project were enriched in sort fractions that preferentially contain cells from that tissue. we identified a significant relationship, consistent with the different fractions having the expected tissue compositions. we identified many anatomy, expression, and gene ontology annotation terms significantly associated with expression in specific sorted fractions. each sorted fraction except for the singlet cells had at least one anatomy term significantly enriched. these were generally consistent with the tissue identities of the cells present in that fraction. similarly, many go terms enriched in particular fractions were consistent with the cell types present in each fraction and in some cases predicted novel gene classes. for example, the mir- fraction, which preferentially contains hypodermal cells, was enriched for the anatomy term “hypodermis” and the go terms “structural constituent of cuticle” and “extracellular region,” consistent with the role of hypodermal cells in secreting the cuticular exoskeleton. similarly, the pha- fraction, which consisted mostly of pharyngeal cells, was enriched for genes associated with metalloendopeptidase activity. such proteases have been implicated in remodeling of extracellular matrix during postembryonic organ growth, and thus may also play a role in the developing pharynx, which undergoes complex morphogenetic changes and extracellular matrix remodeling.fig. annotation of facs-sorted cells. enrichment of chip peaks, motifs, go terms, expression clusters, and anatomy terms associated with genes enriched in each sort fraction. selected pha- and mir- enrichments mentioned in the text are boxed in red and blue, respectively our smallest fraction was the ceh-6;hlh- double-positive cells, which consists of only four cells: the excretory duct and pore cells, and db and db motorneurons. genes preferentially expressed in this fraction were enriched for annotations associated with db neurons. intriguingly, this population was also enriched for genes such as grl- grl- ptr- wrt- grd- and ptc- which are associated with “hedgehog signaling”. this pathway is not thought to be active for signaling in c. elegans, but many genes with homology to the ligands and receptors exist. some hedgehog-related genes have been shown to be expressed in epithelial cell types and consistent with this, a grl- reporter is expressed in the excretory duct and pore cells. genes homologous to the hedgehog receptor patched have been shown to be important for lumen formation in c. elegans glia. in total eleven genes associated with the hedgehog pathway were enriched in the ceh-6;hlh- double-positive cells. this suggests the possibility that many hedgehog-related genes may be involved in lumen formation in the tubular cells of the excretory system. while we chose our sort markers mostly with the goal of maximizing our ability to measure spatial patterns, these reporters may also contain information about the timing of gene expression. we tested this by first identifying a group of “temporally-specific” genes expressed at different times in an rna-seq time-course from whole embryos. we then asked whether these temporally-specific genes were enriched or depleted in each sort fraction. fractions predicted to contain early embryonic cells had higher expression of “early genes,” while we observed higher expression of “late genes” in cells sorted based on reporters expressed later. for example, cnd- is expressed in many cells early in embryogenesis, while pros- is expressed late, in a smaller fraction of cells. genes expressed early in whole embryos tended to be enriched by cnd- sorting, while sorting by pros- depleted for early genes and was enriched for a subset of later genes. each fraction was significantly enriched for specific temporal stages. based on this, we conclude that our expression data includes information about both spatial and temporal expression differences between embryonic cells. motif enrichment predicts regulators acting in each cell population to identify tfs that may regulate genes in each cell population, we searched for tfs that preferentially bind near genes enriched in that population as measured by chip data from modencode, and for tf motifs overrepresented upstream of the genes enriched in each fraction. we compiled a database of chip experiments from c. elegans and tf motifs from multiple species including motifs for c. elegans tfs. this identified motifs and chip signals significantly associated with each of the seventeen facs sorting experiments. these represent candidate regulators of gene expression within the cells in each population. for example, upstream intergenic sequences of genes in the pha-4::gfp fraction are enriched for the foxa motif recognized by pha-, and for binding of pha- as measured by chip, consistent with the known role of pha- in pharynx cell identity and gene expression. the foxa motif, but not pha- chip binding, was also significantly enriched in genes expressed in the pal-1::gfp fraction, which contains a high fraction of rectal cells. since pha- mutants have major rectal defects, pha- may directly regulate many rectal-specific genes, similar to its role in the pharynx, but these genes may be less easily identified by chip on whole embryos because the rectum represents a much small fraction of all embryonic cells than the pharynx. this suggests that the limitations of whole-organism chip in identifying regulators important for expression in small cell populations may be partially overcome by analysis of motif enrichment. clustering overlapping sort fraction expression data identifies genes coexpressed across diverse embryonic cell types since our experiments assayed expression in many partially overlapping populations of cells that collectively cover the full embryo, these data contain information about the expression patterns of every cell type. for example, pharyngeal gland cells are enriched in pha- and pros- fractions. therefore genes preferentially expressed in pharyngeal gland cells should be enriched in these fractions and depleted in other fractions that do not contain these cells. more generally, we predict that genes with similar patterns of enrichment and depletion across sort fractions are expressed in similar tissue-specific patterns. we provide a web-based tool to allow users to find genes with an embryonic expression pattern similar to that of a query gene. we used hierarchical clustering to identify groups of genes with similar expression patterns across all samples, suggesting they are coexpressed in the embryo. we tested different correlation cutoffs for cluster inclusion, and selected a cutoff resulting in clusters that maximized our ability to detect candidate regulators of clusters by motif and chip analysis. we did not use the temporal rna-seq data from whole embryos as part of the clustering, but examining the temporal data for genes within a clusters makes it possible to predict the temporal order of expression for genes within a cluster.fig. clustering of enrichment. a average enrichment for genes grouped into clusters. the timeseries data is from li et al.. b msa lineage, showing expression of pha- and pros-; yellow indicates overlap. pharyngeal gland cells are shown as red rectangles. c cluster enriched with genes known to be expressed in pharyngeal gland cells. d cluster enriched with genes known to be expressed in ciliated neurons. e mean expression, and mean absolute enrichement, for each cluster. clusters with known enriched anatomy annotation are shown in red; selected clusters are labeled. f overlap of expressed and tissue-specific clusters many of the clusters correspond to specific tissues, based on significant enrichment of previously annotated tissue specific genes as curated by wormbase in that cluster. of the clusters had at least one significantly enriched anatomy ontology term at an fdr of. an additional clusters were significantly enriched for tissue-specific genomic expression signatures representing of the embryonic tissues assayed by modencode. finally, clusters were significantly associated with one or more genome-wide expression datasets that did not explicitly assay tissue-specific expression. these and the other clusters represent groups of genes that may be coregulated in distinct patterns not previously assayed by genome-wide methods, since such experiments can contain implicit information about cell type specificity fig. annotation of clusters. a enrichment of chip signals, tf motifs, go terms, expression clusters, and anatomy terms associated with genes in clusters. b expression pattern of hlh- and nhr- in comma-stage embryos, measured by rna-fish. c expression pattern of hlh- and ceh- in a three fold embryo, measured by rna-fish. d enrichment of co-clustered genes in wormnet annotations we tested whether genes in the same cluster are coexpressed across cells by comparing to the epic dataset of cellular resolution expression profiles for genes and to existing larval patterns for genes. for two genes with high-resolution expression data, genes in the same cluster had much more similar expression patterns than genes in different clusters. this similarity was stronger for embryonic than for the larval expression patterns. this consistency is striking given that the rna-seq data includes information about later embryonic stages not assayed in the imaging data. thus, known tissue-specific annotation and expression patterns support the idea that genes which occur in a given cluster are expressed in a similar set of embryonic cells. we further validated the clusters by comparing them with wormnet, which combines many c. elegans genomic resources in a network model. genes in the same cluster were linked by annotations in wormnet 5- fold more often than random, depending on the annotation. this enrichment was strongest for genes whose fly and yeast orthologs undergo protein-protein interactions, consistent with genes in a cluster acting together. we also compared the resolution of our data with the spencer et al. tissue-specific dataset, by seeing if genes which clustered together in one dataset were tightly coexpressed in the other dataset. we clustered the embryonic experiments from spencer et al., and for each cluster, compared the average within-cluster correlation in that dataset, with the same correlation computed in our dataset. we also did the reverse comparison. in each case, many clusters of genes that had very similar expression patterns in one dataset often had different expression patterns in the other dataset. for instance, the cluster of genes shown in additional file 6: figure s8c had a mean within-cluster correlation of in the spencer data, and in our dataset. these genes were enriched in the a-type motor neuron sample in the spencer data, while in the facs data, several of the genes appear somewhat different data. possible explanations include that these genes are expressed in different subsets of a-type motor neurons, or differentially expressed in progenitor cells. conversely, genes in a cluster which is highly correlated in our dataset hlh- fractions; additional file 6: figure s8d) are not highly correlated in the spencer et al. dataset suggesting they correspond to a cell type not assayed in that dataset. this analysis suggests that while our data and the spencer et al. overlap in their coverage, each dataset can find similarities between genes that the other dataset cannot resolve. in some cases, a cluster is enriched for genes known to be expressed in a particular cell type but also predicts novel additional genes to be expressed in those cells. for example, cluster is defined primarily by high expression in the pha-4::gfp and pros- sorted fractions, and the only cells that are included in both of these fractions are the pharyngeal gland cells. furthermore this cluster contains seven genes of the sixteen known to be expressed in the pharyngeal gland cells. however, this cluster also contains an additional genes; we predict that many of these are novel pharyngeal gland-expressed genes. these genes are enriched for transcriptional regulators, especially nuclear hormone receptors, suggesting an important role for these factors in the gland cells. we validated this by using single molecule rna fish to examine the expression of two tfs from this cluster: nhr- and ceh- both of these genes showed expression overlapping with a reporter for the known regulator of gland cell development hlh-, indicating that they are also expressed in gland cells. taken together, our results suggest that membership of a gene in a cluster associated with known anatomy terms is predictive that the gene is expressed in that part of the anatomy. some clusters of the contain mostly genes which are expressed at very low levels, suggesting they may represent molecular or technical noise. other clusters have high expression but little variation between fractions suggesting they contain genes that are more ubiquitously expressed. the clusters also differ in size, from to genes. we predicted the tissue specificity of each cluster using the mean of the absolute value of enrichments across all the sorting experiments. we observed that 86 % of the clusters that are enriched for known tissue-specific annotations had a mean absolute enrichment > and log-expression > based on this cutoff, we estimate that the genes in at least of the clusters have cell type-specific expression. only about half of these cell-type specific clusters were enriched for either anatomy ontology terms or previously described tissue specific expression; this is not surprising, as existing annotations are limited for most cell types. most c. elegans genes’ expression has not been characterized comprehensively across cells, and only a few cell types have been annotated with genome-wide approaches. this suggests that although we only sorted for fourteen markers, the dataset contains information about a much larger number of cell types. enrichment of motifs and tf binding predicts novel regulators of embryonic gene expression if genes coexpressed in a cluster have common upstream regulators, motifs or binding of these regulators should be enriched in that cluster. each cluster thus provides an opportunity to identify cell-specific regulators based on enrichment of regulatory motifs or experimentally defined tf binding. we tested each of the motifs and chip-seq data sets described previously for enrichment within upstream intergenic sequences of genes in each cluster. we refer here to these upstream regions as “promoters”, but they likely include both promoter and enhancer elements. we found tf binding site motifs and tf chip signals enriched in genomic sequence upstream of genes in at least one cluster,. in many cases these motif enrichments were consistent with known regulators. for example, cluster is highly enriched for genes expressed in ciliated neurons. promoters of genes in this cluster are significantly enriched for the x-box homeodomain motif recognized by the c. elegans rfx homolog daf- which is known to regulate expression in ciliated neurons. based on this motif enrichment, we predicted that other genes in this cluster are also regulated by daf- similarly, genes in the pharyngeal gland cluster and six other clusters associated with pharyngeal annotations were highly enriched for pha- chip binding and a forkhead binding motif predicted to be bound by pha- consistent with the broad role of pha-4/foxa in regulating pharyngeal expression. the pharyngeal gland cluster was also enriched for an e-box motif predicted to be bound by hlh factors, likely hlh- which has highly correlated expression enrichments to this cluster centroid and is known to regulate pharyngeal gland fate. we tested these predictions by using qpcr to examine the expression of three genes from each cluster in mutants for the predicted regulator. expression of all six predicted targets was reduced, with 67 % reaching statistical significance, indicating they are regulated directly or indirectly by the predicted factors.fig. predicted regulatory relationships. a enrichment of rfx motif upstream of genes in cluster b enrichment of an eor- motif upstream of genes in cluster c significance of motifs being more or less conserved, or nearer or further from the tss d expression of known and predicted targets of daf- and hlh- when either of those tfs is mutated. e enrichment of tf-cluster pairs in y1h data from reece-hoyes et al. intriguingly, the rfx motif instances in cluster were not uniformly distributed; instead they were highly biased towards positions close to the 5′ end of the annotated transcript, and in conserved sequences, as compared to rfx motifs near genes outside this cluster. the enrichment in conserved sequences is consistent with the known functional importance of daf- in regulating these genes. the enrichment near the transcription start site suggests that daf- primarily acts by binding promoter proximal regulatory elements rather than distal enhancers. in contrast, we identified other cases where an enriched motif was preferentially located further from the 5′ end, suggesting it may act primarily in distal enhancers. based on this, we tested different cutoffs for sequence conservation, and gene-motif distance, and found that tfs differ in those characteristics, with some enriched at specific positions or in conserved sequence, and others more uniformly distributed across upstream sequences. motifs with the highest motif-cluster enrichments tended to be biased for locations further than 1 kb from the 5' end of the annotated transcript, and for higher conservation, although motifs for several other factors in addition to daf- were enriched for proximal locations as well. we expect that some of the regulators of clustered genes will be expressed in similar patterns to their targets. consistent with this, many known tissue identity regulators’ expression was highly correlated with the centroid of a cluster containing genes expressed in that tissue, and also had its predicted binding motif significantly enriched in the same cluster. in total we identified tfs coexpressed with a cluster above a correlation coefficient threshold of and whose predicted binding motif was enriched at fdr corrected p <  with clusters, providing many novel candidate regulators for diverse embryonic cells. for example, genes in the “coelomocyte” cluster were enriched for the presence of a forkhead binding motif in their promoters, and expression of a forkhead tf predicted to bind that motif, let- was highly correlated with that cluster’s centroid. let- is known to be important for development of postembryonic-derived coelomocytes, and our work suggests that it also regulates embryonic coelomocyte development. other prominent examples include daf- in the ciliated neuron cluster, pha- in the pharyngeal gland cluster and hlh- in a presumed body-wall muscle cluster. direct binding of most c. elegans tfs to the promoters of other tfs has been assessed using yeast 1-hybrid interactions. we found that tfs that can bind to the promoters of one or more genes in a cluster were significantly more likely to have their motif enriched in that cluster, compared to random pairs of genes. this enrichment was higher for some classes, such as homeodomain and zinc fingers. this supports the idea that genes in our clusters often share biological functions, and that the motifs we find often correspond to actual regulatory relationships. identification of cell type-specific patterns of noncoding rna expression improved array and sequencing technology have revealed many expressed non-coding transcripts, including long noncoding rnas and rnas that are antisense to protein-coding genes. this noncoding transcription is often tissue-specific. our strand-specific rna sequencing data allowed us to differentiate the expression of non-coding rnas and nearby genes, even if they overlap on opposite strands. consistent with previous studies, we saw higher average expression and facs enrichment for coding transcripts compared with previously annotated lincrna and ancrnas. also consistent with previous studies, we found a positive correlation between expression of lincrnas and that of the nearest annotated coding gene. we found a similar positive expression correlation between antisense “ancrnas” and the overlapping gene. the mean and distribution of coexpression of lincrnas with neighboring genes is similar to that seen between adjacent protein coding genes. we found a similar positive expression correlation between antisense “ancrnas” and the overlapping gene.fig. non-coding rnas. a cluster containing linc- linc- and genes with known neural expression patterns. b-d correlation of ancrnas, lincrnas, and all pairs of genes with their nearest neighboring gene we identified many non-coding rnas that cluster with tissue specific genes. for instance, cluster includes linc- and linc- and is also significantly enriched for genes expressed in pharyngeal ganglia. since these genes are enriched in a similar set of fractions, we expect that they are also expressed in a similar pattern of cells. in total, lincrnas and ancrnas were expressed at levels greater than rpm in one or more samples within clusters for which we have an annotated tissue type, and another lincrnas and ancrnas were expressed in putative cell type-specific clusters for which the tissue type is unknown. we conclude that our data identifies many noncoding rnas likely to be differentially expressed in the embryo. discussion a path towards profiling gene expression genome-wide at cellular resolution across the entire organism while previous studies focused on purifying specific cell types, our “profiling of overlapping populations of cells ” approach to identifying tissue specific expression in principle provides information about all embryonic cell types. using facs sorting and rna-seq, we found groups of genes which are expressed across a panel of partially overlapping cell populations that collectively cover the entire embryo. each embryonic cell type has a specific pattern of enrichment and depletion across these cell populations, allowing us to identify genes expressed in tissues that have not been profiled by genome-wide approaches, as demonstrated for the pharyngeal glands. thus we have substantially expanded our knowledge of spatially patterned gene expression across the embryo. our dataset provides a powerful resource for c. elegans developmental geneticists interested in specific cell types. by identifying genes coexpressed with genes they already know are expressed in their cells of interest, they can identify potential players in the development of those cells. similarly, the tfs for which binding or motifs are enriched near those genes provide candidate regulators of those cells’ development. in our analysis, we identified coexpressed genes by clustering. in principle, because our data include both the expression of each gene in each cell population and the identity of the cells in each population, it should be possible to predict the expression of each gene in each cell. previous deconvolution experiments in other systems such as the arabidopsis root or zebrafish embryo estimated expression essentially on a grid. in contrast, our strategy uses cell populations that share partially overlapping lineage histories to increase our resolution to identify differentially expressed genes. we previously developed a computational unmixing strategy to perform this deconvolution. applying this approach to simulated data suggested that at least fractions are needed for this strategy to yield accurate high-resolution expression patterns, but that smaller numbers of fractions can yield useful information about cell populations. consistent with this, applying these methods to the data in this study gave statistically meaningful predictions. for example, the cell-specific expression predictions resulting from deconvolution could be used to predict the expression within a cell fraction not used for deconvolution. while the accuracy of these predictions is not yet high enough to be useful at resolution approaching single cells, our results suggest that with additional data, the deconvolution approach may allow high-resolution inference of expression genome-wide across all cells. implications for regulatory inference this study identified substantially more motif-cluster enrichments than were observed previously. for example, identified motifs enriched upstream of clusters defined by coexpression across cell types, while we identified motif-cluster associations. this could reflect differences in the information content of the underlying data, the clustering approach, or the motif database. applying our clustering and motif enrichment approach to the data in or identifies a number of motifs comparable to what we found in this study, suggesting that the increased sensitivity reflects differences in the motif association itself. the past few years have seen dramatic growth in our knowledge of experimentally defined tf binding specificities. consistent with this, we identified many more enriched motif-cluster pairs from the experimentally determined binding sites than when we used fire, a de novo motif finding algorithm. clustering the spencer data similarly into clusters, and running fire on the resulting clusters yielded a similar number of motifs this is more motifs than spencer et al. found, but still smaller than the number we found using known motif data. thus, the recent influx of data on tf binding specificity provides a dramatic boost to regulatory inference. our observation that motif enrichment was often biased towards particular levels of conservation or positions relative to the transcription start site suggests that new algorithms to integrate motif enrichment with these and other types of information may further improve regulatory predictions.  in summary, we described here a new resource of genome-wide transcript level measurements for partially overlapping cell populations from c. elegans embryos. we used this dataset to identify thousands of genes differentially expressed across the developmental lineage and to infer new regulators of embryonic development. this resource provides useful data for future studies of embryonic gene regulation in c. elegans and the methods may be broadly useful in other systems where partially overlapping cell populations can be isolated. 
our analysis provides new insight into embryonic gene regulation, and provides a resource for improving our knowledge of tissue-specific expression and its regulation throughout c. elegans development. understanding gene expression across the diverse metazoan cell types during development is critical to understanding their function and regulation. however, most cell types have not been assayed for expression genome-wide. hierarchical clustering analysis identified over sets of coexpressed genes corresponding to distinct patterns of cell type specific expression. we applied a novel approach we term “profiling of overlapping populations of cells ” to assay differential expression across all embryonic cells in the nematode caenorhabditis elegans. in this approach, we use rna-seq to define the transcriptome of diverse partially overlapping facs-sorted cell populations. c. elegansembryonic developmenttissue-specific expressionhttp://dx.doi.org/ /100000002national institutes of healthgm083145gm105676hg000046-13burdick joshua murray john isaac penn genome frontiers institutehttp://dx.doi.org/ /100004897pennsylvania department of healthissue-copyright-statement© the author 2016 we identified thousands of candidate regulators of these clusters based on enrichment of transcription factor motifs and experimentally determined binding sites. this identified thousands of transcripts differentially expressed across embryonic cells.
 understanding gene expression across the diverse metazoan cell types during development is critical to understanding their function and regulation. however, most cell types have not been assayed for expression genome-wide.  we applied a novel approach we term “profiling of overlapping populations of cells ” to assay differential expression across all embryonic cells in the nematode caenorhabditis elegans. in this approach, we use rna-seq to define the transcriptome of diverse partially overlapping facs-sorted cell populations. this identified thousands of transcripts differentially expressed across embryonic cells. hierarchical clustering analysis identified over sets of coexpressed genes corresponding to distinct patterns of cell type specific expression. we identified thousands of candidate regulators of these clusters based on enrichment of transcription factor motifs and experimentally determined binding sites.  our analysis provides new insight into embryonic gene regulation, and provides a resource for improving our knowledge of tissue-specific expression and its regulation throughout c. elegans development. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords c. elegansembryonic developmenttissue-specific expressionhttp://dx.doi.org/ /100000002national institutes of healthgm083145gm105676hg000046-13burdick joshua murray john isaac penn genome frontiers institutehttp://dx.doi.org/ /100004897pennsylvania department of healthissue-copyright-statement© the author 2016 