Skip to main content

Epigenetic adaptations of the masticatory mucosa to periodontal inflammation

Abstract

Background

In mucosal barrier interfaces, flexible responses of gene expression to long-term environmental changes allow adaptation and fine-tuning for the balance of host defense and uncontrolled not-resolving inflammation. Epigenetic modifications of the chromatin confer plasticity to the genetic information and give insight into how tissues use the genetic information to adapt to environmental factors. The oral mucosa is particularly exposed to environmental stressors such as a variable microbiota. Likewise, persistent oral inflammation is the most important intrinsic risk factor for the oral inflammatory disease periodontitis and has strong potential to alter DNA-methylation patterns. The aim of the current study was to identify epigenetic changes of the oral masticatory mucosa in response to long-term inflammation that resulted in periodontitis.

Methods and results

Genome-wide CpG methylation of both inflamed and clinically uninflamed solid gingival tissue biopsies of 60 periodontitis cases was analyzed using the Infinium MethylationEPIC BeadChip. We validated and performed cell-type deconvolution for infiltrated immune cells using the EpiDish algorithm. Effect sizes of DMPs in gingival epithelial and fibroblast cells were estimated and adjusted for confounding factors using our recently developed “intercept-method”. In the current EWAS, we identified various genes that showed significantly different methylation between periodontitis-inflamed and uninflamed oral mucosa in periodontitis patients. The strongest differences were observed for genes with roles in wound healing (ROBO2, PTP4A3), cell adhesion (LPXN) and innate immune response (CCL26, DNAJC1, BPI). Enrichment analyses implied a role of epigenetic changes for vesicle trafficking gene sets.

Conclusions

Our results imply specific adaptations of the oral mucosa to a persistent inflammatory environment that involve wound repair, barrier integrity, and innate immune defense.

Background

Genetic studies of chronic inflammatory diseases focused primarily on identifying DNA variants (e.g., single-nucleotide polymorphisms, SNPs) that confer disease risk through genome-wide association studies (GWASs; see [1, 2]). More recently, studies have also examined differences between patients and controls in patterns of DNA methylation (see [3, 4]). Methylation of cytosine is the most common modification of DNA that occurs at CpG dinucleotides (cytosine followed by guanine). It typically causes chromatin condensation and disruption of interactions between DNA and transcription factors, which are associated with transcriptional regulation [5]. Such epigenetic modifications of the chromatin confer plasticity to the genome and enable flexible and reversible responses of the genetic information to environmental challenges, allowing long-term adaptation and fine-tuning of gene expression levels [6]. Correspondingly, epigenetic changes give direct insight into the genetic information that the tissue uses to adapt to environmental factors. A major focus of research into the cellular mechanisms that underlie chronic disease is the tissue-environment interface. Particularly in the mucosa, a balance must be maintained between host defense and uncontrolled, not-resolving inflammation [7, 8]. Maintaining tissue integrity requires continuous adaptations of immune responses, barrier function and regeneration to a changing external and internal environment. It is considered that the development of adverse immune reactions may be associated with these challenges.

The oral masticatory mucosa locates to the alveolus/tooth complex and the hard palate and comprises all parts of the oral mucosa that are involved in chewing. Its position at the entrance to the gastrointestinal tract and respiratory systems entails its direct exposition to a diverse microbiota and various other environmental stressors such as tobacco smoke. If the mucosal barrier is impaired, pathogens can invade into subjacent tissue layers and promote periodontal inflammation. Persisting inflammation causes gingival bleeding and leads to the loss of connective tissue and alveolar bone with subsequent tooth loss, which determines the clinical characteristics of the common oral inflammatory disease periodontitis. In the etiology of periodontitis, persistent gingival inflammation is the most important risk factor [9], followed by smoking [10, 11]. Inflammation and smoking show strong effects in altering DNA-methylation patterns. This implies that in relation to the masticatory mucosa both can be considered as ‘environmental variables’ with a causal role for the onset and progression of periodontitis.

Recently, an epigenome-wide association study (EWAS) to investigate the specific methylation and expression patterns of the healthy masticatory mucosa in the context of cigarette smoke exposure was performed by our group [12]. In this study, differentially methylated positions (DMPs) were re-discovered to be associated with smoking status at a genome-wide significance level within the genes CYP1B1 (cytochrome P450 family 1 subfamily B member 1) and AHRR (aryl-hydrocarbon receptor repressor). Several EWAS have identified these associations for cells of the alveolar [13] and buccal mucosa [14] before, putting emphasis on the role of epigenetic adaptations of the activity of these genes in barrier tissues that are long-term exposed to tobacco smoke metabolites. Some of these associations were restricted to solid epithelial tissues and were not found in blood, indicating the tissue-specificity of CpG methylation patterns relating to tissue function.

To date, few studies were published on differential DNA methylation in affected solid tissues in the course of inflammatory disease, where differences in cell-type composition are disease-immanent, as cell-type deconvolution algorithms still lack specific reference panels for many cell types. EpiDish provides an accurate reference dataset for generic epithelial tissue [15]. In the present study, we aim to confirm its suitability for the detection of cells specific for gingival tissue and employ it for cell-type deconvolution of our data. As differences in cell-type composition represent a major confounding factor in EWAS of inflammatory diseases, they also hinder the interpretation of effect sizes in affected tissue, which may be crucial for evaluating the relevance of significant associations. By applying our recently published “intercept-method” for the inference of effect sizes adjusted for confounding factors [16] to our data, we provide adapted estimations of the effect sizes at DMPs, which allows us to shed light on the molecular mechanisms altered in the etiology of periodontitis.

The aim of the current study was to investigate epigenetic changes in the gingiva in the course of long-term inflammation that resulted in periodontitis. To this end, the methylation patterns of inflamed and clinically uninflamed gingival tissue biopsies of 60 periodontitis cases were analyzed. To our knowledge, this is the first EWAS that investigated genome-wide epigenetic changes in solid biopsies of the gingiva in response to oral inflammation with adjustment for cell-type heterogeneities.

Results

Pre-processing pipeline

786,547 probes passed the quality control (QC) criteria and were analyzed in ex vivo gingival tissue biopsies from 60 periodontitis patients obtained from an inflamed and an uninflamed site each. Additionally, DNA methylation patterns of cultured gingival epithelial cells (GECs) and human gingival fibroblasts (HGFs), collected from 5 and 4 additional donors, respectively, and of immortalized cell lines of gingival epithelial (OKG4) and fibroblast (ihGF) cells, were analyzed.

Suitability test of the EpiDish algorithm for the detection of gingival epithelial cells and gingival fibroblasts for cell type deconvolution

The cell type deconvolution algorithm EpiDish represents an established method to adjust for cell type heterogeneity in EWAS samples. It uses reference methylomes from various cell types. To test how the implemented reference methylomes of immune cells, epithelial cells, and fibroblasts matched the methylome of cells from gingival tissue, we analyzed the DNA from GECs and HGFs from cell culture on the EPIC BeadChip and performed cell type deconvolution with EpiDish. The primary GECs and both immortalized cell lines showed 100% concordance with the reference datasets (Fig. 1, Additional file 1). Primary gingival fibroblasts showed some deviations from the expected patterns, with 3 of the 4 samples showing an estimated fibroblast fraction of > 94%. One sample implied only 88% fibroblasts, with an infiltration of 11% immune cells and 2% epithelial cells.

Fig. 1
figure 1

