Research Article

# Identification of significantly different modules between permanent and deciduous teeth by network and pathway analyses

Published: November 03, 2016
Genet. Mol. Res. 15(4): gmr15047959 DOI: https://doi.org/10.4238/gmr15047959
J. Kang, R. Bai, K. Liu, X.P. Ji, Y. Li, J.Y. Han (2016). Identification of significantly different modules between permanent and deciduous teeth by network and pathway analyses. Genet. Mol. Res. 15(4): gmr15047959. https://doi.org/10.4238/gmr15047959
2,868 views

### Abstract

The differences in genetic backgrounds between deciduous and permanent teeth might contribute to the differences in developmental processes, histological characteristics, and tooth life cycles. Here, we attempted to identify significantly different modules between permanent and deciduous teeth via network and pathway analyses. We identified 291 differentially expressed genes (DEGs) between permanent and deciduous teeth using significance analysis of microarray methods. Co-expression networks of DEGs were constructed by weighted gene co-expression network analysis (WGCNA). Three pathways with a significant number of DEGs and P value <0.01 were identified. Integrated co-expression network and pathway (pathway and adjacent gene) analyses were used to extract three pathway-related modules: the calcium signaling pathway-related, ECM-receptor interaction pathway-related, and neuroactive ligand-receptor interaction pathway-related modules. We also attempted to analyze the topological centralities (degree, closeness, stress, and betweenness) of co-expression networks and modules. Four genes (TMEM229A, PPAPDC1A, LEPREL1, and GAD1) and three pathway-related modules that were significantly different in the deciduous and permanent teeth showed similar properties and good centralities. The relative expression levels of key genes were validated, and the differential expression of TMEM229A, LEPREL1, and GAD1 was confirmed by reverse transcription-polymerase chain reaction (P < 0.05). In conclusion, the results of this study may provide a greater understanding of the molecular pathogenesis and potential biomarkers of the progression from deciduous to permanent teeth.

### INTRODUCTION

Deciduous teeth, also known as baby teeth or primary teeth, comprise the first set of teeth that grow during early development in humans; however, these are usually lost and replaced by permanent teeth at around 6 years of age (Yamada et al., 2011). Permanent teeth or adult teeth are the second set of teeth that exhibit more pain sensitivity compared to deciduous teeth, because of a difference in the number and innervations of their neural components (Ahovuo-Saloranta et al., 2013). A large number of people around the world face the problem of decay in the deciduous and permanent teeth; preventive treatment options for this include brushing of the teeth, topical fluoride applications, and the use of dental sealants (Oulis et al., 2011). The majority of decay in children and adolescents is concentrated on the biting surfaces of back teeth, which has a serious negative effect on the quality-of-life of the person having a decayed tooth (Masterson et al., 2014). Meanwhile, tooth decays in permanent and deciduous teeth are morphologically, histologically, and functionally different, which may be attributed to differences in their gene expressions (Hunter et al., 2000).

Previous studies have revealed the involvement of several genetic signatures in tooth development and disease (Mu et al., 2013; Zhang et al., 2014). The identification of changes in gene expression facilitate an informative description of the biology of tooth disease, allowing for the generation of new hypotheses and identification of genetic signatures that provide insight into understanding, as well as the diagnosis and treatment of disease (Doniger et al., 2003; Dawson and Kendziorski, 2012). Several effective methods, including a co-expression network and protein-protein interaction (PPI) network, have been designed for this purpose. Networks can be used to detect differentially expressed gene pairs in high-throughput data measured under two or more conditions within a single study or across multiple studies; moreover, these data can be used to mine unknown connections in incomplete networks (Mukhopadhyay et al., 2012).

However, a certain number of significant interactions between key genes in significant pathways remain to be explained (Nibbe et al., 2011). This may be resolved by evaluating sub-networks or modules of complex networks (Wu et al., 2014). Functions of individual genes and gene-gene interactions could be detected and studied in detail in smaller modules. Therefore, in this study, we attempted to identify pathway-related modules based on functional enrichment and sub-network analyses.

