Comparison of gene expression profiles between dental pulp and periodontal ligament tissues in humans.

There are anatomical and functional differences between human dental pulp (DP) and periodontal ligament (PDL). However, the molecular biological differences and function of these tissues are poorly understood. In the present study, we employed a cDNA microarray array to screen for differentially expressed genes (DEGs) between human DP and PDL tissues, and used the online software WebGestalt to perform the functional analysis of the DEGs. In addition, the STRING database and KEGG pathway analysis were applied for interaction network and pathway analysis of the DEGs. DP and PDL samples were obtained from permanent premolars (n=16) extracted for orthodontic purposes. The results of the microarray assay were confirmed by RT-qPCR. The DEGs were found to be significantly associated with the extracellular matrix and focal adhesion. A total of 10 genes were selected to confirm the results. The mRNA levels of integrin alpha 4 (ITGA4), integrin alpha 8 (ITGA8), neurexin 1 (NRXN1) and contactin 1 (CNTN1) were significantly higher in the DP than in the PDL tissues. However, the levels of collagen type XI alpha 1 (COL11A1), aggrecan (ACAN), collagen type VI alpha 1 (COL6A1), chondroadherin (CHAD), laminin gamma 2 (LAMC2) and laminin alpha 3 (LAMA3) were higher in the PDL than in the DP samples. The gene expression profiles provide novel insight into the characterization of DP and PDL tissues, and contribute to our understanding of the potential molecular mechanisms of dental tissue mineralization and regeneration.

Introduction

Dental pulp (DP) tissue, termed as the 'ectomesenchyme', is derived from ectodermal cells that grow on the periphery of the neural tube, migrate to the oral position and then differentiate into cells of the mesenchymal phenotype (1). Epithelial cells form ameloblasts and odontoblasts, DP and periodontal ligament (PDL) (2). DP is responsible for the maintenance and repair of the periodontal tissue and its related immune system, and it has a high regenerative potency and responds to various types of damage (3). PDL connects the tooth and the alveolar jaw bone in the area surrounding the root surface, and functions as continuous support, attachment, proprioception and physical protection for the teeth, and minimizes tissue damage arising from trauma and infection (4). Owing to the anatomical and functional differences between human DP and PDL, it is reasonable to assume that there are also differences in the gene expression profiles of these tissues.

Previous studies have indicated that the gene expression patterns of mesenchymal stem cells (MSCs) derived from dental tissues are different from those of other tissue by employing genome-wide gene expression profiling and gene ontology analysis (5). Differentially expressed proteins have also been demonstrated between dental and non-dental ovine MSC populations from the same donor, which may attribute to their unique growth and capacity to generate structures resembling the specific microenvironments from which they were derived in vivo (6,7).

Recent findings suggest that human DP-derived stem cells (DPSCs) and PDL-derived stem cells (PDLSCs) have the ability to regenerate a dentin/pulp-like and cementum/PDL-like complex, respectively when they are transplanted into the subcutaneous space of immunocompromised mice (8–11). Stem cell-based dentistry has emerged as a promising alternative for the development of regenerative therapies, which have unavoidable limitations and the effects of which have not yet been fully determined (12).

The cDNA microarray technique can provide global profiles of gene expression and facilitate the evaluation of large-scale genes simultaneously. This method has been used in dental studies to compare differentially expressed genes (DEGs) among various types of stem cells (13), or tissues (14–16), or diseases (17–19). However, the differences in gene expression profiles between DP and PDL tissue have not yet been fully elucidated. The use of tissue samples provides more information of the actual situation as the interactions between different cell types can be important for the function of tissues.

Therefore, the present study aimed to evaluate and compare the gene expression patterns in DP and PDL tissues from human permanent teeth, and to identify their molecular biological differences and functions. Furthermore, the results may provide insight into the potential molecular mechanisms of dental tissue regeneration.

Materials and methods

Gene expression data

A gene expression data set (accession no. GSE50639), which included 3 DP and PDL tissues, was downloaded from the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo). Gene expression levels were measured through the Affymetrix Human Gene 1.0 ST Array beadchip platform (Affymetrix, Santa Clara, CA, USA) (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL6244). Platform annotation files were also acquired.

DP and PDL samples

Tissues were obtained from healthy permanent premolars (n=16; from 8 males and 8 females, aged 11–14 years) extracted for orthodontic purposes under approved guidelines set by Nanjing Children's Hospital, Nanjing, China. Written informed consent was obtained from all the respective parents or legal guardians.

The DP and PDL samples used for the experiment were collected according to a previously described procedure (8,9). Briefly, tooth surfaces were cleaned with sterile water, and PDL tissues were gently separated from the middle-third of root with a scalpel. The root was then cut around the cementum-enamel junction using sterilized dental fissure burs, and fractured off along the cutting line with sharp-edged pliers to reveal the pulp chamber. The DP tissues were carefully removed from the crown and root. The extracted PD and PDL samples were then immediately frozen and stored in liquid nitrogen.

Microarray data analysis

