BACKGROUND
salivary gland  morphogenesis requires the complex coordination of cells to orchestrate a number of dynamic cellular processes including cell specification, lineage commitment, cell migration, proliferation and differentiation, all culminating in the formation of this specialized gland .
a network of signaling and regulatory molecules coordinates these vital biological processes, which are accompanied by the dynamic changes in gene expression throughout development.
in the mouse, submandibular salivary gland  morphogenesis occurs over several distinct developmental stages commencing at ~ embryonic day  <dig> .
at this early prebud stage, the primordial sg fate is established with the thickening of the adjoining oral epithelium.
at the next bud stage, which occurs at approximately e <dig> , the thickened epithelium invaginates into the underlying mesenchyme resulting in the formation of a primary bud, which serves as the precursor of the main duct of the salivary gland once the gland reaches maturation.
as the embryo develops to e <dig>  the gland undergoes rapid proliferation and intricate branching morphogenesis, during which the end buds undergo successive rounds of clefting resulting in the generation of multiple epithelial buds.
this pseudoglandular stage also coincides with reorganization of the end buds and the formation of the acini, which are the main secretory units of the salivary gland.
at the canalicular stage , the gland is highly branched with lumenization of the main secretory duct nearing completion.
the onset of gland cyto-differentiation also occurs at this stage – a process that continues to birth.
in the final step of embryonic salivary gland morphogenesis at the terminal bud stage , expansion of acini and lumenization of both ducts and acini is almost finished, resulting in a continuous network of ducts connecting the acini to the oral cavity .
while the gland is functional and able to secrete saliva at birth, further acinar maturation and differentiation continues postnatally, and by puberty, differentiation of the granular convoluted tubules is complete  <cit> .
for many years now, mouse genetic models have been widely utilized to study various facets of salivary gland biology including branching morphogenesis, cleft formation, organ development and differentiation.
while these studies have been instrumental in identifying some of the individual genes and signaling pathways necessary for proper salivary gland function, a limited number of studies have focused on the global examination of salivary gland gene expression in mouse .
in an early study, hoffman et al.
utilized microarrays to examine gene expression profiles at five different stages of mouse submandibular salivary gland embryonic development  <cit> .
in a more recent and complementary study, musselmann et al.
generated a microanatomical atlas of gene expression of embryonic salivary gland by performing laser capture microdissection of distinct epithelial populations obtained from e <dig> -e <dig> mouse embryos  <cit> .
these broad based approaches of global examination of gene expression profiles have been instrumental to deciphering molecular mechanisms of salivary gland morphogenesis and importantly in the discovery of novel signaling pathways such as fibroblast growth factor  and bone morphogenetic protein   <cit>  as well as key signaling molecules such as gsk3β  <cit>  that play important roles in regulating branching morphogenesis in embryonic salivary gland.
it is important to note that whereas gene profiling experiments during embryonic development of salivary gland has received much attention, similar in depth studies in adult has been lacking although one recent work has shed some light on aging associated sg gene signature  <cit> .
while microarray technology have been the application of choice in the past for transcriptome analysis, recent advancements have seen the supplanting of microarrays by genomic methods driven by next-generation sequencing  approaches like rna-sequencing   <cit> .
compared to microarray technology, deep sequencing based methods provide a more sensitive and precise analytical approach that can accurately quantify gene transcript levels and their isoforms across a broad dynamic range  <cit> .
indeed, direct comparisons to rna-seq has revealed the shortcomings of hybridization-based microarray studies such as background noise and saturation and probe set issues such as incorrect annotation and isoform coverage  <cit> .
with the recent explosion in rna-seq based genome wide analysis of the transcriptomic landscapes of multiple cell types and various tissues and organs, the time is ripe to use computational tools and a systems biology based approach to gain novel insights into tissue-specific gene signatures and their potential impact on biological processes.
in the present study, we have performed rna-seq to generate a comprehensive gene expression profile of the mouse submandibular gland at various time points of development during embryogenesis and maturation in adult.
to the best of our knowledge this is the first reported rna-seq based study to examine the transcriptome of the mouse submandibular gland.
an extensive bioinformatics analysis of our datasets has revealed interesting gene regulatory networks and maps that are enriched for and define the various stages of salivary glands.
finally we have leveraged the encode, fantom, human protein atlas as well as other published gene expression datasets to identify a salivary gland specific gene signature that is to a large extent conserved between mouse and human.
collectively our study not only validates existing literature but also provides a wealth of genomic resources that can be harnessed for the discovery of new genes and biologically important pathways in the salivary gland and for formulating testable hypotheses.
RESULTS
transcriptomic map of the mouse salivary gland during development and adult
in order to better define the dynamic changes in global gene expression levels and to identify new tissue-specific and tissue-enriched regulators, we performed rna-seq based expression profiling of the mouse submandibular salivary gland at various key stages of embryonic development in addition to post-natal and adult tissues.
towards this end, we isolated total rna from the mouse submandibular gland dissected from e <dig>  day old embryos , e <dig>   and e <dig>    <cit> .
we also examined submandibular glands isolated from day  <dig> , 4 week , and 12wk old male animals that represented different stages of post-natal maturation  <cit> .
biological replicates of these rna samples were sequenced to high depth using an illumina hiseq platform .
the sequence reads were mapped to the reference genome sequence of mus musculus  using tophat <dig> .
we subsequently performed between-sample normalization using the deseq median normalization method and calculated fragments per kilobase of transcripts per million  mapped reads thereby giving us measurements of relative expression of genes within and between biological samples.fig.
<dig> principal component analysis of the mouse salivary glands at different developmental time points.
a experimental scheme.
we isolated total rna from whole salivary glands ranging from early embryo to adult, and performed rna-seq.
utilizing these datasets, we defined and annotated the salivary gland transcriptional landscape by using various gene ontology  annotation analyses  and pathway analyses .
b proportion of variance in each principle component.
pc <dig>  pc <dig> and pc <dig> represent ~90% of variance in the data.
c projection plots show the pca coordinates for each stage, which are indicated by different colors.
the data indicates that the inherent variations in gene expression between biological samples can distinguish the developing salivary gland in a time dependent manner
in order to better analyze and appreciate the overall gene expression patterns between the various developmental and adult time points, we utilized principal component analysis , a statistical technique that reduces and summarizes large datasets while illustrating relationships between samples based on co-variance of the data being examined  <cit> .
using pca, we found that pc <dig>  pc <dig>  and pc <dig> accounted for approximately 90% of all variations of the original data .
in order to further explore and better depict the major sources of variation, all samples were plotted in a three-dimensional space consisting of pc <dig>  pc <dig> and pc <dig>  interestingly, as demonstrated in fig.
1c, each of the  <dig> representative time points datasets segregated into individual groups demonstrating the highly dynamic variation in gene expression between each sg sample.
indeed, biological replicates cluster tightly together, further highlighting the inherent similarity of the biological samples to one another.
another notable observation was that the embryonic samples clustered more closely to each other and that there is a striking separation between these and the adult samples.
taken together this analysis provided the first hint of a clear dichotomy of gene expression profiles between embryonic and postnatal salivary gland samples.
rna-seq analysis identifies a salivary gland gene signature
to evaluate development-dependent differential gene expression patterns, we next grouped our samples based on  <dig> distinct developmental stages corresponding to embryogenesis , neonatal  and adult .
for this analysis, we selected genes that showed at least a two-fold change in expression between each time point while showing an adjusted p-value of less than  <dig> .
we also considered genes that were expressed at ≥ <dig> fpkm in at least one biological replicate.
using this criterion, we identified a total of  <dig> differentially expressed genes  between the embryonic and neonatal stages with  <dig> genes found to be up-regulated and  <dig> genes showing downregulation .
similar comparison of degs between neonatal and adult developmental time points identified  <dig> genes.
of these degs,  <dig> were upregulated and  <dig> downregulated .
finally, we found  <dig> degs between embryo and adult with  <dig> number of genes showing up-regulation and  <dig> downregulated .
to better appreciate the biological relevance of the global transcriptomic differences between the embryonic and adult mouse salivary gland, we further analyzed the degs using clusterprofiler  <cit>  and identified pathways unique to each state .
interestingly, in all embryo enriched data sets, we observed specific enrichment of biological processes important during organogenesis, such as axon guidance as well as that of wnt and notch signaling pathways, both of which have been shown to be important for proper salivary gland morphogenesis .
in contrast, in the adult enriched data sets we observed enrichment in biological processes that include salivary secretion, protein export, fatty acid metabolism and vesicular signaling .
notably all of these processes have been shown to be critical for proper salivary gland function  and justifiably represent the mature stage of the adult gland.
armed with a global view of the overall transcriptional changes occurring between the  <dig> developmental stages, we next sought to identify a developmental stage-specific global gene expression profile of the salivary gland.
we reasoned that this broad gene profile might be useful in providing new insights into the biological processes and pathways that are important to each developmental stage.
to generate the developmental specific gene expression profile, we focused on genes that a) showed an fpkm of  <dig> or greater in at least one of the sample replicates across all samples and b) were comparatively differentially expressed, with a log2-fold change of  <dig> or greater, at one time point compared to other time points across all  <dig> datasets.
based on the above criteria we identified a total of  <dig> genes that were specifically enriched, across all datasets as demonstrated in fig.
2a.
to identify unique patterns of transcriptomic activities of the developmental stage-specific gene profile, we performed k-means clustering on the standardized log <dig> transformed fpkm values for each of the samples.
this analysis led to the identification of  <dig> different gene clusters that corresponded to genes enriched in a specific developmental time point as well as genes enriched in grouped time points, which we refer to as embryo, neonatal and adult enriched clusters .
interestingly, our analysis revealed that genes that were uniquely expressed during embryogenesis showed a gradual decline in expression as development progressed, while genes enriched in the neonatal stage showed a gradual increase in expression during embryogenesis with peak expression levels observed at the neonatal stage and then showing a gradual decline.
conversely, genes enriched in the adult stage showed a progressive increase during the developmental and neonatal stages with the highest levels of expression in the adult stage .fig.
<dig> cluster analysis of the salivary gland developmental profile.
a z-scores for each of the  <dig> genes were calculated and used as input for gene-wise k-means clustering analysis .
this analysis, visualized as a heatmap, depicts enrichment of genes at both specific development time points and general developmental stages .
b visualization of the gene z-scores from the general developmental stage enriched clusters depicts unique time-dependent patterns of expression.
in general, this analysis has identified genes that have peak expression in the defined developmental stages.
also show are selected enriched go-annotations for each developmental stage gene cluster
gene ontology  analysis of the developmental stage-specific gene profile identified overrepresentation of genes associated with specific functions, indicating the existence of distinct biological processes occurring during embryonic, neonatal and adult stages of salivary gland development.
a closer examination of genes that were uniquely expressed during embryogenesis revealed specific enrichment of biological pathways associated with organ, embryonic and vasculature development, cell differentiation, and transcription - all of these biological processes are associated with organ morphogenesis .
interestingly, we also observed specific enrichment in both nervous system development and nerve growth factor receptor signaling in salivary glands during embryogenesis.
this is in good agreement with studies demonstrating the critical role of the parasympathetic nervous system during salivary gland organogenesis .
genes specifically enriched in the neonatal stage also revealed association with interesting biological terms including regulation of secretion, lipid metabolic processes, catabolic processes, chloride transport, protein modification, and vesicular transport - all of which have been shown to play a role in secretory cell differentiation .
finally, analysis of the adult stage identified enriched biological pathways associated with salivary gland specific functions including saliva secretion, digestion, and ion and carbohydrate transport .
we also identified smooth muscle contraction as an enriched biological pathway in the adult sg , which fits very well with the existence of myoepithelial cell populations that are important for glandular contraction and expulsion of saliva.
taken together, the go analysis, confirms the validity and the possible functional significance of our sg developmental stage-specific gene expression profile .
since our analysis of the broad gene-enriched profile of embryonic, neonatal and adult salivary gland revealed interesting facets of its biology, we wondered if further in-depth parsing of the datasets could provide more detailed insights into salivary gland maturation, particularly in the embryonic stages.
hence, we next evaluated the transcriptomic activities from e <dig>  to e <dig>  to identify genes and pathways that may be critical during the early stages of salivary gland organogenesis and development.
of the  <dig> comparatively differentially expressed developmental specific genes we identified,  <dig> of them were specifically enriched during embryogenesis .
a closer analysis revealed that the greatest number of differentially expressed genes occurred at e <dig> , at a stage when the gland is undergoing rapid proliferation and has embarked upon a dynamic and elaborate program of branching morphogenesis.
while the e <dig>  functional analysis demonstrated broad developmental processes including salivary gland morphogenesis, exocrine system development, epithelium development, axonogenesis and axon guidance, we also observed selective enrichment in biological processes specific to the developmental programs that are at play at this critical time point during morphogenesis.
for instance, in agreement with the burgeoning branching morphogenesis that takes place at e <dig> , genes involved in pathways specific to epithelial development, branching morphogenesis, and salivary gland cavitation were overrepresented .
indeed, both branching and cavitation are critical for the proper development of the tubular network necessary for the transport of saliva to the oral cavity  <cit> .
while some degree of overlap was observed between e <dig>  and e <dig>  in processes involving cell differentiation and organ development, we did identify unique biological functions related to muscle function at e <dig> .
interestingly, we found specific enrichment in pathways involving muscle structure, muscle tissue development and regulation of muscle contraction, properties intimately associated with functions of the myoepithelial cells  <cit>  .
these findings are in good agreement with a previous study demonstrating the emergence of a myoepithelial cell population in the salivary gland at e <dig>   <cit> .
finally, our analysis at e <dig>  uncovered several biological pathways enriched for a variety of salivary gland specific functions including metabolic processes .
more specifically we observed over representation of metabolic processes important in amino acid metabolism and small molecule metabolism.
we also identified pathways involved in transport and proteolysis all of which have been shown to play roles in salivary gland biology and secretion .
taken together, our analysis suggests that the underlying global transcriptomic profile mirrors the distinct developmental stages of the embryonic sg as highlighted by branching morphogenesis , myoepithelial development  and secretory cell development .fig.
<dig> enriched biological process networks during embryogenesis.
a heatmap visualization of the z-scores from the  <dig> genes identified in the embryo specific  clusters generated from the analysis in fig.
<dig>  b–d) network visualization of selected top enriched biological processes  at e <dig> , e <dig>  and e <dig> .
the networks were assembled by cytoscape tool enrichmap, using an organic layout.
the node size represents the number of genes assigned to a biological process and edge width  is proportional to the number of overlapping genes between two nodes
meta-analysis of mouse tissue rna-seq expression data
our salivary gland rna-seq datasets offered a wealth of in-depth information regarding the transcriptomic repertoire of this tissue.
we rationalized that although many of the genes and pathways are common to sg and numerous other tissues and organs, there might exist some genes/pathways that are enriched for in expression and possibly function in the sg.
to probe this further, we examined rna-seq data for  <dig> adult organs and tissues generated by the mouse encode consortium  <cit> , and compared the overall gene expression profiles with those from the adult sg described in this study.
for this purpose, we generated a pairwise pearson correlation matrix using the top  <dig> genes with the highest median absolute deviation.
this, we reasoned, would allow us to identify genes with the most dynamic range of expression across all tissues.
our analysis revealed a clear separation of organs based on gene expression with the salivary gland clustering more closely with exocrine organs such as the pancreas and organs rich in stratified epithelial tissues including the skin, bladder and placenta .
the organ-specific segregation patterns were not a result of using a stringent cutoff of the top  <dig> genes, since a pairwise pearson correlation matrix using a larger set of  <dig>  genes, yielded similar results with the salivary gland continuing to cluster closely with the pancreas, skin, bladder and placenta .
a similar hierarchical clustering analysis also revealed a select list of genes that are similarly expressed between the salivary gland and the closely related four organs identified above  and those that are exclusive to the two exocrine glands; the salivary gland and pancreas .
interestingly, the transcriptomic identity of the mouse salivary gland was quite distinct from the mammary gland, another glandular tissue.
we suspect this is due to the fact the secretory activity of the mammary gland is limited only during parity and lactation rather than the virgin state, from which the rna-seq data was obtained.
in contrast, tissues representing specific brain regions such as the cerebellum and the mouse testis were surprisingly clustered close to the group containing the salivary gland and the four aforementioned organs – the significance of this finding is unclear and worthy of future exploration.
taken together, our findings demonstrate that at the broader global transcriptomic level, the salivary gland most closely resembles tissues and organs with similar morphological and functional characteristics.fig.
<dig> hierarchical clustering of mouse tissues.
fpkm values from the top  <dig> genes with the highest median absolute deviation were used to cluster adult mouse tissues .
the resulting heatmap shows that the salivary gland clusters closely with the pancreas, skin, bladder and placenta
having obtained a global gene expression profile of the salivary gland, we next wondered if we could derive an adult salivary gland specific molecular signature that is constituted by genes that show a distinct expression profile and hence are likely to be more relevant and important in salivary gland biology.
towards this end, we next compared gene expression levels across all adult mouse organs and tissues and identified  <dig> genes that were specifically enriched in the salivary gland and thus represented a potential salivary gland gene signature .
to gain a better understanding of the underlying biological functions and pathways associated with the adult salivary gland gene signature, we utilized cluego  <cit> , which is a functional gene ontology analysis tool that integrates several gene-set enrichment databases, including the kyoto encyclopedia of genes and genomes , reactome pathway database annotations, and the go consortium database, to create a comprehensive go/pathway term network  <cit> .
this analysis not only highlighted the close association of the signature to exocrine and salivary gland development but also identified enrichment in biological processes related specifically to protein export, protein processing and localization to the endoplasmic reticulum, all of which are quite relevant for proper salivary gland function .
interestingly, we also observed pathway enrichment specific to mucin biosynthesis – this is in good agreement with the established roles of mucins in salivary gland biology that includes their role as a first line of defense against microbial infection and their important contributions to various biophysical properties of saliva such as viscoelasticity  <cit> .fig.
<dig> visualization of the tissue specific salivary gland gene signature.
a hierarchical clustering  of the gene expression values selected from the tissue specific salivary gland signature across adult mouse tissues.
b network visualization of enriched pathways  in the gene signature was performed by cluego analysis
we next sought to identify a cohort of transcription factors, which are selectively enriched in our sg gene signature.
based on the crucial regulatory role transcription factors play in driving tissue-specific gene expression, we reasoned that this analysis may aid in identifying novel transcriptional regulators which may be important in salivary gland biology.
to ensure robust coverage, we first compiled a comprehensive mouse transcription factor list using databases from both riken  <cit>  and uniprotkb and then searched for transcription factors which are represented in our mouse adult salivary gland gene signature.
our analysis uncovered  <dig> transcriptional regulators .
of these transcription factors, ascl <dig>  <cit> , bhlha <dig>   <cit> , tfcp2l <dig>   <cit>  and six <dig>  <cit>  have been previously reported to play a role in salivary gland development and function.
our analysis also unearthed elf <dig>  <cit>  and ehf  <cit> , two members of the ets family of transcription factors, which are highly expressed in the salivary gland yet their specific function, if any, in salivary gland biology remains unknown.
moreover, we identified several other factors including eaf <dig> and foxi <dig> neither of which has been previously studied in the context of salivary gland biology.
the high expression of some of these transcription factors in salivary glands was also confirmed using publicly available datasets provided through the salivary gland molecular anatomy project .
in addition to the list of transcriptional regulators that were specifically enriched in our salivary gland gene signature, there were several other genes and gene family members that are both over-represented and likely or known to be associated with salivary gland biology .
this included genes belonging to the secretoglobin  gene superfamily, which encode for small secretory proteins found in high concentrations in secretions from various organs including the lacrimal and salivary glands  <cit> .
in addition, we observed selective enrichment in genes belonging to the kallikrein  family that encode for proteins which function as serine peptidases.
while a major function of this family of proteins have been associated with blood pressure regulation and skin homeostasis, some studies suggest that these enzymes may play a role in cell proliferation, cell survival and wound healing  <cit> .
finally, we also observed elevated levels of genes belonging to the demilune cell and parotid protein  family, which are a family of genes thought to possess antimicrobial activity  <cit> .
we posit that many of these genes unearthed by our systems biology based approach are likely to be strong candidates for various physiological functions associated with mouse salivary glands.fig.
<dig> visualization of selected members of the tissue specific salivary gland gene signature.
a hierarchical cluster analysis of the enriched transcription factors identified in the gene signature.
b hierarchical cluster analysis of enriched genes identified in the signature and which have been shown to play a role in salivary gland biology
evolutionarily conserved expression map between human and mouse salivary gland
mouse models have served a valuable role in elucidating molecular mechanisms of physiological and pathological processes in the salivary gland.
it is likely that the genetic interplays and gene regulatory networks that underlie these processes are mostly conserved between mouse and human.
to examine this we next asked if there were any commonalities in the salivary gland specific gene signature between mouse and human.
towards this end, we mined the publicly available transcriptome database for human tissues and organs generated by the human protein atlas  project  and the fantom <dig> consortium ).
to make proper comparisons, we only chose representative human tissue samples for which the corresponding data was available for the mouse tissues.
upon examination of the human protein atlas  samples, we found that of the  <dig> human genes that were homologous to the mouse,  <dig> genes  showed enrichment in the human salivary gland as compared to the  <dig> other tissues analyzed , thus recapitulating the tissue-restricted patterns of expression observed from our mouse tissue analysis.
one possible caveat to this comparison is that for the hpa dataset, there is limited information about the salivary gland, specifically its subtype or the age of the donor.
the cage tissue transcriptome from riken on the other hand, contained not only representative adult human submandibular glands but also parotid glands thus allowing a more equitable comparison.
upon hierarchical clustering of the processed cage data , the preservation of the mouse tissue specific signature was distinctly more robust, with  <dig> genes  showing cross-species conservation based on shared levels of salivary gland specific gene enrichment .
taken together, our analysis point to a conserved network of crucial genes and pathways that regulate common physiological processes in mouse and human sg and underscores the overall value and usefulness of mouse genetic models to study this tissue.fig.
<dig> preservation of the tissue specific salivary gland gene signature in human tissues.
a hierarchical clustering of human tissues using averaged fpkm values  of the genes from the tissue specific salivary gland gene signature with human-mouse homology.
the red colored dendrogram highlights the genes  that are preserved in the tissue specific enrichment observed in the mouse expression analysis.
b hierarchical clustering of human tissues using the cap analysis gene expression  data, represented by log <dig> transformed deseq <dig> median normalized tss tag counts using the same gene signature as in panel a, with  <dig> genes showing cross-species conservation
discussion
proper organ development requires remarkably intricate and complex biological processes that rely on spatially and temporally controlled gene expression programs.
this paradigm is quite apparent in the case of the submandibular salivary gland, which undergoes a dynamic process of morphogenesis during embryogenesis and further specification and functional maturation after birth.
in order to better understand the molecular underpinnings and identify novel and unique genes that may play important roles in salivary gland biology, here we have performed genome wide expression profiling of the mouse salivary gland during various stages of development and adult.
in addition to generating stage-specific gene expression profiles of the submandibular salivary gland, we have also utilized sophisticated computational tools and leveraged the large encode and fantom based data sets to obtain a systems biology driven perspective on this complex organ.
we have focused our studies on  <dig> different embryonic time points during which the salivary gland undergoes tremendous growth and morphological changes including branching morphogenesis, cyto-differentiation and lumenization.
as expected, these changes are associated with distinct alterations in the gene expression profile of critical regulatory factors and pathways, some of which are likely to be critical drivers of the developmental processes.
our analysis has led to the identification of a stage-specific and more broadly embryonic-specific gene expression profile, which can serve as a starting point in identifying new players that may be important in salivary gland development.
for example after mining our embryo enriched dataset, we identified several genes which showed peak expression levels very early during salivary gland morphogenesis and then declined over the course of late embryonic development and into adult.
while several of the highly differentially expressed genes at e <dig>  have been implicated in axonal guidance such as unc5d  <cit>  and tubb2b  <cit> , we also identified additional genes from this group that may be interesting for future follow up studies including scube <dig>  <cit>  and dlx <dig>  <cit> , both of which have been implicated in craniofacial development but not in salivary gland biology.
needless to say, rna-seq based analysis similar to those described in this work can be further extended to additional time points during mouse embryogenesis, in particular at stages earlier than e <dig>  during which the salivary gland fate is being established.
our rna-seq based experiments not only allowed us to generate a broad transcriptomic map of the postnatal salivary gland, but also yielded interesting insights into the continued salivary gland development and maturation processes that occur after birth.
indeed, our analysis demonstrated elevated expression levels of genes important in secretory cell differentiation during early postnatal development further highlighting the sustained cellular differentiation programs that continue after birth.
in contrast, at ~12 weeks the adult salivary gland is enriched for genes that serve valuable functional roles in this organ.
among these, we identified amy <dig> as one of the most highly expressed genes, which is in good agreement with the known function of this gene in salivary gland biology and oral health  <cit> .
we also observed specific enrichment of the solute carrier family of genes  which have been shown to be important in fluid and electrolyte secretion in the salivary gland  <cit> .
a major innovative aspect of our study is the development of the adult salivary gland gene signature – we posit that the unique combination of genes bestows an organ its tissue identity and can be extremely revelatory about organ development and function.
by leveraging encode data and using robust statistical analyses, we have demonstrated that at the global transcript levels, the adult salivary gland most closely resembles another exocrine organ, the pancreas as well as several other epithelial-rich organs including the bladder, skin and placenta.
although these global transcriptomic comparisons are per se quite insightful and suggest common developmental origins or functional similarities between organs, caveats such as differences in age, sex and the physiological state of the animals and the experimental conditions, is worth keeping in mind.
as an example, unlike the tight clustering of the pancreas and salivary gland, the mammary gland exhibited a slightly distant gene expression signature.
we suspect that this lack of similarity might be in part due to the fact that the secretory activity of the mammary gland is primarily limited to lactation – a physiological state that was not included in our organ comparisons.
our analysis also allowed us to generate a list of genes that are highly enriched for in the mouse salivary glands, a subset of which also show similar selective enrichment in the human salivary gland.
notably these salivary gland-enriched genes hold the key to crucial aspects of adult salivary gland biology as evident from gene ontology terms and thus offer a priority list of candidates that can be leveraged for salivary gland specific functions.
for example in a recent study, maruyama et al.
generated three inducible cre-recombinase mouse strains to examine the roles of various cell types in salivary gland homeostasis - it is interesting to note that all three genes examined in this study are included in our salivary gland gene signature analysis  <cit> .
moreover, the datasets generated by our rna-seq studies can serve as a nice complement to other resources such as the rapidly expanding emage  database of in situ gene expression in the mouse embryo  <cit> .
as an example, a cursory examination of emage showed strong expression of the transcription factor ehf in the submandibular gland primordium, in agreement with its high expression as revealed from our rna-seq analysis.
therefore, our salivary gland gene signature can be a valuable tool for future studies on cell/tissue type specific gene expression mechanisms or delivery tools for salivary gland biologists.
CONCLUSIONS
our work as described here is the most comprehensive systems-wide deep sequencing based transcriptomics study of the submandibular salivary gland to date and illustrates how such studies can shine new light into the rich diversity of genes and pathways that are likely to be functionally important for this organ.
knowing the gene expression profiles of the salivary gland at its major developmental time points can greatly increase our understanding of salivary gland biology, aid with disease state diagnosis, and help identify potential therapeutic targets for regeneration and tissue engineering approaches in the future.

 salivary gland morphogenesis requires the complex coordination of cells to orchestrate a number of dynamic cellular processes including cell specification, lineage commitment, cell migration, proliferation and differentiation, all culminating in the formation of this specialized gland. a network of signaling and regulatory molecules coordinates these vital biological processes, which are accompanied by the dynamic changes in gene expression throughout development. in the mouse, submandibular salivary gland morphogenesis occurs over several distinct developmental stages commencing at ~ embryonic day. at this early prebud stage, the primordial sg fate is established with the thickening of the adjoining oral epithelium. at the next bud stage, which occurs at approximately e, the thickened epithelium invaginates into the underlying mesenchyme resulting in the formation of a primary bud, which serves as the precursor of the main duct of the salivary gland once the gland reaches maturation. as the embryo develops to e the gland undergoes rapid proliferation and intricate branching morphogenesis, during which the end buds undergo successive rounds of clefting resulting in the generation of multiple epithelial buds. this pseudoglandular stage also coincides with reorganization of the end buds and the formation of the acini, which are the main secretory units of the salivary gland. at the canalicular stage, the gland is highly branched with lumenization of the main secretory duct nearing completion. the onset of gland cyto-differentiation also occurs at this stage – a process that continues to birth. in the final step of embryonic salivary gland morphogenesis at the terminal bud stage, expansion of acini and lumenization of both ducts and acini is almost finished, resulting in a continuous network of ducts connecting the acini to the oral cavity. while the gland is functional and able to secrete saliva at birth, further acinar maturation and differentiation continues postnatally, and by puberty, differentiation of the granular convoluted tubules is complete. for many years now, mouse genetic models have been widely utilized to study various facets of salivary gland biology including branching morphogenesis, cleft formation, organ development and differentiation. while these studies have been instrumental in identifying some of the individual genes and signaling pathways necessary for proper salivary gland function, a limited number of studies have focused on the global examination of salivary gland gene expression in mouse. in an early study, hoffman et al. utilized microarrays to examine gene expression profiles at five different stages of mouse submandibular salivary gland embryonic development. in a more recent and complementary study, musselmann et al. generated a microanatomical atlas of gene expression of embryonic salivary gland by performing laser capture microdissection of distinct epithelial populations obtained from e -e mouse embryos. these broad based approaches of global examination of gene expression profiles have been instrumental to deciphering molecular mechanisms of salivary gland morphogenesis and importantly in the discovery of novel signaling pathways such as fibroblast growth factor and bone morphogenetic protein as well as key signaling molecules such as gsk3β that play important roles in regulating branching morphogenesis in embryonic salivary gland. it is important to note that whereas gene profiling experiments during embryonic development of salivary gland has received much attention, similar in depth studies in adult has been lacking although one recent work has shed some light on aging associated sg gene signature. while microarray technology have been the application of choice in the past for transcriptome analysis, recent advancements have seen the supplanting of microarrays by genomic methods driven by next-generation sequencing approaches like rna-sequencing. compared to microarray technology, deep sequencing based methods provide a more sensitive and precise analytical approach that can accurately quantify gene transcript levels and their isoforms across a broad dynamic range. indeed, direct comparisons to rna-seq has revealed the shortcomings of hybridization-based microarray studies such as background noise and saturation and probe set issues such as incorrect annotation and isoform coverage. with the recent explosion in rna-seq based genome wide analysis of the transcriptomic landscapes of multiple cell types and various tissues and organs, the time is ripe to use computational tools and a systems biology based approach to gain novel insights into tissue-specific gene signatures and their potential impact on biological processes. in the present study, we have performed rna-seq to generate a comprehensive gene expression profile of the mouse submandibular gland at various time points of development during embryogenesis and maturation in adult. to the best of our knowledge this is the first reported rna-seq based study to examine the transcriptome of the mouse submandibular gland. an extensive bioinformatics analysis of our datasets has revealed interesting gene regulatory networks and maps that are enriched for and define the various stages of salivary glands. finally we have leveraged the encode, fantom, human protein atlas as well as other published gene expression datasets to identify a salivary gland specific gene signature that is to a large extent conserved between mouse and human. collectively our study not only validates existing literature but also provides a wealth of genomic resources that can be harnessed for the discovery of new genes and biologically important pathways in the salivary gland and for formulating testable hypotheses.  transcriptomic map of the mouse salivary gland during development and adult in order to better define the dynamic changes in global gene expression levels and to identify new tissue-specific and tissue-enriched regulators, we performed rna-seq based expression profiling of the mouse submandibular salivary gland at various key stages of embryonic development in addition to post-natal and adult tissues. towards this end, we isolated total rna from the mouse submandibular gland dissected from e  day old embryos, e and e. we also examined submandibular glands isolated from day, 4 week, and 12wk old male animals that represented different stages of post-natal maturation. biological replicates of these rna samples were sequenced to high depth using an illumina hiseq platform. the sequence reads were mapped to the reference genome sequence of mus musculus using tophat. we subsequently performed between-sample normalization using the deseq median normalization method and calculated fragments per kilobase of transcripts per million mapped reads thereby giving us measurements of relative expression of genes within and between biological samples.fig. principal component analysis of the mouse salivary glands at different developmental time points. a experimental scheme. we isolated total rna from whole salivary glands ranging from early embryo to adult, and performed rna-seq. utilizing these datasets, we defined and annotated the salivary gland transcriptional landscape by using various gene ontology annotation analyses and pathway analyses. b proportion of variance in each principle component. pc pc and pc represent ~90% of variance in the data. c projection plots show the pca coordinates for each stage, which are indicated by different colors. the data indicates that the inherent variations in gene expression between biological samples can distinguish the developing salivary gland in a time dependent manner in order to better analyze and appreciate the overall gene expression patterns between the various developmental and adult time points, we utilized principal component analysis, a statistical technique that reduces and summarizes large datasets while illustrating relationships between samples based on co-variance of the data being examined. using pca, we found that pc pc and pc accounted for approximately 90% of all variations of the original data. in order to further explore and better depict the major sources of variation, all samples were plotted in a three-dimensional space consisting of pc pc and pc interestingly, as demonstrated in fig. 1c, each of the representative time points datasets segregated into individual groups demonstrating the highly dynamic variation in gene expression between each sg sample. indeed, biological replicates cluster tightly together, further highlighting the inherent similarity of the biological samples to one another. another notable observation was that the embryonic samples clustered more closely to each other and that there is a striking separation between these and the adult samples. taken together this analysis provided the first hint of a clear dichotomy of gene expression profiles between embryonic and postnatal salivary gland samples. rna-seq analysis identifies a salivary gland gene signature to evaluate development-dependent differential gene expression patterns, we next grouped our samples based on distinct developmental stages corresponding to embryogenesis, neonatal and adult. for this analysis, we selected genes that showed at least a two-fold change in expression between each time point while showing an adjusted p-value of less than. we also considered genes that were expressed at ≥ fpkm in at least one biological replicate. using this criterion, we identified a total of differentially expressed genes between the embryonic and neonatal stages with genes found to be up-regulated and genes showing downregulation. similar comparison of degs between neonatal and adult developmental time points identified genes. of these degs, were upregulated and downregulated. finally, we found degs between embryo and adult with number of genes showing up-regulation and downregulated. to better appreciate the biological relevance of the global transcriptomic differences between the embryonic and adult mouse salivary gland, we further analyzed the degs using clusterprofiler and identified pathways unique to each state. interestingly, in all embryo enriched data sets, we observed specific enrichment of biological processes important during organogenesis, such as axon guidance as well as that of wnt and notch signaling pathways, both of which have been shown to be important for proper salivary gland morphogenesis. in contrast, in the adult enriched data sets we observed enrichment in biological processes that include salivary secretion, protein export, fatty acid metabolism and vesicular signaling. notably all of these processes have been shown to be critical for proper salivary gland function and justifiably represent the mature stage of the adult gland. armed with a global view of the overall transcriptional changes occurring between the developmental stages, we next sought to identify a developmental stage-specific global gene expression profile of the salivary gland. we reasoned that this broad gene profile might be useful in providing new insights into the biological processes and pathways that are important to each developmental stage. to generate the developmental specific gene expression profile, we focused on genes that a) showed an fpkm of or greater in at least one of the sample replicates across all samples and b) were comparatively differentially expressed, with a log2-fold change of or greater, at one time point compared to other time points across all datasets. based on the above criteria we identified a total of genes that were specifically enriched, across all datasets as demonstrated in fig. 2a. to identify unique patterns of transcriptomic activities of the developmental stage-specific gene profile, we performed k-means clustering on the standardized log transformed fpkm values for each of the samples. this analysis led to the identification of different gene clusters that corresponded to genes enriched in a specific developmental time point as well as genes enriched in grouped time points, which we refer to as embryo, neonatal and adult enriched clusters. interestingly, our analysis revealed that genes that were uniquely expressed during embryogenesis showed a gradual decline in expression as development progressed, while genes enriched in the neonatal stage showed a gradual increase in expression during embryogenesis with peak expression levels observed at the neonatal stage and then showing a gradual decline. conversely, genes enriched in the adult stage showed a progressive increase during the developmental and neonatal stages with the highest levels of expression in the adult stage.fig. cluster analysis of the salivary gland developmental profile. a z-scores for each of the genes were calculated and used as input for gene-wise k-means clustering analysis. this analysis, visualized as a heatmap, depicts enrichment of genes at both specific development time points and general developmental stages. b visualization of the gene z-scores from the general developmental stage enriched clusters depicts unique time-dependent patterns of expression. in general, this analysis has identified genes that have peak expression in the defined developmental stages. also show are selected enriched go-annotations for each developmental stage gene cluster gene ontology analysis of the developmental stage-specific gene profile identified overrepresentation of genes associated with specific functions, indicating the existence of distinct biological processes occurring during embryonic, neonatal and adult stages of salivary gland development. a closer examination of genes that were uniquely expressed during embryogenesis revealed specific enrichment of biological pathways associated with organ, embryonic and vasculature development, cell differentiation, and transcription - all of these biological processes are associated with organ morphogenesis. interestingly, we also observed specific enrichment in both nervous system development and nerve growth factor receptor signaling in salivary glands during embryogenesis. this is in good agreement with studies demonstrating the critical role of the parasympathetic nervous system during salivary gland organogenesis. genes specifically enriched in the neonatal stage also revealed association with interesting biological terms including regulation of secretion, lipid metabolic processes, catabolic processes, chloride transport, protein modification, and vesicular transport - all of which have been shown to play a role in secretory cell differentiation. finally, analysis of the adult stage identified enriched biological pathways associated with salivary gland specific functions including saliva secretion, digestion, and ion and carbohydrate transport. we also identified smooth muscle contraction as an enriched biological pathway in the adult sg, which fits very well with the existence of myoepithelial cell populations that are important for glandular contraction and expulsion of saliva. taken together, the go analysis, confirms the validity and the possible functional significance of our sg developmental stage-specific gene expression profile. since our analysis of the broad gene-enriched profile of embryonic, neonatal and adult salivary gland revealed interesting facets of its biology, we wondered if further in-depth parsing of the datasets could provide more detailed insights into salivary gland maturation, particularly in the embryonic stages. hence, we next evaluated the transcriptomic activities from e to e to identify genes and pathways that may be critical during the early stages of salivary gland organogenesis and development. of the comparatively differentially expressed developmental specific genes we identified, of them were specifically enriched during embryogenesis. a closer analysis revealed that the greatest number of differentially expressed genes occurred at e, at a stage when the gland is undergoing rapid proliferation and has embarked upon a dynamic and elaborate program of branching morphogenesis. while the e functional analysis demonstrated broad developmental processes including salivary gland morphogenesis, exocrine system development, epithelium development, axonogenesis and axon guidance, we also observed selective enrichment in biological processes specific to the developmental programs that are at play at this critical time point during morphogenesis. for instance, in agreement with the burgeoning branching morphogenesis that takes place at e, genes involved in pathways specific to epithelial development, branching morphogenesis, and salivary gland cavitation were overrepresented. indeed, both branching and cavitation are critical for the proper development of the tubular network necessary for the transport of saliva to the oral cavity. while some degree of overlap was observed between e and e in processes involving cell differentiation and organ development, we did identify unique biological functions related to muscle function at e. interestingly, we found specific enrichment in pathways involving muscle structure, muscle tissue development and regulation of muscle contraction, properties intimately associated with functions of the myoepithelial cells. these findings are in good agreement with a previous study demonstrating the emergence of a myoepithelial cell population in the salivary gland at e. finally, our analysis at e uncovered several biological pathways enriched for a variety of salivary gland specific functions including metabolic processes. more specifically we observed over representation of metabolic processes important in amino acid metabolism and small molecule metabolism. we also identified pathways involved in transport and proteolysis all of which have been shown to play roles in salivary gland biology and secretion. taken together, our analysis suggests that the underlying global transcriptomic profile mirrors the distinct developmental stages of the embryonic sg as highlighted by branching morphogenesis, myoepithelial development and secretory cell development.fig. enriched biological process networks during embryogenesis. a heatmap visualization of the z-scores from the genes identified in the embryo specific clusters generated from the analysis in fig.  b–d) network visualization of selected top enriched biological processes at e, e and e. the networks were assembled by cytoscape tool enrichmap, using an organic layout. the node size represents the number of genes assigned to a biological process and edge width is proportional to the number of overlapping genes between two nodes meta-analysis of mouse tissue rna-seq expression data our salivary gland rna-seq datasets offered a wealth of in-depth information regarding the transcriptomic repertoire of this tissue. we rationalized that although many of the genes and pathways are common to sg and numerous other tissues and organs, there might exist some genes/pathways that are enriched for in expression and possibly function in the sg. to probe this further, we examined rna-seq data for adult organs and tissues generated by the mouse encode consortium, and compared the overall gene expression profiles with those from the adult sg described in this study. for this purpose, we generated a pairwise pearson correlation matrix using the top genes with the highest median absolute deviation. this, we reasoned, would allow us to identify genes with the most dynamic range of expression across all tissues. our analysis revealed a clear separation of organs based on gene expression with the salivary gland clustering more closely with exocrine organs such as the pancreas and organs rich in stratified epithelial tissues including the skin, bladder and placenta. the organ-specific segregation patterns were not a result of using a stringent cutoff of the top genes, since a pairwise pearson correlation matrix using a larger set of genes, yielded similar results with the salivary gland continuing to cluster closely with the pancreas, skin, bladder and placenta. a similar hierarchical clustering analysis also revealed a select list of genes that are similarly expressed between the salivary gland and the closely related four organs identified above and those that are exclusive to the two exocrine glands; the salivary gland and pancreas. interestingly, the transcriptomic identity of the mouse salivary gland was quite distinct from the mammary gland, another glandular tissue. we suspect this is due to the fact the secretory activity of the mammary gland is limited only during parity and lactation rather than the virgin state, from which the rna-seq data was obtained. in contrast, tissues representing specific brain regions such as the cerebellum and the mouse testis were surprisingly clustered close to the group containing the salivary gland and the four aforementioned organs – the significance of this finding is unclear and worthy of future exploration. taken together, our findings demonstrate that at the broader global transcriptomic level, the salivary gland most closely resembles tissues and organs with similar morphological and functional characteristics.fig. hierarchical clustering of mouse tissues. fpkm values from the top genes with the highest median absolute deviation were used to cluster adult mouse tissues. the resulting heatmap shows that the salivary gland clusters closely with the pancreas, skin, bladder and placenta having obtained a global gene expression profile of the salivary gland, we next wondered if we could derive an adult salivary gland specific molecular signature that is constituted by genes that show a distinct expression profile and hence are likely to be more relevant and important in salivary gland biology. towards this end, we next compared gene expression levels across all adult mouse organs and tissues and identified genes that were specifically enriched in the salivary gland and thus represented a potential salivary gland gene signature. to gain a better understanding of the underlying biological functions and pathways associated with the adult salivary gland gene signature, we utilized cluego, which is a functional gene ontology analysis tool that integrates several gene-set enrichment databases, including the kyoto encyclopedia of genes and genomes, reactome pathway database annotations, and the go consortium database, to create a comprehensive go/pathway term network. this analysis not only highlighted the close association of the signature to exocrine and salivary gland development but also identified enrichment in biological processes related specifically to protein export, protein processing and localization to the endoplasmic reticulum, all of which are quite relevant for proper salivary gland function. interestingly, we also observed pathway enrichment specific to mucin biosynthesis – this is in good agreement with the established roles of mucins in salivary gland biology that includes their role as a first line of defense against microbial infection and their important contributions to various biophysical properties of saliva such as viscoelasticity.fig. visualization of the tissue specific salivary gland gene signature. a hierarchical clustering of the gene expression values selected from the tissue specific salivary gland signature across adult mouse tissues. b network visualization of enriched pathways in the gene signature was performed by cluego analysis we next sought to identify a cohort of transcription factors, which are selectively enriched in our sg gene signature. based on the crucial regulatory role transcription factors play in driving tissue-specific gene expression, we reasoned that this analysis may aid in identifying novel transcriptional regulators which may be important in salivary gland biology. to ensure robust coverage, we first compiled a comprehensive mouse transcription factor list using databases from both riken and uniprotkb and then searched for transcription factors which are represented in our mouse adult salivary gland gene signature. our analysis uncovered transcriptional regulators. of these transcription factors, ascl, bhlha, tfcp2l and six have been previously reported to play a role in salivary gland development and function. our analysis also unearthed elf and ehf, two members of the ets family of transcription factors, which are highly expressed in the salivary gland yet their specific function, if any, in salivary gland biology remains unknown. moreover, we identified several other factors including eaf and foxi neither of which has been previously studied in the context of salivary gland biology. the high expression of some of these transcription factors in salivary glands was also confirmed using publicly available datasets provided through the salivary gland molecular anatomy project. in addition to the list of transcriptional regulators that were specifically enriched in our salivary gland gene signature, there were several other genes and gene family members that are both over-represented and likely or known to be associated with salivary gland biology. this included genes belonging to the secretoglobin gene superfamily, which encode for small secretory proteins found in high concentrations in secretions from various organs including the lacrimal and salivary glands. in addition, we observed selective enrichment in genes belonging to the kallikrein family that encode for proteins which function as serine peptidases. while a major function of this family of proteins have been associated with blood pressure regulation and skin homeostasis, some studies suggest that these enzymes may play a role in cell proliferation, cell survival and wound healing. finally, we also observed elevated levels of genes belonging to the demilune cell and parotid protein family, which are a family of genes thought to possess antimicrobial activity. we posit that many of these genes unearthed by our systems biology based approach are likely to be strong candidates for various physiological functions associated with mouse salivary glands.fig. visualization of selected members of the tissue specific salivary gland gene signature. a hierarchical cluster analysis of the enriched transcription factors identified in the gene signature. b hierarchical cluster analysis of enriched genes identified in the signature and which have been shown to play a role in salivary gland biology evolutionarily conserved expression map between human and mouse salivary gland mouse models have served a valuable role in elucidating molecular mechanisms of physiological and pathological processes in the salivary gland. it is likely that the genetic interplays and gene regulatory networks that underlie these processes are mostly conserved between mouse and human. to examine this we next asked if there were any commonalities in the salivary gland specific gene signature between mouse and human. towards this end, we mined the publicly available transcriptome database for human tissues and organs generated by the human protein atlas project and the fantom consortium ). to make proper comparisons, we only chose representative human tissue samples for which the corresponding data was available for the mouse tissues. upon examination of the human protein atlas samples, we found that of the human genes that were homologous to the mouse, genes showed enrichment in the human salivary gland as compared to the other tissues analyzed, thus recapitulating the tissue-restricted patterns of expression observed from our mouse tissue analysis. one possible caveat to this comparison is that for the hpa dataset, there is limited information about the salivary gland, specifically its subtype or the age of the donor. the cage tissue transcriptome from riken on the other hand, contained not only representative adult human submandibular glands but also parotid glands thus allowing a more equitable comparison. upon hierarchical clustering of the processed cage data, the preservation of the mouse tissue specific signature was distinctly more robust, with genes showing cross-species conservation based on shared levels of salivary gland specific gene enrichment. taken together, our analysis point to a conserved network of crucial genes and pathways that regulate common physiological processes in mouse and human sg and underscores the overall value and usefulness of mouse genetic models to study this tissue.fig. preservation of the tissue specific salivary gland gene signature in human tissues. a hierarchical clustering of human tissues using averaged fpkm values of the genes from the tissue specific salivary gland gene signature with human-mouse homology. the red colored dendrogram highlights the genes that are preserved in the tissue specific enrichment observed in the mouse expression analysis. b hierarchical clustering of human tissues using the cap analysis gene expression data, represented by log transformed deseq median normalized tss tag counts using the same gene signature as in panel a, with genes showing cross-species conservation discussion proper organ development requires remarkably intricate and complex biological processes that rely on spatially and temporally controlled gene expression programs. this paradigm is quite apparent in the case of the submandibular salivary gland, which undergoes a dynamic process of morphogenesis during embryogenesis and further specification and functional maturation after birth. in order to better understand the molecular underpinnings and identify novel and unique genes that may play important roles in salivary gland biology, here we have performed genome wide expression profiling of the mouse salivary gland during various stages of development and adult. in addition to generating stage-specific gene expression profiles of the submandibular salivary gland, we have also utilized sophisticated computational tools and leveraged the large encode and fantom based data sets to obtain a systems biology driven perspective on this complex organ. we have focused our studies on different embryonic time points during which the salivary gland undergoes tremendous growth and morphological changes including branching morphogenesis, cyto-differentiation and lumenization. as expected, these changes are associated with distinct alterations in the gene expression profile of critical regulatory factors and pathways, some of which are likely to be critical drivers of the developmental processes. our analysis has led to the identification of a stage-specific and more broadly embryonic-specific gene expression profile, which can serve as a starting point in identifying new players that may be important in salivary gland development. for example after mining our embryo enriched dataset, we identified several genes which showed peak expression levels very early during salivary gland morphogenesis and then declined over the course of late embryonic development and into adult. while several of the highly differentially expressed genes at e have been implicated in axonal guidance such as unc5d and tubb2b, we also identified additional genes from this group that may be interesting for future follow up studies including scube and dlx, both of which have been implicated in craniofacial development but not in salivary gland biology. needless to say, rna-seq based analysis similar to those described in this work can be further extended to additional time points during mouse embryogenesis, in particular at stages earlier than e during which the salivary gland fate is being established. our rna-seq based experiments not only allowed us to generate a broad transcriptomic map of the postnatal salivary gland, but also yielded interesting insights into the continued salivary gland development and maturation processes that occur after birth. indeed, our analysis demonstrated elevated expression levels of genes important in secretory cell differentiation during early postnatal development further highlighting the sustained cellular differentiation programs that continue after birth. in contrast, at ~12 weeks the adult salivary gland is enriched for genes that serve valuable functional roles in this organ. among these, we identified amy as one of the most highly expressed genes, which is in good agreement with the known function of this gene in salivary gland biology and oral health. we also observed specific enrichment of the solute carrier family of genes which have been shown to be important in fluid and electrolyte secretion in the salivary gland. a major innovative aspect of our study is the development of the adult salivary gland gene signature – we posit that the unique combination of genes bestows an organ its tissue identity and can be extremely revelatory about organ development and function. by leveraging encode data and using robust statistical analyses, we have demonstrated that at the global transcript levels, the adult salivary gland most closely resembles another exocrine organ, the pancreas as well as several other epithelial-rich organs including the bladder, skin and placenta. although these global transcriptomic comparisons are per se quite insightful and suggest common developmental origins or functional similarities between organs, caveats such as differences in age, sex and the physiological state of the animals and the experimental conditions, is worth keeping in mind. as an example, unlike the tight clustering of the pancreas and salivary gland, the mammary gland exhibited a slightly distant gene expression signature. we suspect that this lack of similarity might be in part due to the fact that the secretory activity of the mammary gland is primarily limited to lactation – a physiological state that was not included in our organ comparisons. our analysis also allowed us to generate a list of genes that are highly enriched for in the mouse salivary glands, a subset of which also show similar selective enrichment in the human salivary gland. notably these salivary gland-enriched genes hold the key to crucial aspects of adult salivary gland biology as evident from gene ontology terms and thus offer a priority list of candidates that can be leveraged for salivary gland specific functions. for example in a recent study, maruyama et al. generated three inducible cre-recombinase mouse strains to examine the roles of various cell types in salivary gland homeostasis - it is interesting to note that all three genes examined in this study are included in our salivary gland gene signature analysis. moreover, the datasets generated by our rna-seq studies can serve as a nice complement to other resources such as the rapidly expanding emage database of in situ gene expression in the mouse embryo. as an example, a cursory examination of emage showed strong expression of the transcription factor ehf in the submandibular gland primordium, in agreement with its high expression as revealed from our rna-seq analysis. therefore, our salivary gland gene signature can be a valuable tool for future studies on cell/tissue type specific gene expression mechanisms or delivery tools for salivary gland biologists.  our work as described here is the most comprehensive systems-wide deep sequencing based transcriptomics study of the submandibular salivary gland to date and illustrates how such studies can shine new light into the rich diversity of genes and pathways that are likely to be functionally important for this organ. knowing the gene expression profiles of the salivary gland at its major developmental time points can greatly increase our understanding of salivary gland biology, aid with disease state diagnosis, and help identify potential therapeutic targets for regeneration and tissue engineering approaches in the future. 