The objective of this study was to identify significantly different modules between permanent and deciduous teeth by integrating network strategies and functional pathway analysis. This was achieved by identifying differentially expressed genes (DEGs) between permanent and deciduous teeth via significance analysis of microarrays (SAM). Subsequently, co-expression networks and significant pathways were obtained based on these DEGs by using weighted gene co-expression network analysis (WGCNA) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional pathway analysis, respectively. Moreover, significant pathway-related modules were extracted from the co-expression network, and topological centralities (degree, closeness, stress, and betweenness) were analyzed to detect the properties of these modules. This would help explain tooth development and disease in the future.

### Data recruitment

The gene expression data obtained from permanent and deciduous tooth samples (three each), E-GEOD-51449 (Kim et al., 2014), were recruited from the ArrayExpress database and presented on the A-AFFY-141-Affymetrix GeneChip Human Gene 1.0 ST Array [HuGene-1_0-stv1] platform (Affymetrix, Santa Clara, CA, USA).

### DEG detection

The changes in gene expressions (DEGs) between permanent and deciduous tooth samples, while accounting for the large volume of data, were identified by SAM (Li and Tibshirani, 2013), wherein a score was assigned to each gene based on the changes in gene expression relative to the standard deviation of repeated measurements. The relative difference d(i) in gene expression was defined as follows:

$d\left(i\right)=\frac{\overline{{x}_{F}}\left(i\right)\text{-}\overline{{x}_{P}}\left(i\right)}{s\left(i\right)+{s}_{o}}$

(Equation 1)

where $\overline{{x}_{F}}\left(i\right)$ and $\overline{{x}_{P}}\left(i\right)$ are the average levels of expression of the gene (i) under experimental and normal conditions, s(i) denotes the standard deviation of repeated expression measurements, and s0 was chosen to minimize the coefficient of variation.

Significant DEGs were further identified by ranking genes in descending order based on d(i) values, such that d(1), d(2), and d(i) were the largest, second-largest, and ith-largest relative differences, respectively. Meanwhile, dp(i) was the ith largest relative difference for permutation p. The expected relative difference, dE(i), was defined as the average of all permutations:

${d}_{E}\left(i\right)=\frac{\sum {}_{p}{d}_{p}\left(i\right)}{n}$

(Equation 2)

For the vast majority of genes, d(i) ≌ dE(i); however, some genes were represented by points displaced from the d(i) = dE(i) line by a distance greater than the threshold Δ. The number of genes identified as significant by SAM increased with the decrease in Δ. In this study, genes that met the criterion (Δ = 2.109) were selected as DEGs.

### Functional enrichment analysis

The biochemical pathways involved in tooth development were analyzed via a KEGG pathway-enrichment analysis of DEGs using the online Database for Annotation, Visualization, and Integrated Discovery (DAVID) tool (Kanehisa and Goto, 2000). In this study, KEGG pathways with P < 0.01 were selected based on the Expression Analysis Systematic Explorer (EASE) test implemented in DAVID (Grossglauser and Vetterli, 2003). The formula calculated by EASE is as follows:

$\text{P}=\frac{\left(\begin{array}{c}a+b\\ a\end{array}\right)\left(\begin{array}{c}c+d\\ c\end{array}\right)}{\left(\begin{array}{c}n\\ a+c\end{array}\right)}$

(Equation 3)

where n = a’ + b + c + d indicated the number of background genes; and a’, a’ + b, and a’ + c denoted the gene number of one gene set in the gene lists, the number of genes in the gene list including at least one gene set, and the gene number of one gene list in the background genes, respectively. a’ was replaced with a = a’ - 1.

### Co-expression network construction

Incorporated co-expression networks are capable of describing gene-gene relationships in detail, and are highly convenient to identify their functions. In this study, the DEG co-expression network was constructed by WGCNA, which describes the correlation structure between gene expression profiles, image data, genetic marker data, proteomics data, and other high-dimensional data (Zhang and Horvath, 2005). The first step was to define a gene co-expression similarity Sij for each pair of genes xi and xj. A signed co-expression measure between xi and xj was used to preserve the sign of correlation, which was defined by a simple transformation of the correlation:

$Sij=\frac{1+cor\left(xi,xj\right)}{2}$

(Equation 4)