EpiDish results for gingival cell cultures. Cell type estimations in the cultured gingival cells inferred by EpiDish as average weight proportions of the major cell types epithelial cells (Epi), fibroblasts (Fib,) and immune cells (IC) for immortalized gingival epithelial cells (A) (n = 1), immortalized gingival fibroblasts (B) (n = 1), primary gingival epithelial cells (C) (n = 5), and primary gingival fibroblasts (D) (n = 4)

Estimation of the immune cell fraction in inflamed and uninflamed gingiva.

Substantial infiltration of immune cells into gingival tissues is disease-immanent in periodontitis. To account for differences in cell type composition, we applied the EpiDish algorithm to our 120 samples of paired inflamed an uninflamed ex vivo biopsies from 60 individuals. In the inflamed samples, half of the cell fraction of inflamed gingival tissues consisted of immune cells (mean 0.52, standard deviation (SD) 0.18). Some of the inflamed biopsies showed estimated immune cell fractions as high as 0.91 (Additional file 1). Uninflamed biopsies showed immune cell infiltration that was significantly lower compared to inflamed tissue (mean 0.28, SD 0.15; p < 10–5, Wilcoxon Test). In both inflamed and uninflamed samples, the amount of estimated immune cells was inversely correlated to the amount of epithelial cells in the samples (Pearson’s correlation coefficient = − 0.93). Estimations for epithelial cells accordingly were higher in uninflamed samples (mean 0.52, SD 0.15) than in inflamed samples (mean 0.32, SD 0.18). The degree of variation of immune cell infiltration was comparable between uninflamed and inflamed tissues, as was the degree of variation of the amount of epithelial cells. For 6 individuals, we observed higher immune cell fractions in clinically uninflamed samples compared to inflamed samples. The variation of the estimated amount of fibroblast cells in inflamed and uninflamed biopsies was lower and comparable between inflamed and uninflamed biopsies, with an estimated fibroblast fraction of 0.20 (SD 0.07) in uninflamed samples and 0.17 (SD 0.08) in inflamed samples.

We furthermore performed a multidimensional scaling analysis to identify the effects of potential confounders on the methylome of gingival biopsies. This analysis revealed a separation according to the amount of immune cells as inferred by EpiDish, leading to the location of some biopsies that were clinically diagnosed as uninflamed in the cluster of inflamed biopsies (Fig. 2A, B). No clustering patterns were associated with the site of biopsy extraction (Additional file 2).

Fig. 2
figure 2

Multidimensional Scaling (MDS) Plots for inflamed and clinically uninflamed samples. A MDS Plot of the 120 samples from the 60 patients initially analyzed for cell-type heterogeneity using EpiDish, colored by diagnosed inflammation status. B MDS Plot of the 120 samples from the 60 patients initially analyzed for cell-type heterogeneity using EpiDish, colored according to the amount of infiltrated immune cells. C MDS Plot of the EpiDish-filtered 96 samples used for the EWAS, colored by diagnosed inflammation status. D MDS Plot of the EpiDish-filtered 96 samples used for the EWAS, colored according to the amount of infiltrated immune cells. In the initial sample panel, clinically uninflamed biopsies with high estimations of infiltrated immune cells clustered together with tissue diagnosed as inflamed (A, B). When including only samples diagnosed as “uninflamed” and “inflamed” that contained ≤ 35% and ≥ 40% immune cells, respectively, inflammation status (“uninflamed” and “inflamed”) clustered distinctly, with immune cell fractions as estimated by EpiDish reflecting this classification (C, D)

Based on the estimations of the cell type deconvolution, we excluded biopsies that were clinically diagnosed as uninflamed but appeared to be infiltrated by high amounts of immune cells, and those that were clinically diagnosed as inflamed but the estimated immune cell fraction was very low. Exclusion of these biopsies reduced confounding effects from putative misclassification or other factors. We set threshold criteria for sample inclusion as follows: samples diagnosed as “uninflamed” and “inflamed” contained ≤ 35% and ≥ 40% immune cells, respectively. In total, 48 samples from clinically uninflamed sites and 48 samples from inflamed sites complied with this criterion and were selected for subsequent association analyses. These 96 samples originated from 57 individuals, 39 of whom donated an uninflamed and an inflamed biopsy each. Thereafter, the mean estimation of immune cells was 0.22 (SD 0.06) for uninflamed and 0.58 (SD 0.13) for inflamed biopsies in the remaining 96 samples (Table 1, Fig. 3).

Table 1 Basic characteristics of the EWAS study population
Fig. 3
figure 3

EpiDish results for the 96 samples selected for EWAS analysis. Cell type estimations in the gingiva inferred by EpiDish as average weight proportions of the major cell types epithelial cells (Epi), fibroblasts (Fib), and immune cells (IC) for uninflamed (A) and inflamed samples (B) (n = 48 each)

After removing individuals with immune cell fractions that contradicted their clinical classification from the analysis, samples clustered according to their clinical classification into inflamed and uninflamed tissue (Fig. 2C, D).

DNA methylation differences of uninflamed and inflamed gingival biopsies

Next, the selected samples were investigated for significant changes in methylation patterns between clinically uninflamed and inflamed samples in the 786,547 CpG sites that passed QC, with adjustment for batch effects and differences in immune cell content. Furthermore, as smoking is one of the major risk factors for periodontitis, 1501 DMPs associated with smoking in oral epithelial cells [14] were removed from the analysis. The quantile–quantile plot revealed global inflation of the test statistics compared to the expected distribution, with an inflation factor of λ = 5.97 (Fig. 4A). After correction for multiple testing, we found 15,507 DMPs with q < 0.05 (Figs. 4B, 5).

Fig. 4
figure 4

Manhattan and quantile–quantile plot for epigenome-wide associations with inflammation in the gingiva. A The quantile–quantile plot showed evidence for inflation of association signals (λ = 5.97). B Manhattan plot showing −log10 transformed p values plotted against the genomic location of the probes. The horizontal lines indicate the genome-wide significance threshold (p < 10–7)

Fig. 5
figure 5

Volcano plot showing methylation differences of clinically uninflamed compared to inflamed samples against −log10 of p values. Depicted in red are the DMPs in periodontal inflammation significant after adjustment for multiple testing. Note that effect sizes were derived from M-values using our intercept-method [16], taking confounding factors, especially cell type heterogeneity, into account

Intra-individual variation in DNA methylation is a strong confounder in epigenetic studies, which can be precluded by analyzing clinically uninflamed and inflamed biopsies of the same patients. Likewise, we observed that of the top 10 associated DMPs (best q = 4.5 × 10–16 at cg19478962, locating to the noncoding RNA gene LOC643339), only three showed nominally significant associations in a sub-analysis of the 39 individuals for whom an uninflamed and inflamed sample each was available for the EWAS. Considering this and the genomic inflation factor of 5.97, a filtering strategy was applied on the associated DMPs to minimize the possibilities of false positive associations due to inter-individual variation.

DMPs in the full EWAS panel and a subset of paired inflamed and clinically uninflamed samples from the same individuals

