Mapping the global mRNA transcriptome during development of the murine first molar

The main objective of this study was to map global gene expression in order to provide information about the populations of mRNA species participating in murine tooth development at 24 h intervals, starting at the 11th embryonic day (E11.5) up to the 7th post-natal day (P7). The levels of RNA species expressed during murine tooth development were mesured using a total of 58 deoxyoligonucleotide microarrays. Microarray data was validated using real-time RT-PCR. Differentially expressed genes (p < 0.05) were subjected to bioinformatic analysis to identify cellular activities significantly associated with these genes. Using ANOVA the microarray data yielded 4362 genes as being differentially expressed from the 11th embryonic day (E11.5) up to 7 days post-natal (P7), 1921 of these being genes without known functions. The remaining 2441 genes were subjected to further statistical analysis using a supervised procedure. Bioinformatic analysis results for each time-point studied suggests that the main molecular functions associated with genes expressed at the early pre-natal stages (E12.5–E18.5) were cell cycle progression, cell morphology, lipid metabolism, cellular growth, proliferation, senescence and apoptosis, whereas most genes expressed at post-natal and secretory stages (P0–P7) were significantly associated with regulation of cell migration, biosynthesis, differentiation, oxidative stress, polarization and cell death. Differentially expressed genes (DE) not described earlier during murine tooth development; Inositol 1, 4, 5-triphosphate receptor 3 (Itpr3), metallothionein 1(Mt1), cyclin-dependent kinase 4 (Cdk4), cathepsin D (Ctsd), keratin complex 2, basic, gene 6a (Krt2-6a), cofilin 1, non-muscle (Cfl1), cyclin 2 (Ccnd2), were verified by real-time RT-PCR.


INTRODUCTION
During murine tooth development substantial changes occur within a time-span of 24 h. From E11.5 up to P7, i.e., in the course of 16 days the tooth germ develops from one layered oral epithelium into a phenotypic molar tooth. The developing murine tooth germ is therefore an excellent model for studying the time-course of gene expression in a rapidly, developing organ.
Differentiation requires participation of a large number of genes, the expression of which is regulated in time and space. The numbers and types of genes involved, however, vary depending on the developmental stage. In addition, miR-NAs and lncRNAs are also known to be involved during tooth development (Jevnaker and Osmundsen, 2008;Michon et al., 2010;Michon, 2011) providing an additional layer of complexity.
Microarray studies generate vast amounts of gene expression data. Genes that exhibit significantly different levels of expression at different developmental stages during tooth development are termed differentially expressed (DE) genes; these must be isolated using appropriate statistical methods (Reiner et al., 2003;Smyth, 2004). Secondly it is of interest to determine biological processes, biochemical functions, and sub-cellular locations significantly associated with DE genes, linking DE genes to alterations in cellular physiology by mapping DE genes to gene ontology (GO) terms (Ashburner et al., 2000;Rhee et al., 2008).
More than 4300 mRNAs and some of their translated proteins have, so far, been detected during tooth germ development by microarray, in situ hybridization and immunocytochemistry (Jevnaker and Osmundsen, 2008;Landin et al., 2012). In a previous study (Landin et al., 2012) we used three pre-selected genes (Ambn, Amelx, and Enam) as starting point for profile search and identified 84 differentially expressed genes with a similar expression pattern as these enamel genes. However, the results of this study only show a tinny frame of the big picture of genetic events that occur during murine tooth development. Mapping global gene expression to capture the majority of genes involved during murine tooth development may provide an overview of the occurring genetic changes.
The global mapping of gene expression for each timepoint studied may also reveal participation of novel genes or transcription factors as well as genes with unknown functions (Etokebe et al., 2009) during murine tooth development. Understanding the molecular cell biology during murine tooth development opens for development in novel bio-therapeutic strategies in dentistry.
In the current study we attempted to map the global gene expression in the molar murine tooth germ at each of 16 timepoints by uploading the 2441 genes differentially (P ≤ 0.05) expressed at every time point studied (E11.5-P7). To interpret the resulting gene-expression data, we used Ingenuity Pathway analysis (IPA) (Kramer et al., 2014).

EXPERIMENTAL DESIGN
The global gene expression in tooth germs from wild type mice was monitored from embryonic day 11 (E11.5) up to 7 days after birth (P7) using a reference design; the 11th embryonic day was considered time point zero (T 0 ) and all time-points studied were compared to E11.5. The sample size for each time point was n = 3-5 embryos/pups from three different mothers.

EXPERIMENTAL ANIMALS
Pregnant Balb-c mice CD-1strain were used in this study. The day of vaginal plug was set to E0.5 Adult mice were sacrificed by cervical dislocation, the embryos or pups by decapitation. Embryos were staged according to the Theiler criteria as described in Landin et al. (2012). Animal housing (Scantainer ventilated cabinet Q-110) had 12 h light/dark cycle. The cabinet temperature was maintained at 21 • C with a relative humidity of 55% (ScanClime plus). Fodder and water were supplied ad lib. The animals were kept according to the regulations of the Norwegian Gene Technology Act of 1994.