The second step was transforming the similarity matrix into an adjacency matrix with the help of an adjacency function. For a weighted network, adjacency was defined by raising the power of the co-expression similarity Sij to β ≥ 1, referred to as soft thresholding. This allows the adjacency to take on values between zero and one in succession to preserve the continuous nature of the co-expression information. The continuous measure to assess gene connection strength is given as ${a}_{ij}={S}_{ij}^{\beta }$, and the weighted adjacency aij between two genes was proportional to their similarity, as seen in the following formula:

$\mathrm{log}\left({a}_{ij}\right)=\beta ×\mathrm{log}\left({S}_{ij}\right)$

(Equation 5)

### Functional module extraction

The detection of functional modules is a major goal of network analysis. In this study, we probed the functional modules from the co-expression network to detect significant genes and pathway-related modules that may play important roles in the progress from deciduous to permanent teeth. This was achieved by mapping the genes in significant pathways onto a global co-expression network; subsequently, pathway-related sub-networks consisting of pathway genes and their adjacent genes, or modules, were extracted.

### Topological centrality analysis

Determining the importance of a particular node or edge in a network is one of the fundamental problems in network analysis; quantifying the centrality and connectivity helps in the identification of portions that may play interesting roles in the network (Bader and Madduri, 2006). In this study, we characterized the biological importance of genes based on the co-expression network and modules using local scale (degree) (Haythornthwaite, 1996) and global scale [stress centrality (Fekete et al., 2005), betweenness centrality (Barthelemy, 2004), and closeness centrality (Wasserman and Faust, 1994)] indices of topological centrality. The genes at the ≥99% quantile distribution in significantly perturbed networks were defined as hub genes.

### Validation by RT-PCR

RT-PCR was used to validate the hub genes in modules. In this study, eight dental pulp samples were obtained from permanent teeth (N = 4) and deciduous teeth (N = 4). Written informed consent was obtained from the children who agreed to participate in this study, as well as their parents. Total RNA was prepared from dental pulp samples using TRIzol regent (Invitrogen, Carlsbad, CA, USA) as described by Kim et al. (2014) and used in first-strand synthesis along with an oligo(dT)18 primer, 2 μL RNasin (40 U/μL), 8.0 μL 5X reverse transcriptase buffer, 8.0 μL dNTPs, and 2 μL avian myeloblastosis virus (AMV) reverse transcriptase (5 U/μL) according to the manufacturer instructions (Invitrogen). cDNA was used as the template, while β-actin was the internal reference for PCR amplification. The primer sequences of four hub genes (TMEM229A, PPAPDC1A, LEPREL1, and GAD1) are listed in Table 1. Each 4-μL reaction mixture used for PCR amplification contained 10 μL 10X PCR buffer I (Invitrogen), 1 μL Taq DNA polymerase High Fidelity (Invitrogen), 3 μL each of the forward and reverse primers, and 8 μL dNTPs. The reaction conditions were set as follows: pre-denaturation at 95°C for 10 min, 40 cycles of denaturation at 95°C for 30 s, annealing at 62°C for 30 s, extension at 72°C for 30 s, and a final extension at 72°C for 10 min. The PCR products were separated by 1.5% agarose gel electrophoresis and analyzed using the Quantity One gel imaging analyzer software (Bio-Rad, Hercules, CA, USA).

### Sequences of primers and lengths of candidate genes.

Genes Primers (5'-3') Length (bp)
Forward Reverse
TMEM229A CTTCGGAATGCACGGCTTTCTG GCTGCCGTACATAAAGAAGGACC 121
PPAPDC1A TGCCCTTGTACTGCGCCATGAT GGAGGATAGTGCTGTCTGTAGC 133
LEPREL1 ATGTGTGAGGGAACTTGCCACC TTGGCACACTCCAGGGCTTTCA 132
β-actin CTCCATCCTGGCCTCGCTGT GCTGTCACCTTCACCGTTCC 268

### Identification of DEGs

SAM analysis of the gene expression profile E-GEOD-51449 identified 291 DEGs between deciduous and permanent tooth samples with a threshold Δ = 2.109. Of these, 24 were upregulated and the remaining were downregulated, which suggested that downregulated genes had a significant effect on the development process (from deciduous teeth to permanent teeth).

