Autophagy-Related Genes Predict the Progression of Periodontitis Through the ceRNA Network

Introduction

Periodontitis, a chronic infectious disease initiated by dental plaque biofilm and characterized by gingival bleeding and destruction of periodontal supporting tissue, is the main cause of tooth loss in adults.¹ It is widely recognized that periodontitis has two-way relationships with a variety of systemic disorders, including cardiovascular diseases and diabetes.² Microbial infection and host immune-inflammatory responses are the key factors affecting the pathogenesis of periodontitis. However, there are still many problems about the molecular mechanism of periodontitis pathogenic processes, which need further researches.

Autophagy is an evolutionarily highly conserved physiological metabolic process, which can keep the intracellular environment in a state of equilibrium.³ The autophagy process necessitates the regulation of a group of autophagy-related genes (ARG) as well as the recognition of corresponding proteins. These proteins form bilayer membrane-structured autophagosomes, which then fuse with lysosomes to remove intracellular invading pathogens and degrade proteins or damaged organelles.⁴ Autophagy dysfunction contributes to the pathological process of many diseases such as neurodegenerative disorders and tumors.^5,6 Furthermore, the role of autophagy in immune defense and inflammatory regulation has received a lot of attention.⁷ Periodontitis is also an inflammatory disease. Autophagy may influence the onset and progression of periodontitis by regulating periodontal pathogen infection, host immune-inflammatory responses, and alveolar bone metabolism.^8,9 However, the details of the interaction between periodontitis and autophagy remain unclear.

In recent decades, sequencing techniques and bioinformatics have been highly developed and widely used to identify biomarkers of diseases. Recent investigations have revealed critical functions of mRNAs, miRNAs, and lncRNAs in the pathogenesis of periodontitis.^10,11 The competing endogenous RNA (ceRNA) hypothesis, suggested by Salmena et al, revealed a new mechanism of interaction between RNAs.¹² For instance, miRNA can lead to gene silencing by binding to the microRNA response elements (MREs) of mRNA, while lncRNA, pseudogenes, and circular RNA can regulate gene expression by competitively sponging miRNA. The ceRNA network has been linked to the onset and progression of a range of diseases, including cancer, cardiovascular disease, and neurodegenerative disease, according to numerous research.¹³ Therefore, the ceRNA network could be a novel way to investigate the pathophysiology of periodontitis.

This research identified differentially expressed ARGs and used bioinformatics tools to construct a ceRNA network. Furthermore, mRNAs in the ceRNA network were validated using quantitative real-time PCR (qRT-PCR). The findings may provide fresh light on the role of autophagy in periodontitis and lay a theoretical foundation for subsequent studies on periodontitis.

Materials and Methods

Data Sources

The sequence data of periodontitis were retrieved from the Gene Expression Omnibus (GEO) in NCBI (http://www.ncbi.nlm.nih.gov/geo/). Three datasets, GSE10334,¹⁴ GSE16134,^15,16 and GSE54710¹⁷ were obtained using “Homo sapiens [porgn:_txid9606] and periodontitis and gingival tissue” and limiting the sample count to more than 100. The characteristics of the selected GEO database data were shown in Table 1. The HADb (Human Autophagy Database, http://www.autophagy.lu/) was used to download 222 autophagy-related genes (ARGs) (Supplementary Table 1).

Table 1 Characteristics of the Selected GEO Database Data

Differential Expression Analysis

GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/) was used to find differentially expressed genes (DEGs) in GSE10334 and GSE16134. The DEGs were given a threshold of adjusted P-value < 0.05 and |log₂ FC (fold change)| > 0.58. Differentially expressed miRNAs (DE-miRNAs) in GSE54710 were defined by the same criterion. To obtain the differentially expressed ARGs, a Venn diagram of the DEGs and ARGs was constructed using the web application “Calculate and draw custom Venn diagrams”.

Functional Analysis of Differentially Expressed ARGs