We tested all 15,507 DMPs that were significant in the full sample panel of 57 EpiDish-filtered individuals (48 uninflamed and 48 inflamed samples, “full analysis”), which gave the largest power and thus reduced the potential of false negative findings, in the subset of paired samples (both inflamed and uninflamed samples from 39 individuals, “paired analysis”), to further reduce the potential of false positive findings. In this sub-analysis, 2347 DMPs showed significant associations with periodontal inflammation in the gingiva after Bonferroni-correction for 15,507 tests (Additional file 3). The 20 most significant DMPs showed p values < 10–15 and adjusted ∆β ≥ 0.05 in the full analysis (Table 2). The most significant DMP (cg23278359, Fig. 6) mapped to the gene PTP4A3 (Protein Tyrosine Phosphatase 4A3) with p = 2.2 × 10–18, q = 1.8 × 10–12, and ∆β = − 0.18. The second most significant DMP mapped to the gene ROBO2 (Roundabout Guidance Receptor 2; cg17282085; q = 4.9 × 10–12, ∆β = − 0.14), which is a suggestive risk gene for periodontitis (age 20–60 years; rs264537-C, p = 3 × 10–6, odds ratio = 1.35 [95% confidence interval 1.19–1.54]) [17]. Of the 2347 DMPs significant in the full and the paired analysis, 63 showed adjusted effect sizes > 0.3, with the 20 most pronounced effect sizes ≥ 0.34 (Table 3). The DMP with the highest effect size mapped to the gene CCL26 (C–C Motif Chemokine Ligand 26; cg11303839, q = 4.0 × 10–6, ∆β = − 0.43). 3 of the top 20 associated DMPs and 2 of the 20 DMPs with the highest effect sizes located < 2 kb to other DMPs with q < 0.05 in the full analysis.

Table 2 Most significant DMPs in full analysis that were also significant in the paired sub-analysis
Fig. 6
figure 6

Methylation at top associated loci. Raw beta values of uninflamed and inflamed EWAS samples (full set, n = 48 each) at cg17282085 (ROBO2), cg11303839 (CCL26), cg23278359 (PTP4A3), cg00320534 (DNAJC1), cg12891342 (LPXN), and cg14991316 (BPI). IC = immune cell fraction in samples

Table 3. 20 DMPs with highest effect sizes

Subsequently, a more stringent filtering process for the 15,507 DMPs that were significant in the full analysis was applied to further reduce the number of putative random associations. These filtering criteria were: (1) ≥ 2 significant DMPs (q < 0.05 in the full analysis) had a maximum distance of 2 kb, (2) ≥ 1 DMP of such a cluster showed effect sizes > 0.1 in the full analysis, (3) ≥ 1 DMP of such a cluster showed padj < 0.05 in the paired analysis. This resulted in 441 DMPs in 193 clusters, with p values < 6.5 × 10–8 in the full analysis (Additional file 4). Of these, 22 DMPs showed padj < 10–5 in the paired analysis (Table 4).

Table 4 DMPs associated with padj < 10–5 in the paired sub-analysis after stringent filtering

The most significant DMP in the paired analysis (cg00320534, adjppaired = 9.1 × 10–6, ∆β = − 0.2, 1 supporting DMP) located to the gene DNAJC1 (DnaJ Heat Shock Protein Family Member C1). The second most significant DMP located to the gene LPXN (Leupaxin) and was supported by 5 additional DMPs that showed q-values < 0.05 in the full analysis (cg12891342, adjppaired = 1.4 × 10–5, ∆β = − 0.2). Subsequently, we ranked the 441 DMPs according to their effect size delta beta (in the full analysis). 7 DMPs showed effect sizes > 0.3 (Table 5). The DMP with the largest effect size from this filtering approach (cg14991316, adjppaired = 1.0 × 10–3, ∆β = − 0.34, 3 supporting DMPs) mapped to the promoter region of BPI (Bactericidal Permeability Increasing Protein).

Table 5 DMPs with effect sizes ≥ 0.3 after stringent filtering

Identification of differentially methylated regions

To identify differentially methylated regions (DMR), we applied the Bump Hunter algorithm [18] to our full dataset of 48 inflamed and 48 uninflamed samples and to the subset of paired samples. In the full dataset, 6 DMR showed a family-wise error rate (FWER) < 0.05, whereas in the paired sub-analysis, 11 DMR showed a FWER < 0.05 (Table 6). 6 DMR were significantly associated in both approaches, with a 151 base pair (bp) region at the transcription start site of CNIH4 (Cornichon Family AMPA Receptor Auxiliary Protein 4) showing the strongest associations (Fig. 7).

Table 6 Summary of significant differentially methylated regions
Fig. 7
figure 7

Methylation at significant differentially methylated regions. Raw beta values of uninflamed and inflamed EWAS samples (full set, n = 48 each) at the genetic regions of CNIH4, LMNB2, NFKBIL1, PRAP1, PRPTN2, and the tRNA pseudogene on chromosome 16

Gene ontology and gene set enrichment analysis

To determine whether the observed differential methylation is enriched in specific gene ontology (GO) terms or Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, we performed an enrichment analysis for the 2347 DMPs that were significantly associated in the full and the paired analysis. No GO term or KEGG pathway showed significant overrepresentation after adjustment for multiple testing (Table 7, Additional file 5). The most significant term (“inactivation of X chromosome by genetic imprinting”, punadj = 4.3 × 10–5) includes only three genes, PCGF3 (Polycomb Group Ring Finger 3), PCGF5 (Polycomb Group Ring Finger 5), and PCGF6 (Polycomb Group Ring Finger 6), which all show significant DMPs in the EWAS and expression in healthy keratinized oral mucosa (transcripts per million (TPM) > 100), and are involved in chromatin modeling.

Table 7 Gene Ontology enrichment analysis

Furthermore, when considering all GO pathways with nominal p values < 0.001, which might indicate a trend, there was an overrepresentation of GO terms corresponding to nervous system-related terms. As an example, we investigated the associations of the second most significant GO term, “synaptic vesicle cycle” (biological process, p = 9.0 × 10–5), more in detail, because this term included more genes compared to the most enriched GO term. 25 genes belonging to this GO term showed DMPs in the full and paired analysis, 10 of which with a q < 0.0001 in the full set (Table 8). From the top 10 differentially methylated genes in this GO term, STRING [20] indicated an interaction of 4 genes, SNAP23 (Synaptosome Associated Protein 23), SYT9 (Synaptotagmin 9), STX1B (Syntaxin 1B), and SYT2 (Synaptotagmin 2) with the periodontitis-associated genes VAMP3 [21, 22] and VAMP8 (Vesicle Associated Membrane Protein 8; Additional file 6) [23]. All of these genes show relatively little to no expression in the healthy keratinized oral mucosa, determined by RNA-Sequencing of 39 samples [12], with exception of CLCN3 and SNAP23 (TPM = 2771 and 615, respectively). Notably, the location of the DMPs within these genes are predicted promotor flanking and transcribed regions as indicated by ENCODE [24], whereas the other, less expressed genes have significant DMPs within regions mapping to sites mainly predicted as repressed or low activity regions. Only 4 KEGG pathways showed nominal significance (p < 0.05), among which was “SNARE interactions in vesicular transport” with p = 0.036 (3rd most significant pathway, Additional file 5). Notably, the functionally related pathway “Mucin type O-glycan biosynthesis” was the 8th most significant pathway out of 343 with p = 0.09.

Table 8 Differentially methylated genes within the gene ontology term “synaptic vesicle cycle” (biological process)

Discussion