The expression data were generated using Affymetrix Expression Console software version 1.4 (Affymetrix). The Robust Multi-array Average (RMA) algorithm implemented through the Affymetrix Expression Console software was used to normalize the data. A one-way ANOVA was performed on the RMA expression values to determine whether genes were differentially expressed between DP and PDL groups. A multiple-testing correction was applied to the p-values of the F-statistics to adjust the false discovery rate (20). Genes with adjusted F-statistic p-values of <0.05 were extracted. Microarray analysis identified 1,405 genes with differences in expression of ≥2-fold, 920 and 485 of which were more abundant in the DP and PDL tissue, respectively. However, only strongly expressed genes in the DP or PDL tissue, which differed by >4- or 2.5-fold from the signal of the DP and PDL tissues, respectively, were selected for further analysis.. In order to classify the co-expression gene group with a similar expression pattern, hierarchical clustering analysis was conducted using Affymetrix Transcriptome Analysis Console (TAC) software. The WEB-based Gene Set AnaLysis Toolkit was performed for the biological interpretation of DEGs (21,22). WebGestalt is a system that facilitates the analysis of sets of genes that can be visualized and organized by a user-selected method. These genes were classified based on data on gene function in the gene ontology of the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database. Protein-protein interaction (PPI) networks represents a significant step in the elucidation of the underlying molecular mechanisms. In our study, PPI networks were constructed for the protein products using information from the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING, version 9.1; http://string-db.org/) (23). Interactions with a score (i.e., required confidence) >0.4 were retained in the network.

Reverse-transcription-quantitative polymerase chain reaction (RT-qPCR)

RNA was isolated from the DP and PDL tissues using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. To determine the number of cDNA molecules in the reverse-transcribed samples, qPCR analyses were performed using the LightCycler system (Roche, Indianapolis, IN, USA). PCR was conducted using 2 µl LightCycler DNA Master SYBR-Green I (Roche), 12.5 µl of reaction mixture, 2 µl of each 5′ and 3′ primer, 2 µl samples. H2O was then added to a final volume of 25 µl. The samples were denatured at 95°C for 10 sec, with a temperature transition rate of 20°C/sec. Four steps were carried out in amplification and fluorescence determination: denaturation at 95°C for 1 sec, with a temperature transition rate of 20°C/sec; annealing for 5 sec, with a temperature transition rate of 8°C/sec; extension at 72°C for 20 sec, with a temperature transition rate of 4°C/sec; and detection of SYBR-Green fluorescence, which reflects the amount of double-stranded DNA, at 86°C for 3 sec. The amplification cycle number was 35. To discriminate specific from non-specific cDNA products, a melting curve was obtained at the end of each run. Products were denatured at 95°C for 3 sec, and the temperature then decreased to 58°C for 15 sec and increased slowly from 58 to 95°C using a temperature transition rate of 0.1°C/sec. To determine the number of copies of the targeted DNA in the samples, purified PCR fragments of known concentrations were serially diluted and served as external standards that were measured in each experiment. Data were normalized to glyceraldehyde 3-phosphate dehydrogenase (GAPDH) levels in the samples. The primer sequences used for PCR are listed in Table I.

Table I

Specific primer sequences used for RT-qPCR.

Name	Primers sequences (5′→3′)	Product size (bp)
ITGA4	F: CACAACACGCTGTTCGGCTA
	R: CGATCCTGCATCTGTAAATCGC	139
ITGA8	F: AGAATGGAGACCTTATTGTGGGA
	R: GAGCCACTTCCGTCTGCTTT	147
NRXN1	F: TAAGTGGCCTCCTAATGACCG
	R: TCGCACCAATACGGCTTCTTT	91
CNTN1	F: CAGCCCTTTCCCGGTTTACAA
	R: TGCTTCTGACCATCCCGTAGT	170
ACAN	F: GTGCCTATCAGGACAAGGTCT
	R: GATGCCTTTCACCACGACTTC	167
COL11A1	F: ACCCTCGCATTGACCTTCC
	R: TTTGTGCAAAATCCCGTTGTTT	128
COL6A1	F: ACAGTGACGAGGTGGAGATCA
	R: GATAGCGCAGTCGGTGTAGG	122
CHAD	F: CGCGGCCTCAAGCAACTTA
	R: TAGGTCAGCTCGGTCAGGTC	95
LAMC2	F: CAAAGGTTCTCTTAGTGCTCGAT
	R: CACTTGGAGTCTAGCAGTCTCT	153
LAMA3	F: CACCGGGATATTTCGGGAATC
	R: AGCTGTCGCAATCATCACATT	165
GAPDH	F: GGAGCGAGATCCCTCCAAAAT
	R: GGCTGTTGTCATACTTCTCATGG	197

ITGA, integrin alpha NRXN1, neurexin 1; CNTN1, contactin 1; ACAN, aggrecan; COL11A1, collagen type XI alpha1; COL6A1, collagen type VI alpha 1; CHAD, chondroadherin; LAMC2, laminin gamma2; GAPDH, glyceraldehyde 3-phosphate dehydrogenase; F, forward; R, reverse.

Statistical analysis

Data from image analysis are presented as the means ± SEM. Statistical comparisons were made using a two-way ANOVA. A value of p<0.05 was considered to indicate a statistically significant difference.

Results

Gene expression profiles in human DP and PDL tissue