To further clarify the potential functional annotation with the differentially expressed ARGs, the Database for Annotation, Visualization, and Integrated Discovery v6.8 (DAVID) (http://david.ncifcrf.gov) was used to perform Gene Ontology (GO) analyses, including biological process (BP), cellular component (CC), and molecular function (MF).¹⁸ Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment was performed based on the KOBased Annotation System (KOBAS) web tool (http://kobas.cbi.pku.edu.cn/index.php).¹⁹ The ClueGO plug-in in Cytoscape v3.8.2 was used to perform interrelation analysis between pathways.²⁰ P<0.05 was regarded as statistically significant.

Construction of the Protein–Protein Interaction (PPI) Network

The PPI network was constructed using the Search Tool for the Retrieval of Interacting Genes v11.0 (STRING) database (http://string-db.org/),²¹ with a confidence score of ≥0.4 as the cutoff. Protein nodes that did not interact with other proteins were deleted from the network. The PPI networks were drawn using Cytoscape v3.8.2.

Construction of the ceRNA Network

Differentially expressed ARGs were performed with miRDB (http://mirdb.org/)²² and StarBase v2.0 (http://starbase.sysu.edu.cn/)²³ to reverse-predict miRNAs. Two databases and DE-miRNAs were used to do intersection. The miRNA-lncRNA interaction was predicted in the StarBase v2.0 and LncBase Predicted v.2²⁴ databases using the miRNAs mentioned above. Only those miRNA-lncRNA relationship pairs found in both databases were chosen to ensure the accuracy of the results. Considering that lncRNAs can serve as nodes of the ceRNA network only in the cytoplasm, we used the lncATLAS database (https://lncatlas.crg.eu/)²⁵ to investigate the intracellular localization of the lncRNAs, excluding those that were only found in the nucleus. The obtained mRNAs, miRNAs, lncRNAs were processed using Cytoscape to construct the ceRNA network.

Validation by qRT-PCR

5 periodontitis patients (average age, 29.8 years old; range 26–33 years old) were recruited from The Second Affiliated Hospital, College of Medicine, Zhejiang University, and were excluded if they had cardiovascular diseases, diabetes, cancer, HIV, or HBV infection. The details of the participants were shown in Table 2. A total of 10 samples were obtained from healthy and diseased gingival tissues in patients. Diseased tissue was collected from periodontal pathology sites (bleeding on probing, probing pocket depth ≥5 mm, and clinical attachment loss ≥3 mm), while healthy tissue was obtained from a neighboring site (no bleeding on probing, probing pocket depth ≤ 3 mm, and no clinical attachment loss). This study was conducted according to the ethical guidelines of the Declaration of Helsinki and was approved by the Research Medical Ethics Committee of The Second Affiliated Hospital, College of Medicine, Zhejiang University. All the participants were made aware of the study and gave their written informed consent.

Table 2 The Details of the Participants

Gingival tissue samples were collected from each participant during surgery. Subsequently, total RNA was extracted using RNAiso Reagent 9180 (TaKaRa, Kyoto, Japan). 4 mRNAs in the ceRNA network were validated by qRT-PCR using TB Green Premix Ex Taq™ RR420A (TaKaRa, Kyoto, Japan) and LightCycler 480 System (Roche, Basel, Switzerland). The primer sequences were designed using Primer 5.0 (Sunya, Hangzhou, China) and reported in Table 3. β-actin was applied as an internal control, and the relative gene expression was determined using the 2−ΔΔCt method.

Table 3 The Primer Sequences for PCR

Statistical Analysis

qRT-PCR results were presented as the mean ± standard deviation. Kolmogorov–Smirnov test was employed to analyze the data’s normality. For data that conformed to the normal distribution, the Student’s t-test was used, and for data that did not adhere to the normal distribution, the Mann–Whitney U-test was used. P<0.05 was considered statistically significant when using GraphPad Prism 9.1.1 (GraphPad Software, San Diego, California USA, www.graphpad.com).

Results

Identification of DEGs, DE-miRNAs, and Differentially Expressed ARGs in Periodontitis

A schematic of the flowchart of this study is shown in Figure 1. The online analysis tool GEO2R was used to evaluate three datasets acquired from the GEO database. Adjusted P-value < 0.05 and |log₂ FC (fold change)| > 0.58 were considered as criteria to screen the DEGs and DE-miRNAs. The volcano plots of DEGs in GSE10334 and GSE16134, as well as DE-miRNAs in GSE54710, were made in Figure 2A–C. 222 ARGs were downloaded from the HADb. The intersections of the two groups of DEGs and ARGs obtained 20 differentially expressed ARGs (Figure 2D), with 10 genes upregulated and 10 genes downregulated (Table 4).

Table 4 Differentially Expressed ARGs of Periodontitis

Figure 1 Study flowchart.

Figure 2 Screening of differentially expressed mRNAs, miRNAs, and ARGs. (A) Volcano plot of differentially expressed mRNAs in GSE10334. Upregulated mRNAs are represented by the red dot, whereas downregulated mRNAs are represented by the green dot. (B) Volcano plot of differentially expressed mRNAs in GSE16134. (C) Volcano plot of differentially expressed miRNAs in GSE54710. (D) A Venn diagram between GSE10334, GSE16134, and ARGs. The coincident section indicates the 20 differentially expressed genes that are shared across the three series.

Enrichment Analysis and PPI Network Construction of the Differentially Expressed ARGs

The DAVID v6.8 online tool was used to conduct GO analyses of differentially expressed ARGs. The enriched GO terms with P < 0.05 were presented in Figure 3A, including apoptotic process, inflammatory response, protein autophosphorylation, negative regulation of the apoptotic process, positive regulation of transcription from RNA polymerase II promoter, neuron apoptotic process, positive regulation of tyrosine phosphorylation of Stat3 protein and ERBB2 signaling pathway in the BP category. For CC, differentially expressed ARGs were mainly clustered in the cytoplasm, endoplasmic reticulum, extracellular space, cytosol, endoplasmic reticulum lumen, extracellular exosome, and lysosome. About MF, differentially expressed ARGs are primarily enriched in processes including protein binding, identical protein binding, protein kinase activity, and transcription regulatory region DNA binding. As demonstrated in Figure 3B, the top 10 KEGG were screened out according to the P-value < 0.05. The interaction network of BP was visualized using the Cytoscape ClueGo plug-in (Figure 4). In addition, a PPI network was built based on the STRING online database and Cytoscape software in Figure 5.

Figure 3 GO and KEGG enrichment analysis of differentially expressed ARGs. (A) Significantly enriched GO terms including BP, CC, and MF. (B) The top 10 KEGG pathways that have been significantly enriched.

Figure 4 Interaction networks of enriched biological processes of differentially expressed ARGs.

Figure 5 PPI network differentially expressed ARGs. The red color symbolizes upregulated genes, while the green color denotes downregulated genes.

Construction of the ceRNA Network

Reverse Prediction of miRNAs Targeted by mRNAs

The miRDB and StarBase v2.0 databases were performed to reverse-predict miRNA of 20 differentially expressed ARGs. The intersecting genes between the predicted miRNAs and 46 DE-miRNAs were then discovered. Finally, mRNA-miRNA interaction pairs involving 4 mRNAs (LAMP2, NFE2L2, NCKAP1, and EGFR) and 3 miRNAs (hsa-miR-140-3p, hsa-miR-142-5p, and hsa-miR-671-5p) were identified. In addition, 4 selected mRNAs were downregulated and 3 miRNAs were upregulated, which agreed with the ceRNA hypothesis.

Prediction of lncRNAs Targeted by miRNAs

We reverse predicted their upstream lncRNAs using 3 miRNAs. A total of 41 lncRNAs was identified using the intersection of the StarBase v2.0 and LncBase Predicted v.2 databases (Table 5). According to the endogenous competition role of lncRNAs is primarily manifested in the cytoplasm, 11 lncRNAs which were located only in the nucleus were excluded. Cytoscape was used to visualize the distribution information for 30 lncRNAs (Figure 6). Finally, a ceRNA of 4 mRNAs, 3miRNAs, 30lncRNAs were identified (Figure 7).

Table 5 The miRNA-lncRNA Network

Figure 6 The subcellular localization of DElncRNAs. Red circles represent the intracellular distribution (nuclear; cytoplasm; both; no data) of 41 lncRNAs (Orange circles).

Figure 7 The ceRNA network. Blue triangles represent 4 mRNAs, yellow circles represent 3 miRNAs and green diamonds represent 30 lncRNAs.

Validation by qRT-PCR

LAMP2, NFE2L2, NCKAP1, and EGFR were validated by qRT-PCR. The result revealed that the expression level of EGFR in diseased gingival tissue of periodontitis patients was significantly lower than that of healthy controls (Figure 8). There was no significant difference in the remaining three mRNAs between the two groups.

Figure 8 The relative expression level of mRNAs. (A) LAMP2. (B) NFE2L2. (C) NCKAP1. (D) EGFR. ***: P < 0.001.

Discussion

Periodontitis, the most common oral infectious disease, not only affects periodontal health and occlusal function by destructing periodontal supporting tissue but also interferes with general health through periodontal pathogens and inflammatory mediators.¹ A great number of studies have shown that genetic factors are strongly linked to periodontitis, and may even affect the speed and severity of disease progression.²⁶ Although simple genetic factors do not cause periodontitis, they can increase the susceptibility of hosts to periodontitis. Therefore, it is meaningful to explore the pathogenesis of periodontitis from a genetic perspective.

Sequencing techniques and bioinformatics have been highly developed and frequently employed to identify biomarkers of diseases in recent decades. The absolute value of log₂FC greater than 1 is commonly used as the screening criterion for differential genes in the literature, whereas periodontitis-related bioinformatics analysis articles frequently choose a lesser value. For instance, Song et al²⁷ chose 0.58 while Jin et al¹¹ selected 0.45 as the screening index. Considering that periodontitis is not only a polygenic disease, but also affected by many factors such as microorganisms, smoking, stress, and systemic diseases.²⁸ If the screening criterion is too stringent, generous potential key genes may be missed, so we also set 0.58 as the screening criterion.

In published bioinformatics articles related to periodontitis, some researchers have analyzed GSE16134 and identified five genes (CTSS, PLEK, IRF-8, PTGS2, and FOSB) that may be implicated in the formation and progression of periodontitis.²⁷ Suzuki et al²⁹ have analyzed three datasets (GSE10334, GSE16134, and GSE23586) and screened several key genes which may be related to chronic periodontitis pathogenesis such as CSF3, CXCL12, IL1B, TAGLN, CD19, IL8, and CD79A. What’s more, upstream genes of biomarker candidates like TNF and FGF2 were predicted. However, there are few bioinformatics analyses on the relationship between periodontitis and autophagy.

It is widely recognized that autophagy has been a research hotspot in recent years. Bullon et al³⁰ found the level of mitochondrial reactive oxygen species and the expression of autophagy gene LC3 increased consistently in peripheral blood mononuclear cells of patients with periodontitis. An et al³¹ also found that the expression of ARGs such as LC3, Beclin1, Atg7, and Atg12 in the periodontal ligament of patients with periodontitis was significantly increased. Indeed, autophagy may influence the onset and development of periodontitis in a variety of ways.⁹ First of all, autophagy can enhance the body’s ability to resist pathogen infection. Secondly, autophagy can regulate inflammation. It is worth noting that moderate inflammation is beneficial to the body, while excessive inflammation can cause damage to the host.³² Last but not least, autophagy can regulate the process of alveolar bone resorption by regulating the differentiation of osteoclasts and osteoblasts.^33,34 Nevertheless, the specific mechanism of autophagy on periodontitis remains to be further studied.

In this work, we found 20 differentially expressed ARGs based on two mRNA microarray datasets and the HADb database, including 10 upregulated (ARSA, CXCR4, DAPK1, DNAJB9, ERN1, FOS, IL24, NAMPT, PRKCQ, and SERPINA1) and 10 downregulated (BAG1, BNIP3, EGFR, HSPB8, LAMP2, NCKAP1, NFE2L2, PTK6, RAB11A, and TP63) genes. To investigate interactions among the differentially expressed ARGs, we used GO analysis, KEGG pathway analysis, and PPI network. GO term analysis revealed these ARGs to be mostly enriched in apoptotic process and inflammatory response (BP), cytoplasm and endoplasmic reticulum (CC), protein binding and identical protein binding (MF). Changes in KEGG were mainly enriched in autophagy – animal and pathways in cancer. We also identified three miRNAs, including hsa-miR-140-3p, has-miR-142-5p, and hsa-miR-671-5p, which were upregulated in the gingival tissue of periodontitis compared to controls. Finally, a ceRNA network consisting of 4 mRNAs (LAMP2, NFE2L2, NCKAP1, and EGFR), 3 miRNAs, and 30 lncRNAs was created after removing those lncRNAs located only in the nucleus.

LAMP2 is critical for the process of fusion between phagosome and lysosome.³⁵ LAMP2-deficient mice showed increased susceptibility to periodontitis due to a reduced maturation process of phagosomes in polymorphonuclear leukocytes.³⁶ NFE2L2 plays a vital role in resisting oxidative stress and maintaining oxidative stress homeostasis. It can not only transcriptionally regulate macroautophagy but also modulate chaperone-mediated autophagy through the regulation of LAMP2A.³⁷ Zhang et al³⁸ showed that NCKAP1 is negatively related to monocyte in periodontitis. EGFR is a transmembrane protein that affects bone homeostasis, proliferation, and differentiation of osteoblasts.³⁹ The inhibitory effects of periodontitis and alveolar bone loss were reversed when silenced EGFR with EGFR siRNA.⁴⁰ In this study, we identified LAMP2, NFE2L2, NCKAP1, and EGFR as key genes between autophagy and periodontitis. These four ARGs were downregulated in periodontitis, which may not be conducive to osteoblast differentiation and increase the host’s susceptibility to periodontitis, thus aggravating periodontitis. Further qRT-PCR results confirmed that the expression of EGFR in the diseased gingival tissue of periodontitis patients was indeed downregulated compared to the healthy controls, which indicated that EGFR may play a significant role in the progression of periodontitis. As a matter of fact, the expression of the other three genes (including LAMP2, NFE2L2, and NCKAP1) in the diseased gingival tissue was also downregulated, but there was no significant difference. The potential reasons need to be taken into account when interpreting these results. On the one hand, the sample size is too small. On the other hand, the diseased tissue and healthy tissue come from the same individual. If tissues are taken from patients with periodontitis and patients without periodontitis, there may be more obvious differences.

Recently, miR-140-3p was manifested as a biomarker for osteoarthritis and osteoporosis.^41,42 Another study found that miR-140-3p mimics significantly reduced protein expression of CXCR4, p-JAK2, and p-STAT3 in IL-1β-stimulated chondrocytes.⁴³ They also suggested that overexpression of miR-140-3p could inactivate JAK2/STAT3 signaling via suppressing CXCR4 in osteoarthritis. Moreover, Mao et al⁴⁴ revealed that overexpression of miR-140-3p inhibited cell survival and differentiation and hastened cell apoptosis in pre-osteoblast. In our research, miR-140-3p was upregulated in periodontitis, which may affect the progress of periodontitis through the above mechanisms. As for mir-142-5p, it has been found that it can promote the osteoclastogenesis of bone marrow-derived macrophages in vitro.⁴⁵ Besides, Sun et al⁴⁶ raised XIST/miR-142-5p/SGTB axis as a novel molecular mechanism in osteoarthritis. We speculated that miR-142-5p may affect bone resorption by promoting osteoclast differentiation, thereby aggravating periodontitis. In addition, a 17-miRNA-based diagnostic signature was identified for periodontitis, which included miR-671-5p.¹¹ A previous study has shown that miR-671-5p played a significant role in inhibiting the characteristics of osteosarcoma cells.⁴⁷ In brief, the above three miRNAs may be linked to the onset of periodontitis.

According to our prediction, lncRNA XIST may be associated with hsa-miR-140-3p, has-miR-142-5p, and hsa-miR-671-5p. LncRNA Xist is the master regulator of X-chromosome inactivation, an epigenetic process that may be important for female development and disease susceptibility.⁴⁸ Furthermore, lncRNA XIST may participate in the development of various cancers and non-cancer diseases by functioning as a ceRNA.⁴⁹ Periodontitis has a gender disparity in incidence rate, with males being at higher risk.⁵⁰ Gender and its related factors, such as genetics, smoking and hormone levels, may affect the incidence rate of periodontitis. It is noteworthy that Chen et al⁵¹ reported that lncRNA XIST inhibited osteogenic differentiation of bone marrow mesenchymal stem cells. Feng⁵² et al also demonstrated that lncRNA XIST participated in osteogenic differentiation of periodontal ligament stem cells. Therefore, lncRNA XIST may be the most influential lncRNA in the ceRNA network we have established.

To be honest, certain limitations still existed in this study. First of all, due to the lack of datasets on the expression of lncRNA in periodontitis, we cannot determine whether the predicted lncRNA was upregulated or downregulated in periodontitis, therefore the established ceRNA network needed further filtering and validation. Secondly, the sample size of our experiment was small, and the expression levels of miRNAs and lncRNAs were not verified. Lastly, non-tumor diseases are often unable to conduct clinical correlation studies and prognosis analysis on account of lacking clinical and prognosis data. Despite these limitations, our study provided more directions to understand the role of autophagy in periodontitis and lay a theoretical foundation for subsequent studies.

Conclusions

In conclusion, we constructed a ceRNA network consisting of 4 mRNAs, 3 miRNAs, and 30 lncRNAs. Combined with qRT-PCR and prediction results, lncRNA XIST - miR-671-5p - EGFR may play a critical role in the ceRNA network. Our study may shed light on the pathogenesis of periodontitis. However, more researches and trials are still needed to clarify the function of these molecules.