In this study, we performed the first cell-type informed EWAS of gingival tissue biopsies to identify CpG methylation differences between the uninflamed and periodontitis-inflamed oral mucosa. EWAS that investigate environment-induced epigenetic changes in tissues of multiple cell types need to address cell type deconvolution and intra-individual variation of methylation patterns. Likewise, we observed substantial immune cell infiltration within the inflamed gingiva leading to substantial effects on the methylome of the oral mucosa, according to which inflamed gingival tissue contained approximately twice as much infiltrated immune cells as clinically uninflamed tissue, with 52% in the total of 60 inflamed samples.

Surprisingly, we furthermore observed considerable amounts of immune cells in samples excised from periodontitis patients at sites diagnosed as uninflamed by specialized periodontologists, sometimes as high as in inflamed samples, with an average of 28% in the total of 60 clinically uninflamed biopsies sampled in this study. Compared to this, in a previous study, we identified immune cell proportions of 16% in oral masticatory mucosa samples from healthy individuals with no diagnosed oral inflammation, taken from the hard palate near the gingival margin [12]. It is conceivable that periodontitis patients show a general immune response in the affected tissue that leads to an increased infiltration of immune cells also at sites where no periodontal tissue inflammation was diagnosed. This underlines the importance of incorporating cell-type deconvolution algorithms in differential expression and methylation studies to avoid type I and type II errors.

Another limitation of EWAS is the complexity of intra-individual variation of methylation patterns, which has a stronger impact if sample size is limited. As a consequence, EWAS have higher genomic inflation factors compared to GWAS, where inter-individual genetic variation generally plays a much lesser role. Accordingly, the current study had a genomic inflation factor of 5.97 despite cell type deconvolution and removal of samples with ambiguous inflammatory data. We aimed to reduce intra-individual methylation variation by analyzing paired biopsies from 60 individuals on the EPIC BeadChip. However, after removing samples with an ambiguous immune state, our study comprised a non-paired analysis panel of 48 clinically uninflamed and 48 periodontitis-inflamed samples from 57 patients in total. Of these, 39 patients both donated an uninflamed and an inflamed biopsy. To increase statistical power and thus, to reduce the potential of false negative findings, the total sample panel of 2 × 48 non-paired samples was analyzed first. Subsequently, a sub-analysis of the significant results in the 2 × 39 paired samples was performed to remove potential false positive findings. This two-step approach in conjunction with cell type deconvolution, an improved intercept method for the estimation of confounder-informed effect sizes, recently developed by our group [16], and stringent filter criteria allowed us to identify several genes that showed robust different methylation between uninflamed and periodontitis-inflamed gingival tissues, independent of immune cell infiltration and inter-individual variation. Some of these genes have plausible roles in the etiology of periodontitis.

The most significant DMP with adjusted effect sizes mapped to the last intron of PTP4A3 (Protein Tyrosine Phosphatase 4A3). PTP4A3 has a role in the positive regulation of endothelium wound repair and angiogenesis [25] and inhibition of the expression of extracellular matrix and adhesion genes, like matrix metalloproteinases (MMPs) and integrins [26, 27]. Substrates of PTP4A3 are e.g. MMP14 [28], integrin β1 [27], and Keratin 8 (KRT8) [29]. The identification of PTP4A3 is a good example of the suitability of our intercept method to provide effect sizes that directly correspond to the p value, which is missing in most EWAS. Currently, p values are usually based on log-transformed M-values because they exhibit better statistical properties, but are biologically meaningless. As a consequence, p values based on M-values are reported together with differences in raw β-values, which do not account for confounding factors. At the PTP4A3 locus, such an approach would result in a highly significant p value with a ∆β of − 0.07 (calculated from the raw β = 0.87 in clinically uninflamed and 0.8 in inflamed tissue). The mean raw β at this locus is 0.8 in cultured GECs and 0.93 in blood cells. As the inflamed biopsies from the EWAS show an average amount of ~ 58% immune cells, 58% of the methylation signals in these samples are derived from cells that presumably exhibit methylation values as high as β = 0.93. Given the periodontal inflammation itself has no impact on methylation here, estimations for these highly admixed samples should be higher than for “pure” epithelial cells (i.e., the GECs). Instead, based on our intercept method, an adjusted ∆β of − 0.18 was estimated, indicating a strong hypomethylation of the gingival epithelial cells at this locus in inflamed samples.

The second most significant DMP mapped to the exonic region of ROBO2 (Roundabout Guidance Receptor 2) that encodes a transmembrane receptor of the immunoglobulin superfamily and has been reported as a suggestive risk gene for periodontitis before [17]. It is involved in epithelial wound repair by promoting extrusion of dying cells from injured tissue [30].

The highest adjusted effect size with a hypomethylation of 43% was observed in the gene CCL26 (C–C Motif Chemokine Ligand 26). CCL26 is a factor for chemotactic eosinophil and basophil attraction in endothelial cells and possesses antimicrobial activity [31, 32]. Notably, CCL26 was the most highly induced gene in patients suffering from eosinophilic esophagitis, a chronic allergic inflammation of the esophageal mucosa, compared with its expression level in healthy individuals [33].

Among the 20 genetic regions with the highest effect sizes was PTP4A2 (Protein Tyrosine Phosphatase 4A2). This gene forms, together with PTP4A3, which harbored the most significant DMP in the full analysis as discussed above, and another phosphatase, PTP4A1, the subfamily of Phosphatase of Regenerative Liver (PRL). Given that two out of three PRL subfamily members showed highly significant DMPs in periodontitis-inflamed gingiva, the PRL subfamily represents a good candidate for evaluating its role in the inflammatory processes of oral epithelial cells.

In a more stringent filtering process, only “clusters” of DMPs that comprised at least 2 DMPs within 2 kb in the full analysis, at least one DMP with an effect size > 0.1, and one DMP significant in the paired analysis were included. In this way, further genes with a suggestive role in periodontal inflammation of the gingiva were identified. The most significant gene in the paired analysis, DNAJC1 (DnaJ Heat Shock Protein Family Member C1), is an endoplasmatic reticulum heat shock protein that binds the molecular chaperone HSPA5 (alias BiP) [34], which is known to exhibit anti-inflammatory and inflammation-resolutory properties. BiP was proposed to form part of “resolution-associated molecular patterns” (RAMPs), intracellular responses to inflammatory and damage-induced stress, mediated by pathogen-associated (PAMPs) and damage-associated (DAMPs) molecular patterns [35]. Furthermore, BiP interacts with SERPINA3 [36], an acute phase protease inhibitor, which is supposed to inhibit angiotensin activation [18]. The second most significant DMP mapped to Leupaxin (LPXN), a focal adhesion protein, which showed a cluster of 6 DMPs. It is preferably expressed in hematopoietic cells with a putative role as an adapter protein in the formation of the adhesion zone in osteoclasts (Uniprot) and functions as a paxillin counterpart that potently suppresses the tyrosine phosphorylation of paxillin during integrin signaling [37].

The DMP with the highest effect sizes from the stringent filtering approach mapped to BPI (Bactericidal Permeability Increasing Protein) and was also among the top 20 significant DMPs (paired analysis) in the stringent filtering approach. BPI is a lipopolysaccharide-binding protein with strong antimicrobial activity and an important part of innate immune response. It is predominantly expressed in hematopoietic cells, but is also found in a variety of other tissues, most notably the epithelial lining of mucous epithelial cells, where it was shown to block endotoxin-mediated signaling and to kill Salmonella typhimurium [38].