Affymetrix Transcriptome Analysis Console software was applied to analyze the cDNA microarray data set (GSE50639). The results indicated that the expression of a total of 1,405 genes was altered by at least 2-fold in one tissue type relative to the other. In the DP tissue, the expression levels of 920 genes were at least double those in the PDL tissues, while in the latter, the expression levels of 485 genes were at least double those in the DP tissue. Ultimately, only strongly expressed genes (529 genes) in the DP or PDL tissue whose expression differed by >4- or 2.5-fold from the signal of the DP or PDL tissue were evaluated further. In the DP tissue, 255 genes were upregulated by at least 4-fold relative to the PDL tissue, while 274 genes were upregulated by at least 2.5-fold in the PDL tissue relative to the DP tissue. Hierarchical clustering and a heatmap of the DEGs in the DP and PDL tissue are shown in Fig. 1, and the up- and downregulated genes in the DP tissues (compared with the PDL tissues) are shown in Table II.

Figure 1

Heatmap of clustering of 1,405 differentially expressed genes in dental pulp (DP) and periodontal ligament (PDL) tissues. Genes expressed with a signal intensity of >529 were included in this analysis. Clustering was carried out using Transcriptome Analysis Console (TAC) software (clustering type: hierarchical clustering; distance metric: Pearson's correlation).

Table II

Up- and downregulated genes in the DP tissues (compared with the PDL tissues).

Upregulated genes in the DP tissues			Downregulated genes in the DP tissue
Name	Gene symbol	Fold change	Name	Gene symbol	Fold change
Anoctamin 1, calcium activated chloride channel	ANO1	6.87	Cathepsin K	CTSK	−6.62
Scinderin	SCIN	28.26	ADAM metallopeptidase with thrombospondin type 1 motif 14	ADAMTS14	−2.64
Bone morphogenetic protein 7	BMP7	5.71	Denticleless E3 ubiquitin protein ligase homolog (Drosophila)	DTL	−3.16
Sorbin and SH3 domain containing 2	SORBS2	9.62	Sushi-repeat containing protein, X-linked 2	SRPX2	−7.94
Leucine-rich repeat containing G protein-coupled receptor 5	LGR5	43.19	Extended synaptotagmin-like protein 3	ESYT3	−2.82
Adherens junctions associated protein 1	AJAP1	10.07	Fibroblast activation protein alpha	FAP	−16.58
Transferrin	TF	53.65	ATPase, class V, type 10A	ATP10A	−2.95
Solute carrier family 16, member 4	SLC16A4	8.58	Collagen type IV alpha 1	COL4A1	−2.87
Family with sequence similarity 134, member B	FAM134B	5.39	Bone gamma-carboxyglutamate (gla) protein	BGLAP; PMF1-BGLAP; PMF1	−3.84
Homer homolog 2 (Drosophila)	HOMER2	4.24	Lysyl oxidase-like 2	LOXL2	−3.11
KN motif and ankyrin repeat domains 1	KANK1	4.27	Integrin alpha 11	ITGA11	−6.81
Chromosome 10 open reading frame 107	C10orf107	8.01	Secreted phosphoprotein 1	SPP1	−14.77
T-cell lymphoma invasion and metastasis 1	TIAM1	7.23	Methylenetetrahydrofolate dehydrogenase (NADP+ dependent) 2	MTHFD2	−3.41
Phosphate regulating endopeptidase homolog, X-linked	PHEX	68.57	Collagen type XVI alpha1	COL16A1	−5.89
Transmembrane protein 156	TMEM156	21.87	Retinol binding protein 4, plasma	RBP4	−7.29
MOB kinase activator 3B	MOB3B	6.42	G protein-coupled receptor 1	GPR1	−3.2
Ras homolog family member U	RHOU	7.76	Nicotinamide N-methyltransferase	NNMT; LOC101928916	−3.94
Retinoic acid receptor responder (tazarotene induced) 1	RARRES1	10.38	Collagen type XVII alpha 1	COL17A1	−2.91
Peroxidasin homolog (Drosophila)-like	PXDNL	10.05	Methylenetetrahydrofolate dehydrogenase (NADP+ dependent) 2	MTHFD2P7; MTHFD2	−2.88
KIAA1324	KIAA1324	10.36	Chondroadherin	CHAD	−4.19
UDP-Gal:betaGlcNAc beta1,3-galactosyltransferase	B3GALT1	10.82	Semaphorin 7A, GPI membrane anchor	SEMA7A	−4.13
Unc-80 homolog (C. elegans)	UNC80	8.46	Dynamin 1	DNM1	−2.91
Abhydrolase domain containing 12B; microRNA 4454	ABHD12B; MIR4454	33.54	Carbonic anhydrase XII	CA12	−2.85
Clusterin	CLU	6.84	Collagen type IV alpha 2	COL4A2	−2.59
WNT1 inducible signaling pathway protein 1	WISP1	5.12	Collagen type VI alpha 1	COL6A1	−3.29
Cadherin 12, type 2 (N-cadherin 2)	CDH12	56.93	microRNA 21; vacuole membrane protein 1	MIR21; VMP1	−5.57
Epithelial cell adhesion molecule	EpCAM	9.81	Matrix metallopeptidase 19	MMP19	−3.26
Sema domain, immunoglobulin domain (Ig),	SEMA3E	46.94	Scavenger receptor class A, member 5 (putative)	SCARA5	−3.77
EPH receptor A5	EPHA5	12.41	MIR155 host gene (non-protein coding); microRNA 155	MIR155HG; MIR155	−2.97
Charged multivesicular body protein 4C	CHMP4C	8.81	CUB domain containing protein 1	CDCP1	−2.71
Bone morphogenetic protein receptor, type IB	BMPR1B	19.06	Tumor necrosis factor receptor superfamily, member 21	TNFRSF21	−3.05
Tripartite motif containing 36	TRIM36	5.13	Myosin ID	MYO1D	−3.13
Aldehyde dehydrogenase 1 family, member A3	ALDH1A3	6.01	Neuroblastoma 1, DAN family BMP antagonist	NBL1; MINOS1-NBL1	−4.74
Nucleic acid binding protein 1	NABP1	−3.17	Gamma-aminobutyric acid (GABA) A receptor beta 2	GABRB2	−3.03
Uridine-cytidine kinase 2; microRNA 3658	UCK2; MIR3658	7.34	RNA, 7SK small nuclear pseudogene 137	RN7SKP137	−4.57
Dystrobrevin binding protein 1	DTNBP1	6.06	Matrix metallopeptidase 3 (stromelysin 1, progelatinase)	MMP3	−5.72
G protein-coupled receptor 63	GPR63	13.14	Coronin 6	CORO6	−2.73
1-Acylglycerol-3-phosphate O-acyltransferase 4	AGPAT4	4.7	CD109 molecule	CD109	−8.66
Distal-less homeobox 3	DLX3	5.61	Parathyroid hormone 1 receptor	PTH1R	−3.29
Cytoplasmic FMR1 interacting protein 2	CYFIP2	12.55	Cadherin-related 23; cadherin-23-like	CDH23; LOC100653137	−3.1
Protein tyrosine phosphatase, receptor type U	PTPRU	4.67	Cholinergic receptor, nicotinic alpha 9 (neuronal)	CHRNA9	−4.68
Nestin	NES	8.01	Acid phosphatase 5, tartrate resistant	ACP5	−7.5
Unc-80 homolog (C. elegans)	UNC80	6.22	Anillin, actin binding protein	ANLN	−3.76
Solute carrier family 2 (facilitated glucose transporter), member 13	SLC2A13	6.24	Aggrecan	ACAN	−5.66
Calpain 6	CAPN6	4.23	Angiopoietin-like 2	ANGPTL2	−7.12
Dehydrogenase/reductase (SDR family) member 3	DHRS3	5.68	Protein tyrosine phosphatase, non-receptor type 22 (lymphoid)	PTPN22	−4.17
Meiotic nuclear divisions 1 homolog (S. cerevisiae)	MND1	6.24	Marker of proliferation Ki-67	MKI67	−4.31
Coiled-coil domain containing 148	CCDC148	4.97	Zinc finger protein 385A	ZNF385A	−2.56
Carboxypeptidase A6	CPA6	22.33	Matrix metallopeptidase 9	MMP9	−22.27
SH2 domain containing 4A	SH2D4A	6.27	Collagen type VI alpha 3	COL6A3	−8.54
Ameloblastin (enamel matrix protein)	AMBN	56.17	Sushi-repeat containing protein, X-linked	SRPX	−5.77
Gastrin-releasing peptide receptor	GRPR	12.65	Collagen type XIII alpha 1	COL13A1	−2.61
Family with sequence similarity 69, member C	FAM69C	6.29	Aldo-keto reductase family 1, member C2	AKR1C2; LOC101930400	−2.9
Scavenger receptor class A, member 3	SCARA3	7.45	Collagen type VI alpha 2	COL6A2	−3.05
Dentin sialophosphoprotein	DSPP	40.4	Collagen type XII alpha 1	COL12A1	−34.05
Myosin XVI	MYO16	9.15	V-maf avian musculoaponeurotic fibrosarcoma oncogene homolog B	MAFB	−2.81
Activating transcription factor 7 interacting protein 2	ATF7IP2	6.18	Absent in melanoma 1	AIM1	−4.65
Serine/threonine kinase 33	STK33	6.99	Dipeptidyl-peptidase 4	DPP4	−5.8
Galactosidase β1-like 3	GLB1L3	4.97	Xg blood group	XG; XGY2	−8.18
Transforming growth factor-beta 1	TGFB1	2.13	Laminin alpha 3	LAMA3	−3.45
Solute carrier family 12 (sodium/potassium/chloride transporter), member 2	SLC12A2	6.02	Cell division cycle 6	CDC6	−2.54
Transforming growth factor-beta 2	TGFB2	3.04	Laminin gamnma 2	LAMC2	−3.29
G protein-coupled receptor 64	GPR64	11.44	Paraneoplastic Ma antigen 2	PNMA2	−3.24
Lipid phosphate phosphatase-related protein type 5	LPPR5	32.35	Glycoprotein (transmembrane) nmb	GPNMB	−3.6
Transmembrane protein 132B	TMEM132B	6.76	X-prolyl aminopeptidase (aminopeptidase P) 2, membrane-bound	XPNPEP2	−3.15
Sodium channel, voltage gated, type VIII alpha subunit	SCN8A	5.21	Collagen type V alpha 1	COL5A1	−6.45
Gamma-aminobutyric acid (GABA) A receptor beta 1	GABRB1	18.02	Asp (abnormal spindle) homolog, microcephaly associated (Drosophila)	ASPM	−2.9
Long intergenic non-protein coding RNA 943	LINC00943	7.72	Solute carrier family 38, member 5	SLC38A5	−2.71
Contactin associated protein-like 4	CNTNAP4	4.1	ATPase, H+ transporting, lysosomal 38 kDa, V0 subunit d2	ATP6V0D2	−7.72
Wingless-type MMTV integration site family, member 5A	WNT5A	4.44	Leucine rich repeat containing 15	LRRC15	−14.31
Gamma-aminobutyric acid (GABA) A receptor gamma 1	GABRG1	11.14	Fibroblast growth factor 9	FGF9	−3.9
Transmembrane protein 38B	TMEM38B	5.8	Syntaxin binding protein 6 (amisyn)	STXBP6	−3.91
Vitrin	VIT	4.77	Mesenchyme homeobox 1	MEOX1	−3.02
Dynein, cytoplasmic 1, intermediate chain 1	DYNC1I1	6.4	Inhibin beta A	INHBA	−2.59
Sparc/osteonectin, cwcv and kazal-like domains proteoglycan (testican) 3	SPOCK3	66.65	Parathyroid hormone-like hormone	PTHLH	−7.65
EF-hand calcium binding domain 1	EFCAB1	4.09	Reticulocalbin 3, EF-hand calcium binding domain	RCN3	−2.65
Phosphorylase, glycogen, liver	PYGL	5.28	Coagulation factor V (proaccelerin, labile factor)	F5	−3.15
Phospholipase C delta 4	PLCD4	7.52	Topoisomerase (DNA) II alpha 170 kDa	TOP2A	−5.27
Fibulin 2	FBLN2	5.8	Protease, serine 35	PRSS35	−4.54
Ras association (RalGDS/AF-6) domain family member 2	RASSF2	4.3	Bone morphogenetic protein 3	BMP3	−15.22
Hyaluronoglucosaminidase 4	HYAL4	9.39	B-cell CLL/lymphoma 11B (zinc finger protein)	BCL11B	−2.84
Integrin alpha 4 (antigen CD49D, alpha4 subunit of VLA-4 receptor)	ITGA4	6.2	Collagen type XI alpha 1	COL11A1	−23.93
Integrin alpha 8	ITGA8	5.52	Wingless-type MMTV integration site family member 2	WNT2	−3.36
Neurexin 1	NRXN1	5.64	Dickkopf WNT signaling pathway inhibitor 2	DKK2	−5.21
Potassium channel, subfamily K, member 2	KCNK2	18.13	Bone morphogenetic protein 8a	BMP8A	−2.22
Contactin 1	CNTN1	4.73	Transforming growth factor beta 3	TGFB3	−2.78
Dentin matrix acidic phosphoprotein 1	DMP1	5.88	Chemokine (C-X-C motif) ligand 13	CXCL13	−3.51
Integrin alpha 2 (CD49B alpha 2 subunit of VLA-2 receptor)	ITGA2	2.35	Chemokine (C-X-C motif) ligand 2	CXCL2	−2.33
Collagen, type XI, alpha 2	COL11A2	4.26	Chemokine (C-X-C motif) ligand 8	CXCL8	−3.19
Reelin	RELN	8.59	Chemokine (C-X-C motif) receptor 1	CXCR1	−2.76
Solute carrier family 1 (neuronal/epithelial high affinity glutamatetransporter, system Xag), member 1	SLC1A1	4.81	Collagen, type XI, alpha 1	COL11A1	−23.93
Potassium intermediate/small conductance calcium-activated channel, subfamily N, member 4	KCNN4	4.15	Tenascin N	TNN	−21.91
ATPase, Na+/K+ transporting, alpha 2 polypeptide	ATP1A2	5	Integrin-binding sialoprotein	IBSP	−5
MET proto-oncogene, receptor tyrosine kinase	MET	6.93	Carboxypeptidase A3 (mast cell)	CPA3	−4.09
Nitric oxide synthase 1 (neuronal)	NOS1	4.43	Insulin-like growth factor 1 (somatomedin C)	IGF1	−5.2
Poliovirus receptor-related 3	PVRL3	5.7	Interleukin 1β	IL1β	−4.78
Claudin 10	CLDN10	5.11	Phospholipase C, beta 1 (phosphoinositide-specific)	PLCB1	−2.68
Claudin 11	CLDN11	6	Cystatin SN	CST1	−5.14
Neuronal growth regulator 1	NEGR1	6.06	Cystatin S	CST4	−2.67
Catenin (cadherin-associated protein), alpha 2	CTNNA2	9.47	Lysozyme	LYZ	−2.86
KIT ligand	KITLG	5.53	Selectin E	SELE	−2.73
Bone morphogenetic protein 5	BMP5	3.92	Peroxisome proliferator-activated receptor gamma	PPARG	−2.56
Bone morphogenetic protein 6	BMP6	3.31	Fas cell surface death receptor	FAS	−2.75
			Matrix metallopeptidase 2 (gelatinase A, 72 kDa gelatinase, 72 kDa type IV collagenase)	MMP2	−3.43
			Cyclin-dependent kinase 6	CDK6	−2.83
			Interleukin 7	IL7	−2.65
			CD4 molecule	CD4	−2.91
			CD1c molecule	CD1C	−2.52
			Alanyl (membrane) aminopeptidase	ANPEP	−4.4
			Leukemia inhibitory factor	LIF	−3.17

DP, dental pulp; PDL, periodontal ligament.

Gene Ontology (GO) analysis

To analyze the specific biological functions and features of the selected genes, an analysis toolkit (WebGestalt) was applied for GO annotation and enrichment analysis. The DEGs were classified according to biological process (BP), molecular function (MF) or cellular component (CC) using the WebGestalt software package on the basis of hypergeometric tests. The Resulting BP, MF and CC networks are shown as directed acyclic graphs (DAG), which are color-coded (red for p-values <0.05) (Fig. 2A). Biological process enrichment was found for genes associated with cell adhesion, extracellular matrix (ECM) organization, regulation of anatomical structure, morphogenesis and system development. Molecular function enrichment was discovered for genes associated with heparin binding, calcium ion binding, integrin binding and metalloendopeptidase activity. Cellular component enrichment was detected for genes associated with the plasma membrane and ECM. Significantly enriched GO categories under biological process, molecular function and cellular component are indicated in Fig. 2B–D. In the biological process category, the GO terms of biological regulation (257 genes), multicellular organismal process (235 genes), response to stimulus (219 genes), related to metabolic process (211 genes) and developmental process were enriched. In the molecular function category, GO terms related to protein binding (213 genes) and ion binding (167 genes) were enriched. In the cellular component category, GO terms related to membrane (254 genes) were enriched.

Figure 2

The identified genes were analyzed according to Gene Ontology (GO) enrichment using WebGestalt (WEB-based GEne SeT AnaLysis Toolkit). (A) Directed Acyclic Graph (DAG) of the enriched GO categories under biological process. Each node shows the name of the GO category, the number of gene in the category and the p-value indicating the significance of enrichment. The red colour represents p-values <0.05. (B–D) GO classification from the biological process, molecular function and cellular component ontology enrichment. The number of proteins enriched in each GO term is shown on the top of each bar.

Interaction network and pathway analysis of DEGs

To further detect the biological function of the DEGs in the DP and PDL tissue, an interaction network was carried out using STRING 9.1 (23), as shown in Fig. 3A. The interaction network was generated on the basis of experimental and database knowledge. Markov cluster algorithm (MCL) was used to find group associations between these DEGs. The inflation factor was set as 1 on a scale of 1–5.

Figure 3

Interaction network and pathway analysis of differentially expressed genes (DEGs). (A) Clustering of the protein-protein interaction network generated by STRING using the Markov cluster algorithm (MCL) algorithm; (B) The significantly enriched KEGG pathways 'ECM-receptor interaction' identified by WebGestalt.

WebGestalt was used to identify the significantly enriched KEGG pathways. The shared enriched pathways, including ECM-receptor interaction, protein digestion and absorption, focal adhesion and cell adhesion molecules (CAMs) (Fig. 3B and Table III) were determined at the significance levels of p<0.05 in WebGestalt. Among the enriched pathways, the upregulated pathways in the DP group included CAMs and salivary secretion pathways and the downregulated pathways in the DP group included ECM-receptor interaction, protein digestion and absorption and focal adhesion. Compared with GO analysis, KEGG pathway analysis provides biological information in a more detailed and specific manner. Furthermore, good concordance was observed from the STRING interaction network and the enriched functional modules identified by the KEGG and GO enrichment analyses.

Table III

List of enriched KEGG pathways of the differentially expressed proteins.

KEGG pathway	No. of hits	Expressed gene participating in the pathway	Statistics for the enrichment of the pathway
ECM-receptor interaction	18	COL11A2, ITGA8, ITGA4, RELN, COL11A1, COL6A3, SPP1, CHAD, LAMA3, COL6A1, COL5A1, ITGA11, COL4A1, TNN, LAMC2, IBSP, COL4A2, COL6A2	C=82; O=18; E=2.12; R=8.47; rawP=2.55e-12; adjP=2.63e-10
Protein digestion and absorption	16	COL6A2, COL17A1, COL4A2, COL11A2, CPA3, COL11A1, SLC1A1, COL12A1, COL6A3, KCNN4, COL6A1, COL5A1, DPP4, XPNPEP2, COL4A1, ATP1A2	C=76; O=16; E=1.97; R=8.12; rawP=8.33e-11; adjP=4.29e-09
Focal adhesion	20	COL6A2, COL4A2, COL11A2, ITGA8, ITGA4, RELN, COL11A1, COL6A3, SPP1, CHAD, LAMA3, COL6A1, COL5A1, IGF1, ITGA11, COL4A1, TNN, MET, LAMC2, IBSP	C=195; O=20; E=5.05; R=3.96; rawP=1.80e-07; adjP=6.18e-06
Amoebiasis	11	LAMA3, COL4A2, COL11A2, COL5A1, IL8, TGFB3, COL4A1, COL11A1, PLCB1, IL1B, LAMC2	C=104; O=11; E=2.69; R=4.08; rawP=8.03e-05; adjP=0.0021
Salivary secretion	8	NOS1, CST1, PLCB1, SLC12A2, ATP1A2, KCNN4, CST4, LYZ	C=78; O=8; E=2.02; R=3.96; rawP=0.0009; adjP=0.0185
Cell adhesion molecules (CAMs)	10	CD4, ITGA8, ITGA4, PVRL3, SELE, NRXN1, CLDN10, CLDN11, NEGR1, CNTN1	C=126; O=10; E=3.26; R=3.06; rawP=0.0016; adjP=0.0265
Pathways in cancer	18	MMP9, FGF9, COL4A2, PPARG, WNT2, CTNNA2, WNT5A, FAS, LAMA3, IL8, MMP2, IGF1, TGFB3, KITLG, COL4A1, MET, LAMC2, CDK6	C=319; O=18; E=8.27; R=2.18; rawP=0.0018; adjP=0.0265
Hematopoietic cell lineage	7	IL7, KITLG, CD4, IL1B, ITGA4, CD1C, ANPEP	C=80; O=7; E=2.07; R=3.38; rawP=0.0047; adjP=0.0543
Cytokine-cytokine receptor interaction	14	IL7, LIF, CXCL13, BMP7, TNFRSF21, IL8, TGFB3, KITLG, MET, IL1B, CXCR1, INHBA, BMPR1B, FAS	C=248; O=14; E=6.43; R=2.18; rawP=0.0054; adjP=0.0543
Rheumatoid arthritis	7	ATP6V0D2, MMP3, IL1B, CTSK, IL8, ACP5, TGFB3	C=82; O=7; E=2.12; R=3.29; rawP=0.0054; adjP=0.0543

KEGG, Kyoto Encyclopedia of Genes and Genomes; ECM, extracellular matrix.

Results of RT-Qpcr

We selected 10 genes [integrin alpha4 (ITGA4), integrin alpha8 (ITGA8), contactin 1 (CNTN1), neurexin 1 (NRXN1), laminin alpha3 (LAMA3), laminin gamma2 (LAMC2), collagen type XI alpha1 (COL11A1), collagen type VI alpha3 (COL6A3), collagen type VI alpha1 (COL6A1) and chondroadherin (CHAD)] to validate the expression data from microarray analysis using SYBR-Green based RT-qPCR. The expression of these genes was significantly altered in various functional modules, as shown in Table III. The 10 genes with expression levels differing by at least 2-fold between the DP and PDL tissues were selected. The results indicated that the mRNA levels of ITGA4, ITGA8, NRXN1 and CNTN1 were significantly higher in the DP compared with the PDL tissues. However, the levels of COL11A1, ACAN, COL6A1, CHAD, LAMC2 and LAMA3 were higher in the PDL tissue compared with the DP tissue (Fig. 4). The mRNA expression levels demonstrated a consistent trend as the cDNA microarray. Taken together, these results suggested that the two types of tissues have a similar mRNA expression profile, paralleling the results determined by RT-qPCR.

Figure 4

RT-qPCR of the relative difference in the mRNA expression extracted from dental pulp (DP) and periodontal ligament (PDL) tissues. mRNA expression levels of integrin alpha 4 (ITGA4), integrin alpha 8 (ITGA8) neurexin 1 (NRXN1), contactin 1 (CNTN1), aggrecan (ACAN), collagen type XI alpha 1 (COL11A1), collagen type VI alpha 3 (COL6A1), chondroadherin (CHAD), laminin gamma 2 (LAMC2) and laminin alpha 3 (LAMA3) were detected. *P<0.05, **P<0.01 and ***P<0.001.

Discussion

With global gene expression profiling, this study uncovered DEGs from DP and PDL tissues. The cDNA microarray results indicated that the expression levels of 1,405 out of 29,096 (4.82%) genes were altered by at least 2-fold in one tissue type relative to the other. Lee et al (14) reported that only 490 out of 33,297 (1.49%) genes were differentially expressed between dental follicle and PDL tissues. This discrepancy may arise from the relative heterogeneity between DP and PDL tissues.

To further analyze the DEGs, functional enrichment analyses were conducted using the WebGestalt tool. Our results revealed that the DEGs were associated with ECM-receptor interaction, protein digestion and absorption and focal adhesion, which were the three most enriched terms. The most enriched term ECM-receptor interaction is a complex network of different combinations of collagens, proteoglycans, hyaluronic acid, laminin, fibronectin and many other glycoproteins, including proteolytic enzymes involved in the degradation and remodeling of the ECM. Dentin ECM proteins play important roles in the dynamics of dentinogenesis (24,25). Detailed studies of dentin matrix proteins may give some insights into unelucidated mechanisms of dentinogenesis.

Focal adhesion formation is initiated upon the binding of adhesion receptors to ECM ligands. The cell adhesion molecule, EpCAM, CNTN1 and NRXN1 had a higher expression in the DP tissues, while CHAD, LAMC2 and LAMA3 were upregulated in the PDL tissues. A recent study presented a close interaction of EpCAM with other cell-cell contact molecules, such as E-cadherin and claudins (26). CNTN1, a prototypical member of the contactin (CNTN) family, is involved, through cis- and trans-interactions with specific cell adhesion molecules, in neural cell migration, axon guidance and the organization of myelin subdomains (27). NRXN functions as synaptic transmission and maturation of contacts (28,29). CHAD, a leucine rich repeat ECM protein with functions in cell to matrix interactions, is associated with both cartilage and bone homeostasis (30). Laminin, one of the major glycoprotein components, may maintain and regulate odontoblast differentiation and enamel crystallization (31).

As the ECM receptors in focal adhesions, integrins are heterodimeric transmembrane proteins that connect the actin cytoskeleton to the extracellular microenvironment and bidirectionally mediate signals across cell membrane (32,33). ITGA2, ITGA4 and ITGA8 were upregulated in the DP compared with the PDL tissues. It has been have found that ITGA2 is positively expressed in DP cells (34) and ITGA4 was stained in the dental neural crest-derived progenitor cells (dNC-PCs) (35). Adherent cells can identify laminin through ITGA2 for adhesion and differentiation (36), and recognize fibronectin via ITGA4 for attachment and migration (37). However, ITGA11 was upregulated in the DPL compared with the PD tissues. A previous study demonstrated that the ITGA11 knockout mice were characterized by a disorganized PDL, which may be attributable to the disturbed matrix metallopeptidase synthesis, and greatly reduced cell adhesion and spreading on collagen I (38). The regulation of matrix metalloproteinase-13 (MMP-13) and cathepsin K is ITGA11-dependent, which is involved in the coordinated extracellular and intracellular collagen proteolysis (39). In addition, integrin-binding sialoprotein (IBSP), which is a major structural protein of the bone matrix (40), was found to be upregulated in the PDL compared with the DP tissue in this study.

Some of the genes that were relatively strongly expressed in PDL tissues were associated with the degradation of the ECM. MMPs, the major players in collagen breakdown, have been identified in periodontal inflammation (41). The robust expressions of collagenases cathepsin K (CTSK) and MMP-3, -9 and -19, which were higher in the PDL tissue than in the DP tissue, are likely to play a part in the turnover of ECM in normal or pathological processes (41–43). A disintegrin and metalloproteinase with thrombospondin motif (ADAMTS14), which is zinc-dependent metalloproteinase and a member of the ADAMTS family of extracellular proteases, is involved connective tissue remodeling and inflammation (44). PDL cells may play a role in both the production and degradation of versican through the secretion of ADAMTS1, ADAMTS4 and ADAMTS5 (45).

In this study, the Wnt pathway member, Wnt family member 2 (WNT2), and Dickkopf-related protein 2 (DKK2), as well as the transforming growth factor (TGF)β/bone morphogenetic protein (BMP) pathway member, BMP3, BMP8A and TGFβ3 were upregulated in the PDL tissues. Cells in the periodontal complex are Wnt responsive, and removing an important member of the Wnt signaling network gives rise to a pathological widening of the PDL space (46). However, other BMP family members (BMP7, BMP6 and BMP5) were upregulated in the DP tissues.

Of note, the PDL expressed more genes associated with inflammation or immune reaction than did the DP tissues. For example, CXCL2, 8 and 13, which are associated with chemotaxis, were upregulated in the PDL tissues. CXCL13, constitutively expressed in secondary lymphoid tissue, is a potent lymphoid chemokine (47). Studies have shown that CXCL13 is associated with B-cell recruitment in chronic inflammatory periodontal lesions (48). The results of this study suggest that anti-CXCL13 may be a promising approach to modulate pathogenic immune responses in PDL tissues.

Dentin sialophosphoprotein (DSPP), the most abundant non-collagenous protein in dentin, is a marker for DPSC differentiation into odontoblasts and is essential for the normal mineralization of dentin (24,25,49,50). DSPP is processed by proteases into three primary domains: dentin sialoprotein (DSP), dentin phosphoprotein (DPP) and dentin glycoprotein (DGP). The dentin matrix proteins (DMPs) lead to tissue calcification due to inherent calcium binding properties in the ECM (50). Both DSPP and DMP1 were upregulated in the DP compared with the PDL tissue in this study.

Genes associated with osteogenic [osteopontin (OPN and osteocalcin (OCN)], osteoclastic [tartrate-resistant acid phosphatase (TRAP)] and chondrogenic (ACAN) functions were more strongly expressed in the PDL than in the DP tissue. OPN, indicated as Spp-1, is a multifunctional sialic acid-rich phosphorylated glycoprotein. OCN, shown as bone Gla protein, is a major non-collagenous protein. TRAP, known as ACP5, is involved in dissolugion of bone mineral though extracellular acidification (51). ACAN, one of the major components of the ECM, is mainly responsible for the high resistance to compression of the load-bearing tissue (52).

RT-qPCR analyses were performed to verify our cDNA microarray results. ITGA4, ITGA8, NRXN1 and CNTN were upregulated in the DP relative to the PDL tissue, while COL11A1, ACAN, COL6A1, CHAD, LAMC2 and LAMA3 were overexpressed in the PDL tissue compared with the DP tisue. These findings are in line with the microarray results.

In conclusion, this study compared gene expression profiles between DP and PDL tissues from human permanent teeth. Although only the RNA from the entire PD and PDL tissues was detected in this study, and not from the individual cell that constitutes these tissues, we consider that our results provide some novel insight into the characterization of DP and PDL tissues, and provide the potential molecular mechanisms concerning dental tissue mineralization and regeneration. The two types of tissue expressed specific genes related to their functions. The knowledge generated from this study demonstrated the differences between the DP and PDL tissues at the molecular biological level and may narrow down the field of potentially important signaling pathways for clinically relevant tissue regeneration.