System-level analyses of keystone genes required for mammalian tooth development

Abstract

When a null mutation of a gene causes a complete developmental arrest, the gene is typically considered essential for life. Yet, in most cases, null mutations have more subtle effects on the phenotype. Here we used the phenotypic severity of mutations as a tool to examine system-level dynamics of gene expression. We classify genes required for the normal development of the mouse molar into different categories that range from essential to subtle modification of the phenotype. Collectively, we call these the developmental keystone genes. Transcriptome profiling using microarray and RNAseq analyses of patterning stage mouse molars show highly elevated expression levels for genes essential for the progression of tooth development, a result reminiscent of essential genes in single-cell organisms. Elevated expression levels of progression genes were also detected in developing rat molars, suggesting evolutionary conservation of this system-level dynamics. Single-cell RNAseq analyses of developing mouse molars reveal that even though the size of the expression domain, measured in the number of cells, is the main driver of organ-level expression, progression genes show high cell-level transcript abundances. Progression genes are also upregulated within their pathways, which themselves are highly expressed. In contrast, a high proportion of the genes required for normal tooth patterning are secreted ligands that are expressed in fewer cells than their receptors and intracellular components. Overall, even though expression patterns of individual genes can be highly different, conserved system-level principles of gene expression can be detected using phenotypically defined gene categories.

Keywords: essential genes; keystone genes; single-cell RNAseq; tooth development; transcript abundance; transcriptomes.

Figure 1

Keystone gene categories of tooth development. Mouse molar development progresses from initiation and patterning to formation of the hard tissues and eruption. These steps are mediated by reciprocal signaling between epithelium (pink) and mesenchyme (magenta). A central step in the patterning is the formation of the epithelial signaling center, the primary enamel knot (blue oval inside the cap stage tooth). Several genes are known to be required for the developmental progression and regulation of the shape around the time of cap stage, and here we focused mainly on transcriptomes in the bud, and cap stage molars. Expression of progression and shape category genes were compared to tissue and dispensable genes, as also to other developmental process genes. Fewer initiation and eruption category genes are known, and they were excluded from the analyses. For listing of the genes, see Appendix S1 and Table S1 [Color figure can be viewed at wileyonlinelibrary.com]

Figure 2

Bud and cap stage mouse molars show elevated expression of progression and shape genes. (a) Microarray and (b) RNAseq analyses of the transcriptomes of E13 and E14 molars show highly elevated expression of progression category and moderately elevated shape category genes (for tests, see Table S3). Tissue category genes are involved in dentine and enamel formation that begin at birth, around 6 days from the E14 cap stage. The number of genes having RNAseq expression data in each category are 15, 28, 27, 100, 4106, and 16,165 for progression, shape, tissue, dispensable, developmental process, and other, respectively. The corresponding numbers for the microarray data are 15, 28, 27, 98, 3983, and 14,825. Boxes enclose 50% of observations; the median and mean are indicated with a horizontal bar and diamond, respectively, and whiskers extend to last values within 1.5 interquartiles. Individual data points are shown for the smaller categories [Color figure can be viewed at wileyonlinelibrary.com]

Figure 3

Bud and cap stage rat molars show elevated expression of keystone genes. RNAseq analyses of the transcriptomes of E15 and E17 rat molars show highly elevated expression of progression category and moderately elevated shape category genes (for tests, see Table S3). The numbers of genes having RNAseq expression data in each category are 15, 28, 27, 95, 3843, and 12,473 for progression, shape, tissue, dispensable, developmental process, and other, respectively. Boxes enclose 50% of observations; the median and mean are indicated with a horizontal bar and diamond, respectively, and whiskers extend to last values within 1.5 interquartiles. Individual data points are shown for the smaller categories [Color figure can be viewed at wileyonlinelibrary.com]

Figure 4

Single‐cell RNAseq data reflect the bulkRNA analyses. (a) Bulk and scRNAseq expression levels show overall correspondence in tissue level expression in E14 mouse molar. (b) Progression category genes show the strongest upregulation whereas shape category is intermediate between the progression and other categories (for tests, see Table S3). The number of genes having expression data in each category are 15, 28, 25, 99, 3771, and 16,362 for progression, shape, tissue, dispensable, developmental process, and other, respectively. Boxes in (b) enclose 50% of observations; the median and mean are indicated with a horizontal bar and diamond, respectively, and whiskers extend to last values within 1.5 interquartiles. Individual data points are shown for the smaller categories [Color figure can be viewed at wileyonlinelibrary.com]

Figure 5

Expression domain size and high cell‐level transcript abundance of progression genes. (a) The size of the expression domain, measured as the number of cells, is the key driver of organ level expression (scRNAseq on E14 mouse molar, plotted for progression, shape, tissue, and dispensable category genes). (b) Progression category shows high cell‐level expression or transcript abundance, indicating high expression relative to the expression domain size. The p values for progression, shape, and tissue categories compared to dispensable category are .0008, .7009, .3413, respectively (one‐tailed significance levels obtained using 10,000 permutations). (c) tSNE plots showing diverse expression patterns of three different progression category genes in the scRNAseq data. Clusters containing epithelial and mesenchymal cells are marked (note the limited presence of mesenchymal marker Msx1 transcripts in the epithelium, agreeing with previous reports by Coudert et al., 2005). Boxes in (b) enclose 50% of observations; the median and mean are indicated with a horizontal bar and diamond, respectively, and whiskers extend to last values within 1.5 interquartiles [Color figure can be viewed at wileyonlinelibrary.com]

Figure 6

Restricted expression of ligands explains differences in upregulation. (a) Compared to all the developmental process genes, the pathways containing progression and shape category genes are generally upregulated (p < .0001, black line, random sampling using median expression levels and 10,000 permutations). (b) Median expression of progression genes (black line) shows further upregulation compared to all the genes in the corresponding pathways (p = .0004). (c) Median expression of shape category genes (black line) shows comparable expression with the genes in the corresponding pathways (p = .5919), but excluding ligands makes the shape category genes highly expressed (p = .0154, black dashed line). (d) Of the progression and shape category genes, ligands are expressed in fewer number of cells compared to the non‐secreted proteins. (e) Ligands are also expressed in relatively few cells among all the pathway genes. Analyses using E14 mouse bulk RNAseq (a–c) and E14 scRNAseq (d, e). The p value for ligands in (d) and (e) are p < .0001 (random sampling compared to all the other gene types). Boxes in (d, e) enclose 50% of observations; the median and mean are indicated with a horizontal bar and diamond, respectively, and whiskers extend to last values within 1.5 interquartiles