### Pathway-enrichment analysis

A pathway-enrichment analysis was performed based on the DEGs between deciduous and permanent teeth samples; the results within a threshold P (P < 0.01) are listed in Table 2. The calcium signaling (P = 2.65 x 10-3), extracellular matrix (ECM)-receptor interaction (P = 5.70 x 10-3), and neuroactive ligand-receptor interaction (P = 7.92 x 10-3) pathways were significantly enriched, with 9, 6, and 10 enriched DEGs, respectively.

### Pathways with a significant number of differentially expressed genes.

Pathway P value Count Genes
Calcium signaling pathway 2.65 x 10-3 9 EGFR, PLCB4, HTR4, PDE1C, PLCG2, GRPR, SLC25A4, GRIN2A, PLCD4
ECM-receptor interaction 5.70 x 10-3 6 IBSP, CD36, ITGA10, ITGA2, ITGB3, SV2A
Neuroactive ligand-receptor interaction 7.92 x 10-3 10 GPR83, GABRG1, GLRB, FPR3, GABRA4, GABRB1, GRPR, HTR4, GRIN2A, GRIK1

### Co-expression network construction and analysis

Many genes play important roles in performing a single biological function and highly co-expressed genes participate in similar biological processes and pathways. Here, WGCNA analysis generated 477 interactions; subsequently, a co-expression network comprising 282 nodes and 477 edges was built using Cytoscape (Figure 1).

Co-expression network comprising 282 nodes and 477 edges constructed by WGCNA. Nodes represent genes and edges denote the interactions between genes. Rhombus nodes indicate genes enriched in significant pathways. The color of a node was decided by the logFoldChange (logFC) value of the gene. The color of the node became darker with the increase in |logFC| value. The yellow nodes indicating TMEM229A, PPAPDC1A, LEPREL1, and GAD1 are hub genes with a higher degree of distribution. Red edges stand for interactions between pathway-enriched genes.

Topological centrality revealed the correlation of a gene as functionally capable of holding communicating nodes together in a biologically complex network. In addition, genes at the top 1% distribution in the co-expression network were defined as hub genes, and were obtained by calculating four types of centralities (degree, closeness, betweenness, and stress) based on the complex network. This revealed that hub genes (or the top 1% distribution genes) were not consistent among the analyses of different centralities in the same gene (Table 3). PPAPCD1A, TMEM229A, and LEPREL1 were identified by the four methods, while GAD1 was identified by the analysis of three centralities; therefore, PPAPCD1A, TMEM229A, LEPREL1, and GAD1 were considered as hub genes of the co-expression network.

### Top 1% ranked genes in co-expression networks based on centrality analyses.

No. Degree Closeness Stress Betweenness
Terms Value Terms Value Terms Value Terms Value
1 PPAPCD1A 206 PPAPCD1A 0.789 PPAPCD1A 72468 PPAPCD1A 0.714
2 TMEM229A 140 TMEM229A 0.658 TMEM229A 38204 TMEM229A 0.379
4 GAD1 33 LEPREL1 0.498 LEPREL1 9582 LEPREL1 0.052

### Pathway module extraction and analysis

The co-expression network comprised of 20 pathway-enriched genes (GRPR, GLRB, GABRA4, GRIK1, SV2A, PDE1C, EGFR, IBSP, ITGA10, CD36, PLCG2, GABRG1, PLCD4, HTR4, ITGA2, PLCB4, GRIN2A, GPR83, GABRB1, and ITGB3). We identified no interactions among these genes, that is, each gene enriched in the calcium signaling (Figure 2A), ECM-receptor interaction (Figure 2C), and neuroactive ligand-receptor interaction (Figure 2E) pathways acted independently. Extraction of the first neighbors of the pathway genes from the co-expression network revealed that the genes formed a module via nodes in the co-expression network (Figure 2B, D, and F). We identified 12, 10, and 14 genes in the calcium signaling pathway, ECM-receptor interaction, and neuroactive ligand-receptor interaction modules, respectively. Moreover, the neighbor genes of all three pathways were the same (TMEM229A, PPAPDC1A, LEPREL1, and GAD1), and were denoted as key genes. Therefore, the three modules could be merged and constructed into a global module (Figure 3).