to better appreciate the wide spectrum of gene expression profiles, we isolated rna from mouse submandibular salivary glands at different embryonic and adult stages. compared to microarray, rna-sequencing offers unbiased detection of novel transcripts, broader dynamic range and high specificity and sensitivity for detection of genes, transcripts, and differential gene expression. these rna-seq datasets will complement pre-existing microarray based datasets, including the salivary gland molecular anatomy project by offering a broader systems-biology based perspective rather than the classical gene-centric view. mouse models have served a valuable role in deciphering various facets of salivary gland biology, from normal developmental programs to diseased states. to identify molecular players and pathways likely to be relevant for sg biology, we performed functional gene enrichment analysis, network construction and hierarchal clustering of the rna-seq datasets obtained from different stages of sg development and maturation, and other mouse organs and tissues. to facilitate such studies, gene expression profiling maps have been generated for various stages of sg organogenesis. our rna-seq based atlas has revealed a high-resolution cartographic view of the dynamic transcriptomic landscape of the mouse sg at various stages. in parallel, we processed rna-seq data for organs and tissues obtained from the mouse encode consortium and calculated the average gene expression values. although rna-seq data, particularly under the auspices of the encode project, have covered a large number of biological specimens, studies on the sg have been lacking.
 mouse models have served a valuable role in deciphering various facets of salivary gland biology, from normal developmental programs to diseased states. to facilitate such studies, gene expression profiling maps have been generated for various stages of sg organogenesis. however these prior studies fall short of capturing the transcriptional complexity due to the limited scope of gene-centric microarray-based technology. compared to microarray, rna-sequencing offers unbiased detection of novel transcripts, broader dynamic range and high specificity and sensitivity for detection of genes, transcripts, and differential gene expression. although rna-seq data, particularly under the auspices of the encode project, have covered a large number of biological specimens, studies on the sg have been lacking.  to better appreciate the wide spectrum of gene expression profiles, we isolated rna from mouse submandibular salivary glands at different embryonic and adult stages. in parallel, we processed rna-seq data for organs and tissues obtained from the mouse encode consortium and calculated the average gene expression values. to identify molecular players and pathways likely to be relevant for sg biology, we performed functional gene enrichment analysis, network construction and hierarchal clustering of the rna-seq datasets obtained from different stages of sg development and maturation, and other mouse organs and tissues. our bioinformatics-based data analysis not only reaffirmed known modulators of sg morphogenesis but revealed novel transcription factors and signaling pathways unique to mouse sg biology and function. finally we demonstrated that the unique sg gene signature obtained from our mouse studies is also well conserved and can demarcate features of the human sg transcriptome that is different from other tissues.  our rna-seq based atlas has revealed a high-resolution cartographic view of the dynamic transcriptomic landscape of the mouse sg at various stages. these rna-seq datasets will complement pre-existing microarray based datasets, including the salivary gland molecular anatomy project by offering a broader systems-biology based perspective rather than the classical gene-centric view. ultimately such resources will be valuable in providing a useful toolkit to better understand how the diverse cell population of the sg are organized and controlled during development and differentiation. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords salivary glandrna-sequencinggene signaturedevelopmentstate university of new york at buffaloschool of dental medicine start-up fundsromano rose-anne issue-copyright-statement© the author 2016 