Abstract
Objective. Peri-implantitis (PI) is one of the main reasons for dental implant failure. Until now, the etiology and pathogenesis of PI remain unclear. Methods. In this study, we used differentially expressed genes (DEGs) analysis and gene function enrichment analysis to assess the expression profile of peri-implant bone tissue and gingiva in PI public data from the Gene Expression Omnibus (GEO) database. Then, we used gingival tissues from patients with PI and healthy individual to construct gene coexpression networks to reveal the biological functions of the genes in PI using RNA sequencing data. Afterward, key gene modules were selected to reveal the critical biological process or signaling pathway using Hallmark’s gene enrichment and expression analysis of the related pathway members in PI. Results. DEGs were enriched in the formation of cellular responses to external stimuli in bone tissue. Cytokine production, lymphocyte activation, immune response-regulating signaling pathway, and blood vessel development were the top GO biology process or pathways of the DEGs in gingival tissue. Weighted gene coexpression network analysis (WGCNA) of RNA-seq data was used to assess the results of correlation analysis between modules and traits and correlation analysis between modules and functions. kMEpurple, kMEgreen, and kMEred modules were selected as the key gene modules. Signaling pathways and gene expression analysis were performed on selected modules, such as IL2/STAT5 signaling pathway, TNFα signaling pathway via NFκB, and angiogenesis were enriched in kMEpurple module. Hedgehog signaling pathway, Wnt β-catenin signaling pathway, and IL2/STAT5 signaling pathway were enriched in kMEgreen module. Peroxisome, IL2/STAT5 signaling pathway, and epithelial-mesenchymal transformation process were enriched in kMEred module. All the enrichment results of key modules contained IL2/STAT5 signaling pathway. Conclusion. Differential gene and enrichment analysis based on public data showed differences in gene expression patterns and biological process between bone and gingival tissues in PI. This spatial-temporal heterogeneity is reflected in the formation of cellular responses to external stimuli, which was enriched in bone tissue, but cytokine production, lymphocyte activation, immune response regulating signaling pathway, and blood vessel development were enriched in gingival tissue. WGCNA and Hallmark gene sets enrichment analysis of the gingival tissue expression profile and showed that IL2-mediated activation of immune cells could be a critical mechanism in PI. As a new clinical treatment alternative, we suggest that IL2/STAT5 pathway blockers could be helpful in the treatment of PI.
1. Introduction
Dental implant, because of its advantages, such as perfect retention, less damage to adjacent teeth, and less foreign body sensation, has been widely used to reconstruct aesthetic and functional problems that result from teeth loss in clinic [1]. However, peri-implantitis (PI) characterized by infection of soft tissue and bone resorption is considered to be the result of an imbalance between the bacterial challenge, and the host response can affect the long-term success rate of dental implant, which is one of the main reasons for the failure of implant [2].
Bacterial invasion in the tissue surrounding implants can trigger an immune response. This response can remove harmful substances such as bacteria and toxins. However, the cytokines, proteases, and prostaglandins produced during this process can accelerate the destruction of tissue around the implants [3]. With the growing popularity of dental implants, PI has attracted considerable attention, but the etiology and pathogenesis are still unclear [4].
Weighted gene coexpression network analysis (WGCNA) [5] aims to find gene modules for coexpression and to explore the relationship between gene networks and phenotypes of interest, as well as the hub genes in the network. WGCNA could avoid the extensive false-positive and false-negative results of prior biological methods and exclude unreasonable statistical filtering in differential gene analysis. WGCNA has been widely used in cancer research, developmental biology, and systems biology [6, 7]. However, WGCNA is rarely used in the study of oral diseases.
In this study, we analyzed the expression profiles of peri-implant bone tissue and gingiva in PI from the GEO database using differentially expressed gene analysis and gene function enrichment analysis. Then, WGCNA was used to reveal the biological functions of the genes in PI. Owing to the similarity between genes and genes in the expression profiles data of probe-based PCR microarrays, the WGCNA of microarray data was failed to screen out gene modules and gene module members with biological significance. Therefore, we collected gingival tissues from patients with PI and healthy individual to construct gene coexpression network by RNA sequencing data. Next, we selected key modules to reveal the critical biological process or signaling pathway by Hallmark gene sets enrichment analysis and expression analysis of related pathways members in PI. In this study, we wish that high-throughput sequencing technology can be used to analyze the core issues that are plaguing the study of PI. Further, avoiding problems caused by defective technical means within basic medical research of peri-implantitis.
2. Material and Methods
2.1. Microarray Data Acquisition
Gene expression data were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds/). The datasets GSE57631 and GSE33774 were queried in the database using “peri-implantitis” as the search term. GSE57631contains 2 healthy peri-implantation bone tissues (BT_HI) and 6 peri-implantitis bone tissues (BT_PI) of expression profiling by array. GSE33774 contains 8 healthy peri-implantation gingival tissues (GT_HI) and 7 peri-implantitis gingival tissues (GT_PI) of expression profiling by array.
2.2. Microarray Data Analysis
-score was used to normalize the data. Then we used the differential genes analysis of expression matrix GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/). All the genes expression profiles were normalized using the R software package. Principal component analysis (PCA) and clustering analysis of expression patterns were performed using omicshare (http://www.omicshare.com).
2.3. Functional Enrichment Analysis of Microarray Data
Functional enrichment analysis was performed using the Metascape database (http://metascape.org/). KEGG pathway, GO biological processes, reactome gene sets, canonical pathways, and CORUM ontology were selected as sources. Terms with a value <0.01, a minimum count of 3, and an enrichment were collected and grouped into clusters based on their similarities. Kappa scores were used as a metric of similarity when performing hierarchical clustering on the enriched terms, and sub-trees with a similarity of >0.3 were considered a cluster. The most statistically significant term within a cluster was chosen to represent the cluster. All the visualizations were performed using R software.
2.4. Weighted Gene Coexpression Network Analysis of Microarray Data
WGCNA is an algorithm used in gene coexpression network identification in profiles with different traits. In this section, R software package WGCNA was used to construct weighted coexpression networks to find key gene modules of interested within different traits.
2.5. Sample Collection
Gingival tissue samples of PI and healthy individual admitted to the Affiliated Stomatological Hospital of Guangxi Medical University were collected from December 2017 to December 2018. Inclusion criteria of gingival samples were described in the consensus report of the workgroup 4 in the 2017 World Workshop on the Classification of Periodontal and Peri-Implant Diseases and Conditions [2].
To obtain PI samples, the inflamed soft tissues around the implants were removed during an open surgical debridement following currently approved protocols. Gingival tissues from patients without any clinical infection were used a control. Control patients were operated due to wisdom teeth removals or teeth needed to be removed for orthodontics.
All the patients aged 18 years and above who did not smoke or drink alcohol were included, while we excluded patients with systemic diseases and other oral diseases, such as common mucous membrane disease, jaw cyst, and tumors. This study has been reviewed and approved by the Ethics Committee of Guangxi Medical University (2017-No. 165). All the patients agreed to participate in the study and signed informed consent before surgical intervention.
2.6. RNA Sequencing
Gingival tissue gene expression profiling was performed using RNA sequencing. Library construction was performed following the manufacturer’s instructions provided by Illumina (San Diego, CA, USA). Samples were sequenced on an Illumina HiSeq 2500 instrument.
2.7. Weighted Gene Coexpression Network Analysis and Functional Enrichment Analysis of RNA-Seq Data
Same as WGCNA microarray data. The number of genes arranged in the constructed modules and functional enrichment analysis was performed on genes in these modules. The corresponding genes’ information was mapped to Metascape (http://metascape.org/).
Hallmark gene sets were used as gene annotation source. Since the construction of the Molecular Signatures Database (MSigDB), it has been widely used as biological processes and diseases databases in metabolic disease and cancer. However, the increasing heterogeneity within gene sets is harmful to the utility of the database. Concerned with this situation, the hallmark gene sets were created as a part of MSigDB [8]. Each hallmark conveys a specific biological process and displays a coherent expression, which provides refined inputs for gene enrichment analysis. value ≤0.05 after the correction was used as a threshold. The modules of interest were visualized using R software.
2.8. Statistics
The GraphPad Prism (Prism 8 for Windows, GraphPad Software Inc., San Diego, CA, USA) software was used for statistical analysis. Data obtained from the experiments are reported as the (SD). The difference between the two groups was determined using Student’s -test. A value of less than 0.05 was considered statistically significant.
3. Results
3.1. Differential Expression Gene Analysis from Microarray Data
We analyzed the expression profiles of genes in 2 healthy peri-implant bone tissue samples, 6 peri-implantitis bone tissues samples, 8 healthy gingival tissue samples, and 7 peri-implantitis gingival tissues samples. The results showed that 930 upregulated mRNAs and 1189 downregulated genes were identified in BT_HI-vs-BT_PI. Then 1735 significantly upregulated genes and 715 downregulated mRNAs were identified in GT_HI-vs-GT_PI. The list of top 20 DEGs is showed in Table 1.
To test the quality of the two-trait sample groups within expression profiles, the principal component analysis was used, as showed in Figure 1(a). The results of PCA analysis showed that the healthy bone tissue expression profiles and PI bone tissue of patients’ expression profiles could be well distinguished, but there was a significant overlap between PI gingival tissue and healthy gingival tissue.
(a)
(b)
(c)
The expression pattern clustering the two-trait sample groups were shown in Figure 1(b). Clustering analysis of the expression patterns of genes with significant differences can adequately find the common points of expression among different genes and infer the similarity of gene functions according to the similarity of expression patterns. According to the results of the clustering analysis of expression patterns, the expression trends of gene groups with similar expression patterns in each sample can be expressed by curves. The distance calculation algorithm was used, the sample was Spearman, the gene was Pearson, and the clustering method was Hcluster.
3.2. Functional Enrichment Analysis of DEGs from Microarray Data
To test the biological function of the identified genes, information from differentially expressed genes were applied to the enrichment analysis using Metascape. Reactome gene sets, canonical pathways, CORUM, gene ontology (GO), and Kyoto encyclopedia of genes and genomes (KEGG) were used to obtain comprehensive functional annotations from multiple gene repositories for enrichment analysis.
The top 10 GO terms and pathways with the lowest value of each group were shown in Table 2, and the results of the two groups enrichment were shown in Figure 1(c). DEGs were enriched in the formation of cellular responses to external stimuli in bone tissue. Moreover, cytokine production, lymphocyte activation, immune response regulating signaling pathway, and blood vessel development were the top GO biology process or pathways of DEGs in gingival tissue. It is suggested that the differential biological processes involved in the expression of gingival tissue genes and bone tissue genes may be the mechanism of spatiotemporal heterogeneity in peri-implantitis.
3.3. Weighted Gene Coexpression Network Analysis of Microarray Data
The individual gene variance in each sample was calculated according to the normalization of the expression profile. The unsigned network was constructed, selecting genes with a standard deviation higher than 1.2. The expression profiles and traits data consisted of 23 samples, 18357 genes, and 4 traits. Cluster analysis of all the samples was shown in Figure 2(a).
(a)
(b)
(c)
To ensure that the network is unsigned, the soft threshold value was chosen. The expression profiles were transformed into the adjacency matrix and later transformed into the topological matrix. Based on the topological overlap measure (TOM), we used the average-linkage hierarchical clustering method to cluster genes. According to the standard of a hybrid dynamic cut tree, the minimum number of bases for each gene network module was 30. After determining the gene module with a dynamic splicing method, we calculate the eigengenes of each module, then cluster the modules, merge the nearer modules into new modules, and set the . Only three modules were obtained, as shown in Figure 2(b), in which the grey module is unable to aggregate into the gene set of other modules.
The gene significance of the members of module blue and turquoise was shown in Figure 2(c). The scatter plot of module kME value and gene significance value shows that higher the core value, smaller the value, and the module members are more representative of the module characteristics.
3.4. Weighted Gene Coexpression Network Analysis of RNA-Seq Data
To test the quality of RNA sequencing within expression profiles, the principal component analysis was used, as shown in Figure 3(a). The results of PCA analysis showed that healthy gingival tissue expression profiles and PI gingival tissue of patients’ expression profiles could be well distinguished, compared with the microarray data. Cluster analysis of all the samples was shown in Figure 3(b), and the gene dendrogram with traits was shown in Figure 3(c). It could be seen that genes were allocated to 11 modules, which could be used for function and module correlation analysis. According to the eigengenes of each module, the correlation between these modules and each trait was calculated, as shown in Figure 3(d).
(a)
(b)
(c)
(d)
(e)
kME was used to evaluate the value of effective connectivity between hub genes and to identify module members. Selecting the as the members of the modules, also named hub genes, can represent better the expression trend of the entire module. To reveal the functional correlation with gene module members, all the modules hub members were representative members of the module for gene enrichment analysis; the results were shown in Figure 3(e) and Table 3. Considering the results of the correlation analysis between modules and traits, and correlation analysis between modules and functions, the kMEpurple, kMEgreen, and kMEred modules were selected as the key gene modules. The gene significance of the members of the three modules was shown in Figure 4(a).
(a)
(b)
3.5. Signaling Pathway Analysis and Gene Expression Analysis of Key Gene Modules
The Hallmark gene sets were used as a source to reveal the critical biological process or signaling pathway of key gene modules using gene enrichment analysis. As shown in Figure 4(b), IL2/STAT5 signaling pathway, TNFα signaling pathway via NFκB, and angiogenesis were enriched in kMEpurple module. Hedgehog signaling pathway, Wnt β-catenin signaling pathway, and IL2/STAT5 signaling pathway were enriched in kMEgreen module. Peroxisome, IL2/STAT5 signaling pathway, and epithelial-mesenchymal transformation process were enriched in kMEred module. All the enrichment results of key modules contained IL2/STAT5 signaling pathway. It is suggested that IL2-mediated activation of immune cells could be a critical mechanism in PI.
The gene expression of all the pathway members in different modules was plotted in Figure 5. The expression of some genes was inconsistent in bone and gingival tissues. More gene expressions were inconsistent in microarray and RNA sequencing data. These results confirm once again that difference in the expression of gingival tissue genes and bone tissue genes could be the mechanism of spatiotemporal heterogeneity in PI. The accuracy of gene expression profiles detected with microarray is inconsistent with RNA sequencing-based on next-generation sequencing.
4. Discussion
The development of molecular biology and bioinformatics has revolutionized pathology. Diseases are no longer considered to be caused by abnormal expression or single genes structural changes. Dynamic network relationships between genes and multiple negative feedbacks of signal pathways are considered important regulatory models of homeostasis against pathological factors [9]. To analyze the dynamic changes of gene expression profiles, weighted gene coexpression regulation analysis was created to explore the relationship between gene networks and phenotypes. It includes three steps: calculation of correlation coefficient between genes, constructed coexpression network, and determination of gene module with traits and function [5]. In this study, we firstly used WGCNA to reveal the PI mechanism and then identified IL2/STAT5 signaling pathway as a critical element in three key modules.
Previous studies demonstrated that the expression of cytokines stimulated by exogenous factors, degradation of the extracellular matrix [10], and cellular oxidative stress [11] is PI basic biological processes. This study confirmed these conclusions through the analysis of public databases. Moreover, through the comprehensive analysis of expression profiles, we found differences in gene expression patterns between bone and gingival tissues in PI. The differential genes of bone tissue were related to vascular growth and protein translation, processing, and transportation, while the differential genes in gingival tissues are involved in immune stress response. These differences in biological processes and cellular behavior expose the role of spatiotemporal heterogeneity in PI development. Previous studies have not considered the effects of different tissue behavior on disease occurrence. The gingival tissue is the first barrier against exogenous stimulants such as food debris, dental plaque [12], and implant dissolution [13]. It also plays a pioneering role in local tissue inflammation induced by host stress. Moreover, fibroblasts have been identified to be involved in PI pathogenesis by enhancing vascular and matrix degradation [14]. In clinical practice, uncontrolled PI has a poor prognosis leading to bone resorption, implant loosening, or loss. The proliferation of osteoclasts and apoptosis of osteoblasts is the cellular behaviors leading to this outcome. At the molecular level, active protein translation, processing, modification, and transport are the factors needed to complete this process.
The results of functional module gene enrichment analysis showed that IL2/STAT5 signaling pathway, TNFα signaling pathway via NFκB, and angiogenesis were enriched in kMEpurple module. Hedgehog signaling pathway, Wnt β-catenin signaling pathway, and IL2/STAT5 signaling pathway were enriched in kMEgreen module. Peroxisome, IL2/STAT5 signaling pathway, and epithelial-mesenchymal transformation process were enriched in kMEred module. All the enrichment results of key modules contain IL2/STAT5 signaling pathway, has been suggested that IL2-mediated activation of immune cells play a critical mechanism in PI. IL2 receptor-dependent nuclear transcription factor STAT5 plays a key role in activating Treg cells, and Treg cells negatively regulate the body’s immune response in vivo [15]. They usually play an important role in maintaining self-tolerance and avoiding body-injury by the immune response, but they also participate in immune surveillance and chronic infection [16, 17]. As an important pathway of host stress, the activation of Treg cells mediated by IL2/STAT5 signaling pathway is activated in PI gingival tissue. To some extent, this blocked the cascade amplification of inflammation signal and alleviated local tissue necrosis caused by inflammation and infection. However, this response is also an important way to induce host immune tolerance, causing persistent infection, and repeated clinical conditions [18]. This provides us with some clinical implications: IL2/STAT5 pathway blockers such as CMD178 [19] or pimozide [20] could be helpful in the PI treatment and inhibit Treg cell activation at the pathway and molecular level.
In summary, there were differences in gene expression patterns and enriched biological process between bone and gingival tissues in PI, which means that biological processes and cellular behavior reveal the spatiotemporal heterogeneity in PI development. WGCNA and Hallmark gene enrichment analysis of the gingival tissue expression profile showed that IL2-mediated activation of immune cells could be a critical PI mechanism. As a new clinical treatment alternative, it is suggested that IL2/STAT5 pathway blockers could be helpful in PI treatment.
5. Conclusion
Differential gene and enrichment analysis based on public data showed differences in gene expression patterns and biological process between bone and gingival tissues in PI. This spatial-temporal heterogeneity is reflected in the formation of cellular responses to the external stimuli, which was enriched in bone tissue. In contrast, cytokine production, lymphocyte activation, immune response regulating signaling pathway, and blood vessel development were enriched in gingival tissue. WGCNA and Hallmark gene enrichment analysis of the gingival tissue expression profile showed that IL2-mediated activation of immune cells could be a critical PI mechanism. As a new clinical treatment, it is suggested that IL2/STAT5 pathway blockers could be helpful in PI treatment.
Data Availability
The datasets used and analyzed during the current study are available from the corresponding author on reasonable request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Authors’ Contributions
Li Tang and Hailun Zhou contributed equally to this work.
Acknowledgments
This study was funded by National Natural Science Foundation of China (Grant No.81460108).