Modules extracted from co-expression network. Individual pathway-enriched genes in the (A) calcium signaling, (C) ECM-receptor interaction, and (E) neuroactive ligand-receptor interaction pathways; modules of the (B) calcium signaling, (D) ECM-receptor interaction, and (F) neuroactive ligand-receptor interaction pathways. Nodes represent genes and edges denote the interactions between genes. Rhombus nodes indicate genes enriched in significant pathways. Color of each node was decided based on the logFoldChange (logFC) value of the gene. The color of the node was darker with a higher |logFC| value. The circle nodes, representing TMEM229A, PPAPDC1A, LEPREL1, and GAD1, are the key genes of modules.

Global module for three significant modules. Nodes represent genes and edges denote the interactions between genes. Rhombus nodes indicate genes enriched in significant pathways. Color of each node was decided based on the logFoldChange (logFC) value of the gene. The color of the node was darker with a higher |logFC| value. The circle nodes, representing TMEM229A, PPAPDC1A, LEPREL1, and GAD1, are the key genes of modules.

The functions and biological properties of the three pathway modules were further investigated by a topological centrality analysis (Table 4). This comparison revealed a similarity in the betweenness (same) and degree and stress (similar) centralities between the three modules; the closeness centrality of the neuroactive ligand-receptor interaction module was the highest. TMEM229A, PPAPDC1A, LEPREL1, and GAD1 were key genes present in each individual module, which was consistent with the hub genes of co-expression networks and confirmed the validation of our module analysis method.

### Topological properties of modules.

Modules Topological centrality analysis
Degree Closeness Stress Betweenness
Calcium signaling pathway 453 5.834 125058 1.165
ECM-receptor interaction 452 5.168 125058 1.165
Neuroactive ligand-receptor interaction 456 6.842 125048 1.165

### Validation by RT-PCR

The results of the network-based analysis were validated and the key genes were investigated in permanent tooth samples and the relative expression levels of TMEM229A, PPAPDC1A, LEPREL1, and GAD1 in the calcium signaling, ECM-receptor interaction, and neuroactive ligand-receptor interaction pathway modules were analyzed by RT-PCR (Figure 4). We observed that the relative expression levels of the TMEM229A, PPAPDC1A, LEPREL1, and GAD1 genes in permanent teeth were lower than those in deciduous teeth (downregulated DEGs), which was consistent with our outcomes regarding DEGs. The significance of a gene in permanent teeth compared to the same in deciduous teeth was measured by its P value. PPAPDC1A was not differently expressed (P > 0.05), while TMEM229A and LEPREL1 were differently expressed (0.001 < P < 0.05) and GAD1 was significantly differently expressed (P < 0.001). These results indicated that the calcium signaling, ECM-receptor interaction, and neuroactive ligand-receptor interaction pathway modules were significantly different between permanent and deciduous teeth, because of the differently expressed key genes and other DEGs in these modules.

RT-PCR results and relative expression levels of (A) TMEM229A, (B) PPAPDC1A, (C) LEPREL1, and (D) GAD1. D and C represent the permanent and deciduous tooth groups. Expression of a gene in a permanent tooth compared to that in a deciduous tooth was analyzed based on its P value. P > 0.05 indicated that the gene was not significantly differentially expressed; on the contrary, P < 0.05 indicated that the gene was significantly differentially expressed. *0.001 < P < 0.05. **P < 0.001.

### DISCUSSION

We obtained 20 enriched pathways in this study, and extracted three pathway-related modules: calcium signaling, ECM-receptor interaction, and neuroactive ligand-receptor interaction pathways, from the hub genes of the co-expression network and the key genes of modules TMEM229A, PPAPDC1A, LEPREL1, and GAD1. The expression of key genes was validated and the differences between modules were detected by RT-PCR. The results revealed that TMEM229A, LEPREL1, and GAD1 were differentially expressed between permanent and deciduous teeth; therefore, the three modules were significantly different.