Enrichment analysis for the 2347 DMPs that were significant in the full as well as in the paired analysis showed no significant enriched GO terms or KEGG pathways after correction for multiple testing. However, GO terms with nominal p values < 0.001 were enriched for vesicle trafficking. Specifically, while the strongest enriched GO term, “inactivation of X chromosome by genetic imprinting”, encompassed only 3 genes, for the second most significant GO term, “synaptic vesicle cycle”, we found 25 out of 229 genes to be differentially methylated in periodontitis-inflamed gingiva. Likewise, the KEGG pathway “SNARE interactions in vesicular transport” was the 3rd strongest nominally significantly enriched pathway (p = 0.036) and the functionally related pathway “Mucin type O-glycan biosynthesis” was the 8th most enriched pathway out of 343 (p = 0.09). Interestingly, several genes among the most significant differentially methylated genes that were part of the vesicle trafficking gene sets are known to interact with the genes VAMP3 and VAMP8, which are both suggestive risk genes of periodontitis [21,22,23]. VAMP8 and VAMP3 are highly expressed in keratinized oral mucosa, with TPM values of 2404 and 1839, respectively [12], indicating a functional role of these genes in epithelial and connective oral tissues. For example, in colonic epithelial cells, VAMP8 is required for the secretion of the main mucin of the mucin layer of the colon, Mucin-2, which has an important role in the maintenance of the mucosal barrier integrity. However, the observed trend of gene sets enriched for pathways of vesicle trafficking and mucin biosynthesis was not significant after correction for multiple testing, why these data need to be carefully interpreted.

In this context, a limitation of the present EWAS was the small sample size that impeded clear statistical evidence for enriched gene sets and a putatively more comprehensive set of DMPs. However, to our knowledge, it currently represents the largest EWAS that used gingival tissue biopsies and the only EWAS that analyzed paired inflamed and uninflamed samples from the same patients. Two previous EWAS that investigated differential methylation of periodontitis-inflamed gingiva [39, 40] included only 19 and 12 patients, respectively, and similar numbers of healthy controls. Moreover, the previous studies did not adjust for differences in cell type composition. Given the substantial immune cell infiltration, it is questionable whether findings from studies that ignore cell type deconvolution approaches point towards differential methylation patterns that reflect disease-relevant cellular processes or rather reflect differences in cell type composition.

In conclusion, this EWAS identified several genes that are significantly differentially methylated between periodontitis-inflamed compared to uninflamed gingiva. The strongest differences were observed for genes with roles in wound healing (ROBO2, PTP4A3), cell adhesion (LPXN) and innate immune response (CCL26, DNAJC1, BPI). Functional enrichment analyses implied that differentially methylated genes were overrepresented in gene sets annotated to vesicle trafficking. These results propose that the oral mucosa responds to a long-term inflammatory environment with specific adaptations in wound repair, barrier integrity, and innate host defense.

Methods

Study sample

A total of 60 periodontitis patients of Caucasian descent, aged 24–65, with ≥ 30% bone loss from the cement-enamel junction to the root apex, documented by a full-mouth set of radiographs or orthopantomographs, at ≥ 2 teeth were enrolled in this study. Written informed consent was obtained from all subjects according to the Declaration of Helsinki. All participants completed a detailed questionnaire to provide general personal information (e.g., sex, age, geographical and family descent), information on general and oral health, and smoking habits. This study was conducted in accordance with the principles of Good Clinical Practice and approved by the Independent Ethics Committee of each participating University Hospital.

Collection of ex vivo tissue samples from clinically uninflamed and inflamed gingiva

Clinically uninflamed gingival tissue samples were collected during routine tooth extraction, surgical tooth lengthening, reopening of a submerged healing implant as anyhow performed during the intervention, with a tissue puncher with a diameter of 3 mm or from the tissue margins with a scalpel. A tissue sample of periodontitis-inflamed gingiva was taken from excised tissue due to modified Widman flap, osseous resective surgery [41], or with a disposable 3 mm gingival tissue puncher. To stabilize DNA and RNA, the biopsies were stored in the AllProtect reagent (Qiagen, Hilden, Germany) immediately after punching at 4 °C for 24 h, and subsequently transferred to − 20 °C.

Cell culture of gingival epithelial cells and gingival fibroblasts

Solid ex vivo biopsies of the masticatory oral mucosa of the hard palate were taken from additional healthy subjects of Caucasian descent with 3 mm tissue punchers. Biopsies were dissected enzymatically to separate the epithelial cells from the fibroblasts by overnight incubation in 5 mg/mL dispase II (Sigma Aldrich) diluted in cell growth medium (DMEM, 1% Pen/Strep) at 4 °C. Immortalized and primary gingival fibroblast cells (iHGFs and pHGFs) were cultured in DMEM cell growth medium with 1% Amphotericin B, 1% Pen/Strep, and 1% non-essential amino acids. Immortalized gingival epithelial cells (OKG4) [42, 43], kindly provided by James Rheinwald, Boston, MA, USA, and primary GECs were cultured in DermaLife K medium with 1% penicillin/streptomycin in the presence of 60 μM or 1.4 mM Ca2+ (CellSystems, Troisdorf,Germany). Gingival fibroblasts and epithelial cells were cultivated in collagen-coated flasks at 37 °C and 5% CO2 until reaching 100% and 80% confluence, respectively, before DNA extraction.

DNA extraction

The conserved ex vivo tissue samples were manually broken up into small pieces with a scalpel and subsequently homogenized using the automated mixer mill MM 400 (Retsch, Haan, Germany) with frozen beads (3 mm, Retsch) for 90 s at 30 Hz before further processing. DNA was extracted using the AllPrep DNA/RNA/miRNA Universal Kit (Qiagen, Hilden, Germany). DNA integrity was subsequently verified with a 2% agarose gel electrophoresis. Concentrations were measured using the Multiskan GO Microplate Spectrophotometer (Fisher Scientific, Hampton, USA). DNA samples were stored at − 80 °C.

Bisulfite conversion and hybridisation to the Infinium MethylationEPIC BeadChips

500 ng DNA per sample was bisulfite converted with the EZ-96 DNA Methylation Kit (Zymo Research, Irvine, USA) and hybridized to the Infinium DNA MethylationEPIC BeadChip (Illumina, San Diego, USA) on an iScan Microarray Scanner (Illumina) at the Institute for Clinical Molecular Biology, Christian-Albrechts-University Kiel, Germany.

Pre-processing and normalization