DISSECTION OF TOOTH GERMS AND RNA EXTRACTION
Dissection, homogenization of whole tooth germs and total RNA extraction was carried out as previously described in Osmundsen et al. (2007) and Landin et al. (2012). At the pre-natal stages total RNA was isolated from 3 to 5 tooth germs. At post-natal stages batches of at least three tooth germs were used at each time point. RNA concentration was measured at 260 nm in a Nanodrop ND 1000 spectrophotometer. RNA fractions with the ratio of absorbance 260 and 280 nm around 2.0 and with RINvalues higher than 8.5 as measured using an Agilent Bioanalyzer (Agilent, Palo Alto, CA, USA) were used for analysis of gene expression using deoxyoligonucleotide microarrays and real-time RT-PCR.

COMPLEMENTARY DNA SYNTHESIS AND LABELING
Complementary DNA (cDNA) was synthetized and indirect labeled from 1 μg total RNA using Genisphere Array 900™. Indirect labeling was used to avoid bias associated with differences in molecular size of the indicator molecules as previously described (Osmundsen et al., 2007;Landin et al., 2012).

MICROARRAY ANALYSIS OF mRNAs ISOLATED FROM TOOTH GERMS
Murine deoxyoligonucleotide (30 k)-microarrays were purchased from the NTNU Microarray Core Facility, Trondheim, Norway. The slides had been printed using the Operon murine v.3 oligo set (Qiagen GmbH, Hilden Germany). The arrays included probes for 10 different mRNAs from Arabidopsis thaliana. A spike mixture of 10 different mRNAs from A. thaliana (purchased from Stratagene, La Jolla, CA, USA) mixed in pairs at 10 different ratios, ranging from 0.1 to 5 was used to monitor the quality of the hybridisation. Each time point 3-6 biological replicates were subjected to independent microarray analysis. Each microarray was scanned in a Packard Bioscience Scanarray Lite microarray scanner (Perkin-Elmer). The Cy3 (543 nm) and Cy5 (645 nm) fluorescence signals were quantified by using the ScanArrayExpress v.3.0 software (Perkin-Elmer).

ACCESSION NUMBER
Raw and normalized microarray data have been deposited in the ArrayExpress database with reference E-MEX-3581.

MICROARRAY DATA PREPARATION AND NORMALIZATION
For each microarray, measured net fluorescence intensities (median values, with background subtracted) were Lowess normalized. Spots with net intensity less than 200 across the entire time-course were filtered away. The filtered data were log 2 transformed and subjected to median subtraction and z-score normalization (Quackenbush, 2002;Cheadle et al., 2003).

STATISTICAL ANALYSIS OF MICROARRAY DATA
To find non-constant (non-zero) genes expressed during murine tooth development (E11.5-P7), statistical analysis of microarray data was carried out, using Spotfire v. 9 Microarray Analysis Software (TIBCO Software Inc, Palo alto, CAL, USA). Microarray data was derived from sets of three to six arrays at each time point. Data from a total of 58 arrays were combined into a single data file and treated as single color data to facilitate statistical analysis of time-courses. False discovery rate (FDR; 0.05) of Benjamini and Hochberg (Benjamini et al., 2001), was used to correct selection of genes for false positives. The ANOVA facility of the Spitfire program was used to select genes which exhibited statistically significant differences in levels of expression (P < 0.05) between the various developmental stages.