TMEM229A, PPAPDC1A, LEPREL1, and GAD1 actively connected common genes together and promoted the formation of small networks using individual genes. Moreover, the roles of TMEM229A, LEPREL1, and GAD1 have been validated in our research. Transmembrane protein 229A (TMEM229A) plays a major role in the binding of zona pellucida and participates in Ca2+ signaling-associated acrosomal exocytosis (Lin et al., 2013a). Glutamate decarboxylase 1 (GAD1) encodes a glutamic acid decarboxylase that catalyzes the production of gamma-aminobutyric acid from l-glutamic acid (Lin et al., 2013b). Moreover, a study has demonstrated that dental pulp from exfoliated deciduous teeth expresses a variety of neural cell markers including β III-tubulin and GAD, and that dental pulp stem cells lose their fibroblastic morphology and develop multi-cytoplasmic processes correlated with β III-tubulin/GAD (Casagrande et al., 2011).

Leprecan-like 1 (LEPREL1), which is mainly concentrated in the endoplasmic apparatus and Golgi complex of the cell and is abundant in the basement membranes, encodes a member of the prolyl 3-hydroxylase subfamily of 2-oxoglutarate-dependent dioxygenases, which plays a critical role in collagen chain assembly (Järnum et al., 2004). Moreover, downregulation of this gene may affect the progression of permanent teeth (Mordechai et al., 2011), which was in accordance with our results. These three genes may not be directly correlated to dental development; however, their functions could affect and regulate the activities of permanent teeth. This is the first study revealing the relationship between TMEM229A, PPAPDC1A, and LEPREL1 expression and tooth progression.

Three modules that were significantly different between deciduous and permanent teeth were identified in this study: the calcium signaling, ECM-receptor interaction, and neuroactive ligand-receptor interaction pathway modules. Calcium binds to calmodulin and stimulates the activity of a number of enzymes, including calcium-calmodulin kinases and calcium-sensitive adenylate cyclases, in the calcium signaling pathway module (Xu et al., 2011). These enzymes transduce the calcium signal and induce short-term biological responses, such as the modification of synaptic proteins and long-lasting neuronal responses that require changes in gene expression (Wang et al., 2015). Moreover, calcium signals are crucial for diverse cellular functions including adhesion, differentiation, proliferation, effector functions, and gene expression (Bonnefond et al., 2015). Fuks (2000) reported that the application of calcium hydroxide to permanent teeth as part of the medication for direct pulp capping induces the deposition of hard tissue as reparative dentin in permanent, but not deciduous, pulp. Recent studies of calcium signal-transduction mechanisms have revealed that, depending on the route of entry into a neuron, calcium differentially affects processes that are central to the development and plasticity of the nervous system, such as the growth of teeth (Ghosh and Greenberg, 1995). Therefore, the calcium signaling pathway module is closely related to permanent teeth.

Teeth develop as epithelial appendages, and their morphogenesis is regulated by epithelial-mesenchymal interactions and conserved signaling pathways that are common to many developmental processes (Mustonen et al., 2002), such as the ECM-receptor and neuroactive ligand-receptor interaction pathways. For example, the ECM consists of a complex mixture of structural and functional macromolecules and plays an important role in tissue and organ morphogenesis and in the maintenance of cell and tissue structure and function (Bonnans et al., 2014). Specific interactions between cells and the ECM results in a direct or indirect control of cellular activities such as adhesion, migration, differentiation, proliferation, and apoptosis and provides an opportunity for ECM to function with others (Jones and Jones, 2000). The ECM-receptor interaction pathway indicates the relationship between structural proteins and exerts a significant effect on some processes, including the progression from deciduous to permanent teeth (Pfefferle and Wray, 2013). Therefore, this module might play an important role in tooth development.

In conclusion, we successfully identified three hub genes (TMEM229A, LEPREL1, and GAD1) and three significantly different modules (calcium signaling pathway module, ECM-receptor interaction module, and neuroactive ligand-receptor interaction module) between deciduous and permanent teeth, which might provide an insight into the potential mechanisms involved in the differentiation, growth, and histological characteristics during tooth development. However, there are some limitations in this analysis. In this study, microarray data were measured by an external agency, while the sample size was quite small. Moreover, these results must be confirmed by experimental and clinical validation.