For analysis and quality control, the R environment (Version 4.0.3) and the R package minfi (Version 1.36.0) [44] were used, including the R package limma (Version 3.36.0) for the differential methylation analysis [45]. The Red/Green channel of the array were converted into one methylation signal without any normalisation. Using the function plotQC, sample-specific QC was estimated and no sample was removed. The QC criteria for probes were filtering of probes with median detection p values > 0.05, probes that lay within 5 bp of a single nucleotide polymorphism (SNP) with > 5% minor allele frequency, probes on the sex chromosomes and cross-reactive probes using the R package maxprobes (https://rdrr.io/github/markgene/maxprobes) according to Pidsley et al. [46] and McCartney et al. [47]. In total, 786,547 probes complied with the QC criteria and were included in the analysis. Methylation status was estimated according to the fluorescent intensity ratio, as any value between β = 0 (unmethylated) and 1 (completely methylated), and log2-transformed into M-values, which are considered a more statistically valid estimator [48]. After quality control, a functional normalisation was performed using the preprocessFunnorm function in minfi [49].

Cell type deconvolution of EWAS samples

To identify the presence of non-epithelial cells in the oral mucosa samples, the EpiDISH algorithm (Version 2.6.0) [15] was used, applying the centEpiFibIC.m.rda reference dataset on our raw beta-values, which includes centroids for epithelial cells, fibroblasts, and immune cells (IC).

Differential DNA methylation analysis

To identify DMPs correlating significantly with periodontal inflammation of the gingiva, we performed a linear mixed model (R package lme4, Version 1.1-27.1) [50] on the approximatively normally distributed M-values. We included affection status (inflamed vs. uninflamed) and IC-content (derived from the EpiDish algorithm) as fixed effects and the batch effect, i.e., slide and position on slide, as random effect with a constant intercept, using the R packages Combat and BatchQC [51]. Therefore, we adjusted the effect of the affection status, i.e., the delta M, for the confounder IC-content. Afterwards, we transformed the differences in M-values into differences in β-values by our recently developed intercept method [16]. By doing so, we were able to include the batch effect as well as the confounding effect of immune cell infiltration in gingival tissue.

To exclude potential stratification by the effects of smoking, which represents a major risk factor for periodontitis, 1501 probes that corresponded to DMPs in buccal mucosa of smokers were removed from the analysis beforehand [14]. Correction for multiple testing was carried out using the method by Benjamini and Hochberg [52]. CpGs were annotated to genes according to GRCh37/hg19 as provided in the MethylationEPIC BeadChip manifest.

In general, reported effect sizes Δβ were adjusted for IC-content and batch effects as described above, unless otherwise stated (i.e., referred to as “raw beta” values). To compare raw beta values of gingival cells with raw beta values of hematopoietic cells at significant DMPs, Illumina EPIC BeadChip data of 69 blood samples taken from 34 healthy female individuals at two different time points, was downloaded from the NCBI GEO database [53], accession GSE123914 [54].

Stringent filtering approach

A more stringent filtering approach was performed for DMPs that were significant in the full analysis. As methylation is spatially correlated [55], next to the information on adjusted effect sizes, we incorporated information on nearby co-methylation in this filtering approach.The applied filtering criteria were: (1) ≥ 2 significant DMPs (q < 0.05 in the full analysis) had a maximum distance of 2 kb, (2) ≥ 1 DMP of such a cluster showed effect sizes > 0.1 in the full analysis, (3) ≥ 1 DMP of such a cluster showed padj < 0.05 in the paired analysis. We chose a maximum gap between DMPs of 2 kb as co-methylation was found to occur within this distance [56] and allowed for DMPs only associated in the full analysis as “supporting DMPs” to circumvent the problem of false negative findings in the paired analysis (which had less statistical power due to its limited size).

Identification of differentially methylated regions

In addition to the stringent filtering process described above, which filtered also according to a spatial pattern of DMPs, we aimed to identify differentially methylated regions (DMR) using a distinct straightforward algorithm, tailored to find regions where methylation patterns deviate from the expected distribution. While this approach cannot account for our specific experimental setup with a subset of paired samples, it provides a high accuracy in finding true positive associations in “standard” case control settings, also controlling for confounding variables. To this end, we applied the Bump Hunter package, Version 1.32.0 [18], to our full dataset of 48 inflamed and 48 uninflamed samples and to the subset of paired samples. The cut-off value, which is a user-defined numeric value that determines the upper and lower bounds of the genomic profiles that are used as candidate regions, was set to 0.02 and the number of permutations was set to 1000.

Gene ontology and gene set enrichment analysis of differentially methylated genes

For gene set and gene ontology enrichment analysis, we used the r packages missMethl and gometh [19]. The R package IlluminaHumanMethylationEPICanno.ilm10b4.hg19 (Version 3.13) was used for annotation. 344 KEGG terms were gathered using limma’s internal function getGeneKEGGLinks with default options for gene set enrichment analysis. 22,746 GO terms were derived fromorg.Hs.eg.db after cleaning. Gene ontology terms and gene sets with q < 0.05 were considered as being significantly enriched.

Availability of data and materials

Methylation data from the EWAS experiment will be available from the corresponding author on reasonable request.

References

  1. 1.

    Buniello A, MacArthur JAL, Cerezo M, Harris LW, Hayhurst J, Malangone C, et al. The NHGRI-EBI GWAS Catalog of published genome-wide association studies, targeted arrays and summary statistics 2019. Nucleic Acids Res. 2019;47(D1):D1005–12.

    CAS  PubMed  Article  Google Scholar 

  2. 2.

    Visscher PM, Wray NR, Zhang Q, Sklar P, McCarthy MI, Brown MA, et al. 10 years of GWAS discovery: biology, function, and translation. Am J Hum Genet. 2017;101(1):5–22.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  3. 3.

    Flanagan JM. Epigenome-wide association studies (EWAS): past, present, and future. Methods Mol Biol. 2015;1238:51–63.

    PubMed  Article  PubMed Central  Google Scholar 

  4. 4.

    Rakyan VK, Down TA, Balding DJ, Beck S. Epigenome-wide association studies for common human diseases. Nat Rev Genet. 2011;12(8):529–41.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  5. 5.

    Rothbart SB, Strahl BD. Interpreting the language of histone and DNA modifications. Biochim Biophys Acta. 2014;1839(8):627–43.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  6. 6.

    Bayarsaihan D. Epigenetic mechanisms in inflammation. J Dent Res. 2011;90(1):9–17.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  7. 7.

    Meyle J, Chapple I. Molecular aspects of the pathogenesis of periodontitis. Periodontol 2000. 2015;69(1):7–17.

    PubMed  Article  PubMed Central  Google Scholar 

  8. 8.

    Hajishengallis G, Darveau RP, Curtis MA. The keystone-pathogen hypothesis. Nat Rev Microbiol. 2012;10(10):717–25.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  9. 9.

    Lang NP, Schatzle MA, Loe H. Gingivitis as a risk factor in periodontal disease. J Clin Periodontol. 2009;36(Suppl 10):3–8.

    PubMed  Article  PubMed Central  Google Scholar 

  10. 10.

    Eke PI, Dye BA, Wei L, Thornton-Evans GO, Genco RJ. Prevalence of periodontitis in adults in the United States: 2009 and 2010. J Dent Res. 2012;91(10):914–20.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  11. 11.

    Eke PI, Dye BA, Wei L, Slade GD, Thornton-Evans GO, Borgnakke WS, et al. Update on prevalence of periodontitis in adults in the United States: NHANES 2009 to 2012. J Periodontol. 2015;86(5):611–22.

    PubMed  PubMed Central  Article  Google Scholar 

  12. 12.

    Richter GM, Kruppa J, Munz M, Wiehe R, Hasler R, Franke A, et al. A combined epigenome- and transcriptome-wide association study of the oral masticatory mucosa assigns CYP1B1 a central role for epithelial health in smokers. Clin Epigenetics. 2019;11(1):105.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  13. 13.

    Stueve TR, Li WQ, Shi J, Marconett CN, Zhang T, Yang C, et al. Epigenome-wide analysis of DNA methylation in lung tissue shows concordance with blood studies and identifies tobacco smoke-inducible enhancers. Hum Mol Genet. 2017;26(15):3014–27.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  14. 14.

    Teschendorff AE, Yang Z, Wong A, Pipinikas CP, Jiao Y, Jones A, et al. Correlation of smoking-associated DNA methylation changes in buccal cells with DNA methylation changes in epithelial cancer. JAMA Oncol. 2015;1(4):476–85.

    PubMed  Article  PubMed Central  Google Scholar 

  15. 15.

    Zheng SC, Webster AP, Dong D, Feber A, Graham DG, Sullivan R, et al. A novel cell-type deconvolution algorithm reveals substantial contamination by immune cells in saliva, buccal and cervix. Epigenomics. 2018;10(7):925–40.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  16. 16.

    Kruppa J, Sieg M, Richter G, Pohrt A. Estimands in epigenome-wide association studies. Clin Epigenetics. 2021;13(1):98.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  17. 17.

    Teumer A, Holtfreter B, Volker U, Petersmann A, Nauck M, Biffar R, et al. Genome-wide association study of chronic periodontitis in a general German population. J Clin Periodontol. 2013;40(11):977–85.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  18. 18.

    Jaffe AE, Murakami P, Lee H, Leek JT, Fallin MD, Feinberg AP, et al. Bump hunting to identify differentially methylated regions in epigenetic epidemiology studies. Int J Epidemiol. 2012;41(1):200–9.

    PubMed  PubMed Central  Article  Google Scholar 

  19. 19.

    Phipson B, Maksimovic J, Oshlack A. missMethyl: an R package for analyzing data from Illumina’s HumanMethylation450 platform. Bioinformatics. 2016;32(2):286–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  20. 20.

    Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019;47(D1):D607–13.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  21. 21.

    Divaris K, Monda KL, North KE, Olshan AF, Lange EM, Moss K, et al. Genome-wide association study of periodontal pathogen colonization. J Dent Res. 2012;91(7 Suppl):21S-S28.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  22. 22.

    Bochenek G, Hasler R, El Mokhtari NE, Konig IR, Loos BG, Jepsen S, et al. The large non-coding RNA ANRIL, which is associated with atherosclerosis, periodontitis and several forms of cancer, regulates ADIPOR1, VAMP3 and C11ORF10. Hum Mol Genet. 2013;22(22):4516–27.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  23. 23.

    Munz M, Richter GM, Loos BG, Jepsen S, Divaris K, Offenbacher S, et al. Genome-wide association meta-analysis of coronary artery disease and periodontitis reveals a novel shared risk locus. Sci Rep. 2018;8(1):13678.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  24. 24.

    Consortium EP. An integrated encyclopedia of DNA elements in the human genome. Nature. 2012;489(7414):57–74.

    Article  CAS  Google Scholar 

  25. 25.

    Zimmerman MW, McQueeney KE, Isenberg JS, Pitt BR, Wasserloos KA, Homanics GE, et al. Protein-tyrosine phosphatase 4A3 (PTP4A3) promotes vascular endothelial growth factor signaling and enables endothelial cell motility. J Biol Chem. 2014;289(9):5904–13.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  26. 26.

    McQueeney KE, Salamoun JM, Ahn JG, Pekic P, Blanco IK, Struckman HL, et al. A chemical genetics approach identifies PTP4A3 as a regulator of colon cancer cell adhesion. FASEB J. 2018;32(10):5661–73.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  27. 27.

    Foy M, Anezo O, Saule S, Planque N. PRL-3/PTP4A3 phosphatase regulates integrin beta1 in adhesion structures during migration of human ocular melanoma cells. Exp Cell Res. 2017;353(2):88–99.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  28. 28.

    Maacha S, Anezo O, Foy M, Liot G, Mery L, Laurent C, et al. Protein tyrosine phosphatase 4A3 (PTP4A3) promotes human uveal melanoma aggressiveness through membrane accumulation of matrix metalloproteinase 14 (MMP14). Investig Ophthalmol Vis Sci. 2016;57(4):1982–90.

    CAS  Article  Google Scholar 

  29. 29.

    Mizuuchi E, Semba S, Kodama Y, Yokozaki H. Down-modulation of keratin 8 phosphorylation levels by PRL-3 contributes to colorectal carcinoma progression. Int J Cancer. 2009;124(8):1802–10.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  30. 30.

    Iida C, Ohsawa S, Taniguchi K, Yamamoto M, Morata G, Igaki T. JNK-mediated Slit-Robo signaling facilitates epithelial wound repair by extruding dying cells. Sci Rep. 2019;9(1):19549.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  31. 31.

    Shinkai A, Yoshisue H, Koike M, Shoji E, Nakagawa S, Saito A, et al. A novel human CC chemokine, eotaxin-3, which is expressed in IL-4-stimulated vascular endothelial cells, exhibits potent activity toward eosinophils. J Immunol. 1999;163(3):1602–10.

    CAS  PubMed  Google Scholar 

  32. 32.

    Kitaura M, Suzuki N, Imai T, Takagi S, Suzuki R, Nakajima T, et al. Molecular cloning of a novel human CC chemokine (Eotaxin-3) that is a functional ligand of CC chemokine receptor 3. J Biol Chem. 1999;274(39):27975–80.

    CAS  PubMed  Article  Google Scholar 

  33. 33.

    Blanchard C, Wang N, Stringer KF, Mishra A, Fulkerson PC, Abonia JP, et al. Eotaxin-3 and a uniquely conserved gene-expression profile in eosinophilic esophagitis. J Clin Investig. 2006;116(2):536–47.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  34. 34.

    Chevalier M, Rhee H, Elguindi EC, Blond SY. Interaction of murine BiP/GRP78 with the DnaJ homologue MTJ1. J Biol Chem. 2000;275(26):19620–7.

    CAS  PubMed  Article  Google Scholar 

  35. 35.

    Shields AM, Panayi GS, Corrigall VM. Resolution-associated molecular patterns (RAMP): RAMParts defending immunological homeostasis? Clin Exp Immunol. 2011;165(3):292–300.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  36. 36.

    Kroczynska B, Evangelista CM, Samant SS, Elguindi EC, Blond SY. The SANT2 domain of the murine tumor cell DnaJ-like protein 1 human homologue interacts with alpha1-antichymotrypsin and kinetically interferes with its serpin inhibitory activity. J Biol Chem. 2004;279(12):11432–43.

    CAS  PubMed  Article  Google Scholar 

  37. 37.

    Tanaka T, Moriwaki K, Murata S, Miyasaka M. LIM domain-containing adaptor, leupaxin, localizes in focal adhesion and suppresses the integrin-induced tyrosine phosphorylation of paxillin. Cancer Sci. 2010;101(2):363–8.

    CAS  PubMed  Article  Google Scholar 

  38. 38.

    Canny G, Levy O, Furuta GT, Narravula-Alipati S, Sisson RB, Serhan CN, et al. Lipid mediator-induced expression of bactericidal/ permeability-increasing protein (BPI) in human mucosal epithelia. Proc Natl Acad Sci U S A. 2002;99(6):3902–7.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  39. 39.

    Planello AC, Singhania R, Kron KJ, Bailey SD, Roulois D, Lupien M, et al. Pre-neoplastic epigenetic disruption of transcriptional enhancers in chronic inflammation. Oncotarget. 2016;7(13):15772–86.

    PubMed  PubMed Central  Article  Google Scholar 

  40. 40.

    Kim H, Momen-Heravi F, Chen S, Hoffmann P, Kebschull M, Papapanou PN. Differential DNA methylation and mRNA transcription in gingival tissues in periodontal health and disease. J Clin Periodontol. 2021;48(9):1152–64.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  41. 41.

    Carnevale G, Kaldahl WB. Osseous resective surgery. Periodontol. 2000;2000(22):59–87.

    Article  Google Scholar 

  42. 42.

    Lindberg K, Rheinwald JG. Three distinct keratinocyte subtypes identified in human oral epithelium by their patterns of keratin expression in culture and in xenografts. Differentiation. 1990;45(3):230–41.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  43. 43.

    Dickson MA, Hahn WC, Ino Y, Ronfard V, Wu JY, Weinberg RA, et al. Human keratinocytes that express hTERT and also bypass a p16(INK4a)-enforced mechanism that limits life span become immortal yet retain normal growth and differentiation characteristics. Mol Cell Biol. 2000;20(4):1436–47.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  44. 44.

    Aryee MJ, Jaffe AE, Corrada-Bravo H, Ladd-Acosta C, Feinberg AP, Hansen KD, et al. Minfi: a flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays. Bioinformatics. 2014;30(10):1363–9.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  45. 45.

    Smyth GK, Michaud J, Scott HS. Use of within-array replicate spots for assessing differential expression in microarray experiments. Bioinformatics. 2005;21(9):2067–75.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  46. 46.

    Pidsley R, Zotenko E, Peters TJ, Lawrence MG, Risbridger GP, Molloy P, et al. Critical evaluation of the Illumina MethylationEPIC BeadChip microarray for whole-genome DNA methylation profiling. Genome Biol. 2016;17(1):208.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  47. 47.

    McCartney DL, Walker RM, Morris SW, McIntosh AM, Porteous DJ, Evans KL. Identification of polymorphic and off-target probe binding sites on the Illumina Infinium MethylationEPIC BeadChip. Genom Data. 2016;9:22–4.

    PubMed  PubMed Central  Article  Google Scholar 

  48. 48.

    Du P, Zhang X, Huang CC, Jafari N, Kibbe WA, Hou L, et al. Comparison of Beta-value and M-value methods for quantifying methylation levels by microarray analysis. BMC Bioinformatics. 2010;11:587.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  49. 49.

    Fortin JP, Labbe A, Lemire M, Zanke BW, Hudson TJ, Fertig EJ, et al. Functional normalization of 450k methylation array data improves replication in large cancer studies. Genome Biol. 2014;15(12):503.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  50. 50.

    Bates D, Mächler M, Bolker B, Walker S. Fitting linear mixed-effects models using lme4. J Stat Softw. 2015;67(1):1–48.

    Article  Google Scholar 

  51. 51.

    Manimaran S, Selby HM, Okrah K, Ruberman C, Leek JT, Quackenbush J, et al. BatchQC: interactive software for evaluating sample and batch effects in genomic data. Bioinformatics. 2016;32(24):3836–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  52. 52.

    Benjamini Y, Hochberg Y. Controlling the false discovery rate—a practical and powerful approach to multiple testing. J R Stat Soc Ser B Methodol. 1995;57(1):289–300.

    Google Scholar 

  53. 53.

    The ARIC investigators. The Atherosclerosis Risk in Communities (ARIC) Study: design and objectives. Am J Epidemiol. 1989;129(4):687–702.

    Article  Google Scholar 

  54. 54.

    Zaimi I, Pei D, Koestler DC, Marsit CJ, De Vivo I, Tworoger SS, et al. Variation in DNA methylation of human blood over a 1-year period using the Illumina MethylationEPIC array. Epigenetics. 2018;13(10–11):1056–71.

    PubMed  PubMed Central  Article  Google Scholar 

  55. 55.

    Guo S, Diep D, Plongthongkum N, Fung HL, Zhang K, Zhang K. Identification of methylation haplotype blocks aids in deconvolution of heterogeneous tissue samples and tumor tissue-of-origin mapping from plasma DNA. Nat Genet. 2017;49(4):635–42.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  56. 56.

    Lister R, Pelizzola M, Dowen RH, Hawkins RD, Hon G, Tonti-Filippini J, et al. Human DNA methylomes at base resolution show widespread epigenomic differences. Nature. 2009;462(7271):315–22.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

Download references

Acknowledgements

We thank James Rheinwald (Boston, MA, USA) for providing us with the cell line OKG4/bmi1/TERT [42, 43]. Raw beta values from whole blood cells as a proxy for infiltrating immune cells were derived from [54], downloaded from the GEO database (GSE123914).

Funding

Open Access funding enabled and organized by Projekt DEAL. This study was funded by a research grant from the Deutsche Forschungsgemeinschaft (DFG), RI 2827/1-1, the Bundesministerium für Bildung und Forschung (01DL15002), and a grant of the DG PARO/CP GABA-Forschungsförderung.

Author information

Affiliations

Authors

Contributions

GR and AS conceived and designed the study, interpreted the data and obtained third-party funding. ETAD, CB, HD, CG, GH, OM, RN, NP, CR, YJS, and IS recruited the study participants and collected the biopsies. GR processed the biopsies in the laboratory. AF provided the infrastructure and supervised the laboratory work for the methylation arrays. JK performed the statistical analyses of the EWAS. GR analyzed the data. JK and GR created the figures for the manuscript. GR and AS wrote the manuscript and supervised the study. All authors helped draft the manuscript, participated in the manuscript editing and read/approved the final manuscript.

Corresponding author

Correspondence to Gesa M. Richter.

Ethics declarations

Ethics approval and consent to participate

This study was conducted in accordance with the principles of Good Clinical Practice and approved by the Independent Ethics Committee of the participating Universities. Written informed consent was obtained from all subjects according to the Declaration of Helsinki.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1

(xlsx). EpiDish results for cell culture and biopsies. Given are the estimations for the proportions of the major cell types epithelial cells (Epi), fibroblasts (Fib), and immune cells (IC) per sample.

Additional file 2

(png). Multidimensional Scaling (MDS) Plot for extraction sites. MDS Plot of the initial 120 samples, colored by site of extraction.

Additional file 3

(xlsx). Significant DMPs in full set. Results for all 15,507 DMPs significant in the full set of 48 samples of clinically uninflamed and 48 samples of inflamed gingiva.

Additional file 4

(xlsx). 441 DMPs from stringent filtering step. Results for the 441 DMPs in 193 clusters that were derived from a more stringent filtering process (q < 0.05, effect size ≥ 0.1, distance to nearest significant (q < 0.05) DMP \(\le\) 2 kb, padj < 10–5 in paired analysis (n = 39) for ≥ 1 DMP in the cluster. TPM = transcripts per million, mRNA-sequencing in healthy keratinized oral mucosa with 16mio reads/sample [12].

Additional file 5

(xlsx). Enrichment analysis of GO terms and KEGG pathways. Gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways for genes corresponding to the 2347 DMPs significant in the full and paired sub-analysis. BP = biological process, MF = molecular function, CC = cellular component.

Additional file 6

(png). Interaction networks of the 10 most significant genes from the GO term “synaptic vesicle cycle”. Interaction networks from the STRING database [20] for the top 10 differentially methylated genes in the GO term “synaptic vesicle cycle”.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Richter, G.M., Kruppa, J., Keceli, H.G. et al. Epigenetic adaptations of the masticatory mucosa to periodontal inflammation. Clin Epigenet 13, 203 (2021). https://doi.org/10.1186/s13148-021-01190-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13148-021-01190-7

Keywords

  • EWAS
  • Methylation
  • Periodontitis
  • Gingiva
  • Inflammation
  • Cell type deconvolution
  • ROBO2
  • PTP4A3