BIOINFORMATICS ANALYSIS OF DIFFERENTIALLY EXPRESSED GENES
The 2441 differentially expressed genes with known Entrez Gene ID at all time-points, were uploaded onto IPA (Ingenuity Systems
Network analysis (Thomas and Bonchev, 2010) was used to create graphical representations of molecules interacting at each time-point. The network size was set to 35 nodes/molecules. Network analysis also predicted the upstream and downstream effects of activation or inhibition on other molecules by applying expression values from the dataset (Supplementary data).

GENES WITH EXPRESSION LEVEL CHANGING OVER TIME
Microarray results showed that a total 4362 of non-constant (non-zero) genes are expressed during murine tooth development from E11.5 up to P7; 1921 genes with unknown function and very highly expressed at all time-points with net fluorescence intensities above 10000 (results not shown) and 2441 differentially expressed genes at all time-points with known Entrez Gene ID. These 2441 genes were subjected to further bioinformatics analysis.

TIME COURSES OF EXPRESSION FOR SELECTED GENES AS ASSAYED USING REAL-TIME RT-PCR
Time-course of expression of selected genes (Figure 1) was also monitored using real-time RT-PCR. The results suggest that time-courses assayed by real-time RT-PCR show a similar trend to expression data obtained using microarrays (Figure 2).

Transcription factor analysis
Transcription factor analysis suggested that 19 transcription factors are involved in the transcription of 23 genes during the invagination of the epithelium into the mesenchyme at E12.5 ( Table 2). The transcription factor (TF) Huntingtin (Htt) in the tooth germ regulated other transcription factors e.g., Hif1a, Purb, and C/ebp ( Figure 3A).
During the formation of the enamel knot (E14.5) the number of transcription factors decreased to 3 compared to E12.5-E13.5 ( Figure 3C). At this embryonic stage of tooth development, myc seems to regulate the transcription factors Hif1a, Smarcc1, and Mycn. Mycn regulates Mxl1 (Figure 3C).

Network analysis
Genes expressed during placode formation (E12.5) were associated with network functions such as post-translational modification, cellular growth and proliferation, cell cycle and cell-to-cell signaling ( Table 4 and Supplemental data pages 3-6).

www.frontiersin.org
February 2015 | Volume 6 | Article 47 | 3    During bud formation (E13.5) the genes expressed at this stage were associated with the following network functions: Posttranscriptional modification, protein synthesis and folding, cellular compromise, cell death (Table 4 and Supplemental data pages 7-23).
Genes expressed at bell stages E17.5-E18.5 were significantly associated with networks connected to lipid metabolism, nucleic acid and carbohydrate metabolism, small molecule biochemistry, senescence and energy production ( Table 4 and Supplemental data pages 53-82).
Genes expressed at the post-natal stages (P0-P7) were associated with networks involved in DNA/RNA replication, cell-to-cell signaling, protein synthesis, cell death and free radical scavenging ( Table 4 and Supplemental data pages 84-197).
The canonical pathways associated with the 2441 DE genes are listed in Table 3.

DISCUSSION
Transcription factor and network analysis of the 2441 differentially expressed genes suggests that during the embryonic stages of murine tooth development (E12.5-E13.5) cell proliferation and cell death are highly regulated e.g., expression of transcription regulators like Huntingtin (Htt), Hif1a, Purb, Cbp1, C/ebp, Myc, and Id1. The epithelial and mesenchymal cells during bud-and bell-stages, proliferate, migrate, adhere and communicate through tight junction signaling as shown by the expression and regulation of Ccnd2 on the data set. Ccnd2 also plays a role in Wnt/β-catenin signaling pathway (Liu and Millar, 2010). Wnt/βcatenin is known to play a central role for the morphogenesis of ectodermal appendages such as teeth, hair and exocrine glands (Haara et al., 2012).
During early murine tooth development, ectodermal placodes proliferate and invaginate into the mesenchyme to form an epithelial bud (E12.5-E13.5). The proliferating epithelial cells form a controlled invasion front into the mesenchyme and hold a correct temporal and spatial pattern in the growing bud; this process requires genes that control migration and adhesion between cells e.g., Ctnnb1. Migration of epithelial cells into the underlying mesenchyme requires changes in the cytoskeleton enabling the epithelial cells to migrate and invade the mesenchyme. The cytoskeleton provides both structural scaffolding for cells and a transportation network, where particles move along microtubule and actin highways powered by molecular motors that burn the cellular fuel, ATP (Ikuta et al., 2014).

Frontiers in Genetics | Systems Biology
February 2015 | Volume 6 | Article 47 | 14 2010). Zn deficiency causes growth retardation, reduced bone volume, dental decay. Zinc equilibrium is also required for odontoblast differentiation and dentin formation during dentinogenesis (Lin et al., 2011). During differentiation stages (Wang et al., 2013) (P0-P5) and early mineralization stages (P6-P7) in addition to cell migration apoptosis and cell death, mesenchymal cells facing the epithelial cells at the basal side of the growing bell start to differentiate, elongate following by secretion of dentin around P4-P5, triggering secretion and mineralization initiation by the epithelial cells now differentiating into ameloblasts (P6-P7). Some transcription factors like Myc, Mycn, Htt Stat6, Mxi1, and Drap1 seem to regulate different genes at the post-natal stages (P0, P2, P4, P5, and P6) compared to the embryonic stages (E12.5, E13.5, and E16.5). Esr1 (P4) encode for an estrogen receptor. Estrogen and its receptors are essential for epithelial cell development, proliferation, suggesting that not only genes and their products but also hormones play an important during murine tooth development (Wang et al., 2013). Low estrogen production is associated with increased production of Tgfa, interleukins 1, 6, 8, and 10 leading to periodontal disease (Tezal et al., 2000;Dvorak et al., 2009;Zhang et al., 2011).
It is important to point out that levels of mRNA do not necessarily correlate with levels of translated protein (protein biosynthesis) from the same mRNA (Pascal et al., 2008). The quantitation of levels of mRNA and protein are complementary and also necessary for a complete understanding of how altered gene expression affects cellular physiology (Greenbaum et al., 2003), in addition to miRNAs, long non-coding RNAs (lncRNAs) (Okazaki et al., 2002;Batista and Chang, 2013) that do not code for functional proteins (Rinn and Chang, 2012) may play an important role in spatial positioning of the epithelial/mesenchymal cells during murine tooth development adding layers of complexity to the study of murine tooth development.

AUTHOR CONTRIBUTIONS
ML performed microarray, RT-PCR and bioinformatics analysis. Prof. HO provided exceptional scientific resources, guidance and incredible scientific inspiration. MS helped with qPCR. SN contributed with bioinformatics expertise. EB performed qPCR oligoprimers synthesis and analysis. JR helped with fruitfully discussions editing the manuscript.