Gene Expression Profiles in Paired Gingival Biopsies from Periodontitis-Affected and Healthy Tissues Revealed by Massively Parallel Sequencing

Periodontitis is a chronic inflammatory disease affecting the soft tissue and bone that surrounds the teeth. Despite extensive research, distinctive genes responsible for the disease have not been identified. The objective of this study was to elucidate transcriptome changes in periodontitis, by investigating gene expression profiles in gingival tissue obtained from periodontitis-affected and healthy gingiva from the same patient, using RNA-sequencing. Gingival biopsies were obtained from a disease-affected and a healthy site from each of 10 individuals diagnosed with periodontitis. Enrichment analysis performed among uniquely expressed genes for the periodontitis-affected and healthy tissues revealed several regulated pathways indicative of inflammation for the periodontitis-affected condition. Hierarchical clustering of the sequenced biopsies demonstrated clustering according to the degree of inflammation, as observed histologically in the biopsies, rather than clustering at the individual level. Among the top 50 upregulated genes in periodontitis-affected tissues, we investigated two genes which have not previously been demonstrated to be involved in periodontitis. These included interferon regulatory factor 4 and chemokine (C-C motif) ligand 18, which were also expressed at the protein level in gingival biopsies from patients with periodontitis. In conclusion, this study provides a first step towards a quantitative comprehensive insight into the transcriptome changes in periodontitis. We demonstrate for the first time site-specific local variation in gene expression profiles of periodontitis-affected and healthy tissues obtained from patients with periodontitis, using RNA-seq. Further, we have identified novel genes expressed in periodontitis tissues, which may constitute potential therapeutic targets for future treatment strategies of periodontitis.


Introduction
Periodontitis is a chronic inflammatory disease characterized by the destruction of periodontal tissue. This common disease, primarily initiated by periodontal pathogens, is an outcome of a complex interaction between periodontal microorganisms and the host inflammatory response [1]. The host response involves proinflammatory cytokines, chemokines, prostaglandins, Toll-like receptors and proteolytic enzymes, which have all been demonstrated to play an important role in the pathogenesis of periodontitis [2,3].
Studies have been performed combining in vivo and in vitro approaches to identify genes responsible for periodontitis. To date, there are a few published microarray studies investigating the gene expression profile in periodontits. One microarray study reported no significant differences in gene expression at different pathological sites in patients with chronic and aggressive periodontitis [4], whereas Kim et al. [5] and Demmer et al. [6] showed a number of genes that were upregulated in periodontitis compared to healthy controls. In addition, Beikler et al. [7] demonstrated that in periodontitis sites, the expression of immune and inflammatory genes was down-regulated following non-surgical therapy. With regard to in vitro studies, gene expression profiling has been performed on gingival fibroblasts from inflamed and healthy gingival tissues, for a limited number of inflammatory markers, such as interleukin (IL)-1, IL-6, IL-8, tumor necrosis factor-a (TNF-a) and CD14 [8]. Furthermore, microarray analysis has also been performed on periodontal ligament cells and gingival keratinocytes [9,10]. With regard to disease susceptibility at a genomic level, one genome-wide association study (GWAS) has been conducted in patients with aggressive periodontitis showing an association between aggressive periodontitis and intronic single nucleotide polymorphism rs1537415, which is located in the glycosyltransferase gene GLT6D1 [11].
Despite research investigating periodontitis gene expression profiles through microarray analysis, specific genes responsible for the disease have not yet been found. However, the recent development of massively parallel sequencing has provided a more comprehensive and accurate tool for gene expression analysis through sequenced based assays of transcriptomes, RNA-Sequencing (RNA-Seq). This method enables analysis of the complexity of whole eukaryotic transcriptomes [12] and studies comparing RNA-Seq and microarrays have shown that RNA-Seq has less bias, a greater dynamic range, a lower frequency of false positive signals and higher reproducibility [13,14]. The aim of the present study was to investigate the general pattern of the gene expression profile in periodontitis using RNA-Seq. We also aimed to investigate the local variation in gene expression at site level, comparing periodontitis-affected and healthy gingival tissues obtained from the same patient.

Ethics Statement
The study was performed in accordance with the Declaration of Helsinki and the current legislation in Sweden and after approval from the Karolinska Institutet Ethical Research Board. The Regional Ethics Board in Stockholm approved the collection of the biopsies and informed consent was obtained from all patients.

Collection of gingival tissue samples
A total of 10 nonsmoking individuals (20 biopsies), were included in the study. Four patients in the study group had other types of diseases: patient 2 was undergoing investigations for the disease sarcoidosis, patient 3 had diabetes type-2, patient 7 had a history of osteoarthritis and patient 10 was diagnosed with asthma. All participants were examined for periodontal disease and those with a tooth site demonstrating a probing depth $6 mm, clinical attachment level $5 mm and bleeding on probing were included in the periodontitis-affected group, according to the clinical parameters previously used as indicators of periodontitis [15,16,17]. During flap surgery, two adjacent gingival biopsies with identical clinical status were harvested from a periodontal pocket affected by periodontitis. The sizes of the specimens were approximately 262 mm, and included the connective tissue and the epithelium. In the same subjects, two adjacent gingival biopsies with identical clinical status and of about the same size were also obtained from a clinically healthy gingival pocket. Clinically healthy pockets were defined as sites with no gingival/periodontal inflammation, no bleeding on probing, a probing depth #3.5 mm and a clinical attachment level #3.5 mm. One of the biopsies from each site was stored in RNA Later (Applied Biosystems, USA) overnight at 4uC and thereafter stored at 280uC for subsequent RNA isolation. The second biopsy from each site was used for histological and immunohistochemical analysis.

Hematoxylin-Eosin staining
Deparaffinized serial sections of gingival tissues were formalin fixed (4% neutral buffered formalin) and paraffin embedded. For assessment of orientation of the epithelium and connective tissue as well as the degree of inflammation, deparaffinized serial sections (4 mm) were prepared and sections of each biopsy were stained with Hematoxylin-Eosin (H&E). The degree of inflammatory cell infiltration was evaluated by three blinded observers, using a relative scale from 0 to 3, and statistical differences between periodontitis-affected and healthy sites were tested using the Wilcoxon signed-rank test.

Immunohistochemical stainings in gingival tissue
For staining of the T cell marker CD3, interferon regulatory factor 4 (IRF4) and chemokine (C-C motif) ligand 18 (CCL18), gingival tissues were rinsed in phosphate buffered saline (PBS) with 0.1% Saponin (PBS-Saponin buffer) for 10 min. After an antigen retrieval procedure, 10 mM Tris, 1 mM EDTA (pH 9.0) for CD3 and 0.01 M Citrate acid (pH 6.0) for IRF4 and CCL18, sections were blocked in 1% H 2 O 2 in PBS-Saponin for 60 min at room temperature (RT) for CD3 and for 45 min at RT for IRF4 and CCL18. Subsequently, tissues were rinsed in PBS-Saponin for 10 min and further treated with 3% bovine serum albumin (BSA) diluted in PBS-Saponin for 30 min at RT. The expression of CD3, IRF4 and CCL18 was investigated using CD3 polyclonal rabbit antibody (1 mg/ml, PBS-Saponin) from Dako Sweden AB (Stockholm, Sweden), IRF4 polyclonal rabbit antibody (0.5 mg/ml, PBS-Saponin) from Atlas antibodies (Stockholm, Sweden) and CCL18 polyclonal rabbit anti-human antibody (0.5 mg/ml, PBS-Saponin) from Sigma-Aldrich (St. Louis, MO, USA). Normal rabbit IgG from R&D systems (MN, USA) was used as negative control. After incubation with primary antibody, sections were blocked with 1% normal goat serum in PBS for 15 min. Afterwards, sections were incubated with a biotinylated secondary antibody provided in the Vectastain ABC-Elite Complex Kit (Vector labs, Burlingame, CA, USA) followed by application of the Elite ABC solution for 40 min at RT in the dark. Thereafter, sections were washed with PBS and the peroxidase activity was visualized with 0.3% (v/v) in DAB buffer containing 0.1% (v/v) H 2 O 2 . Finally, the slides were washed with distilled water, dehydrated through an ethanol series (70%, 95%, 99.9%) into xylene, mounted, and photographed using a light microscope. For CD3 stainings, the amount of positive cells was evaluated by three blinded observers, using a relative scale from 0 to 3, and statistical differences between periodontitisaffected and healthy biopsies were tested using the Wilcoxon signed-rank test.

RNA extraction
RNA was extracted from gingival biopsies using steel-bead matrix tubes and a tabletop Fast-Prep homogenizer by two sequential centrifugations for 20 s at speed 6.5 (Qbiogene, Irvie, CA, USA). The RNA was purified on RNeasy Spin Columns (Qiagen, Valencia, CA, USA), treated with DNAse H to ensure degradation of DNA, and thereafter eluated in RNase-free water. The average RNA yield was 15.6 mg. RNA quality was assessed using the RNA 6000 NanoLabChip Kit of the Bioanalyzer system from Agilent Technologies (Santa Clara, CA, USA).

Transcriptome sample preparation for sequencing
A total amount of 2-3 mg per sample was used as input material for the RNA sample preparations. All samples had RIN values above 8. The samples were bar-coded and prepared according to the protocol (Cat# RS-930-1001) from the manufacturer (Illumina, San Diego, CA, USA), as previously described by Stranneheim et al. [18]. All sample preparation reagents were taken from the Illumina mRNA Sample Preparation Kit or ordered from vendors specified in the mRNA sample preparation protocol, except for automation specific reagents: carboxylic acid beads used for precipitation; the ethanol and tetraethylene glycol (EtOH/TEG) and the Polyethylene Glycol and sodium chloride (PEG/NaCl) precipitation buffers.

Clustering and sequencing
The clustering of the bar-coded samples was performed on a cBot Cluster Generation System using an Illumina HiSeq Single Read Cluster Generation Kit according to the manufacturer's instructions. The library preparations were sequenced on an Illumina HiSeq 2000 as single-reads to 100 bp. Two sequencing runs were performed according to the manufacturer's instructions where two and three lanes were used in the first sequencing and second sequencing run, respectively (Table S1). The runs generated a total of 402 million reads with an average of 15 million reads per sample that passed the Illumina Chastity filter; these reads were included in the study.

Sequence analysis
All sequences were aligned to the human genome reference hg19 with TopHat [19,20] version 1.1.4 and Samtools [21] version 0.1.8 using TopHat standard parameters except for parameters -solexa1.3-quals -p 8 -GTF Homo_sapiens.GRCh37.59.gtf. Annotations from Ensembl and RefSeq, downloaded from University of California, Santa Cruz (UCSC) Genome Browser, were used to assign features to genomic positions. Sequences aligned to the human genome were assigned to features and counted using HTSeq version 0.4.6 with parameters -m intersection-strict -s no -t exon. The R/Bioconductor package DESeq [22] was used to call differential gene expression on read counts generated by HTSeq and to perform hierarchical clustering of samples. All biological replicates for healthy and periodontitis-affected had R 2 (Spearman) correlation of gene expression (read counts) above 0.92.

Functional analyses of gene lists using WebGestalt
Analyses of gene categories and pathways was performed using the WEB-based Gene Set Analysis Toolkit v2 (WebGestalt) [23] with parameters: Id Type: Ensembl_gene_stable_id, Ref Set: Entrez Gene, Significance Level: p,0.05, Statistics Test: Hypergeometric, MTC: BH, Minimum: 2. KEGG analysis was used for pathway enrichment analysis and the Gene ontology (GO)  category Biological process was used for the functional annotation analysis.

Patients and gingival tissues
A total of 10 patients, six males and four females, with a mean age of 5068, were included in the study. For each patient, a total of four gingival biopsies of about the same size were obtained from periodontitis-affected and healthy gingiva, with two biopsies from each site. Bleeding status, probing depth and degree of inflammation in the gingival tissues for each of the two gingival sites was recorded (Table 1). To assess the degree of gingival inflammation in the periodontitis-affected and healthy tissues, histological and immunohistochemistry staining was performed using H&E and anti-CD3 (Fig. 1). Scoring of the degree of inflammatory cell infiltration, assessed by H&E staining, and the amount of CD3 positive cells showed significantly higher inflammation in tissue from periodontitis-affected sites (p,0.01 for H&E and p,0.05 for CD3; Table 1).

RNA-Sequencing
We sequenced cDNA from 10 periodontitis-affected and 10 healthy gingival tissues, with an average of 15 million reads of 100 bp in length per sample. A pairwise approach, where each periodontitis-affected biopsy had a healthy counterpart from the same individual, was used to eliminate the background noise of individual-specific gene transcription, enabling acquisition of more relevant data from the cohort. Aligning the sequence reads against the human genome yielded a median of 68% of uniquely aligned reads across all samples. The expression pattern, based on RNA-Seq reads, of well-known inflammatory mediators IL-1b, IL-6, IL-8, TNFa, Regulated upon Activation, Normal T-cell Expressed, and Secreted (RANTES) and Monocyte Chemotactic Protein-1 (MCP-1) were analyzed in all the tissue samples. The expression (log 2 fold change) of these mediators was shown to be higher in the majority of the periodontitis-affected gingival tissue compared to healthy gingival tissue from the same patient (Fig. 2).  Distribution of gene transcripts between periodontitisaffected and healthy gingival tissues A total of 22 122 different mRNA transcripts were expressed in the periodontitis-affected and healthy gingival tissue samples. Among these transcripts, 1375 were unique to the periodontitisaffected tissue samples whereas 511 genes were uniquely transcribed in healthy gingival tissues (Fig. 3). KEGG enrichment analysis using WebGestalt [24] was performed among the unique genes for the periodontitis-affected and healthy tissues which revealed several regulated pathways indicative of inflammation for the periodontitis-affected condition ( Table 2 and Table S1). In contrast, in the healthy gingival tissues, regulated pathways indicated a non-inflammatory profile among the unique genes, as demonstrated in Table 3 and Table S1.

Clustering of biopsies
Unsupervised hierarchical clustering was performed on all gene transcripts having a median read count above a cutoff level set to 0.3 read counts per feature, to exclude expression due to spurious transcription (Fig. 4). The gingival tissues from periodontitis-affected sites from different patients showed a more similar gene expression pattern than healthy gingival tissues from the same patient. Clustering according to individual, where the paired healthy and periodontitis-affected biopsies cluster together, was only observed for patient 6 and 7. However, the biopsies showed a general trend of clustering according to the degree of inflammation as assessed by H&E staining (Table 1), except for sample 7H, sample 2H and an outlier sample 1H, which clustered separately. There was also a trend of forming larger clusters depending on sequence run, but paired biopsies (periodontits-affected and healthy) from each patient were always analyzed in the same sequence run.

Differential gene expression between periodontitisaffected and healthy gingival tissues
Differential gene expression between periodontitis-affected and healthy gingival tissues was analyzed using read counts for each gene with the DeSeq package [22]. The analysis revealed a total of 453 significantly (adj p,0.01) differentially expressed genes. Additional analyses of genes expressed in periodontitis-affected gingiva, showed that 381 genes were upregulated, whereas 72 genes were shown to be down-regulated (Fig. 5, Table S2).

Gene Ontology enrichment analysis of differentially expressed genes
Investigation of functional associations of gene expression changes in the tissue samples was performed using WebGestalt. Gene ontology (GO) Biological process was used for enrichment analysis. Significant gene enrichments (p,0.05) as well as their parent terms are demonstrated in Fig. 6. Several GO categories were over-represented among genes differentially expressed in periodontitis-affected versus healthy gingival tissues. The categories were mainly indicative of immune and inflammatory responses. Further enrichment analysis regarding Molecular function and Cellular components are provided in the supplementary data (Table S3).

Top 50 upregulated genes in periodontitis-affected gingival tissue
The top 50 significantly upregulated genes in periodontitisaffected gingival tissue with Unigene entry are displayed in Table 4 together with Ensemble ID, gene symbol, fold change, log 2 fold Figure 6. Gene ontology (GO) analysis of differentially expressed genes. All significant (p,0.05) Biological processes (GO categories) and their parent terms are shown. The color of each node illustrates the significance and can be interpreted using the scale bar, which displays the p value. Each node is also marked with the number of significantly regulated genes mapped to the GO category. doi:10.1371/journal.pone.0046440.g006  change and p value. We investigated whether there were any available reports on the involvement of these genes in periodontitis or other chronic inflammatory conditions. Among the top 50 upregulated genes, we identified a number of candidate genes, which were not previously demonstrated to be involved in periodontitis but have been shown to be associated with other chronic conditions such as rheumatoid arthritis (RA). These candidate genes included FCRL5, adenosine monophosphate deaminase 1 (AMPD1), CCL18, tumor-necrosis factor receptor superfamily 17 (TNFRSF17) and leukocyte immunoglobin-like receptor, subfamily A (without TM domain) member 3 (LILRA3), and IRF4 which has shown to be involved in chronic inflammatory diseases such as RA and inflammatory bowel disease (IBD), ( Table 5).

The protein expression of IRF4 and CCL18 in periodontitis-affected tissue
The expression of two of the top 50 differentially upregulated genes, IRF4 and CCL18 where further investigated at the protein level in gingival tissue samples from five additional patients with periodontitis. Immunohistochemical analysis showed that the transcription factor IRF4 and the chemokine CCL18 were expressed at the protein level in gingival tissue from patients with periodontitis (Fig. 7). IRF4 protein was expressed in cells including fibroblasts and inflammatory cells in the gingival connective tissue, as shown by morphology. For the chemokine CCL18, cellular staining of fibroblasts and inflammatory cells was observed, as well as some diffuse extracellular staining, consistent with chemokine secretion.

Discussion
This study provides a novel quantitative comprehensive mapping of gene expression in gingival tissues from patients diagnosed with periodontitis, using RNA-Seq.
We first confirmed that the degree of inflammation was higher in periodontitis-affected gingival tissue compared to healthy tissues obtained from the same individual. Our results were based on immunohistological staining of CD3 positive cells, and further verified by RNA-Seq quantification of gene expression of the established inflammatory markers IL-1b, IL-6, IL-8, TNFa, RANTES and MCP-1. These inflammatory mediators have previously been reported to be elevated in patients with periodontitis [25,26,27].
Next, we performed unsupervised clustering of the gingival tissues to get an overview of the data generated from the RNA-Seq analysis. Cluster analysis revealed that the majority of periodontitis-affected clustered together and the majority of the healthy gingival tissues also clustered together, which is in line with our results regarding inflammation in the tissues. The degree of inflammation, rather than the individual, seemed to affect the clustering, indicating a common gene expression profile for periodontitis. Our results, based on the gene expression pattern of the inflammatory markers (IL-1b, IL-6, IL-8, TNFa, RANTES and MCP-1) and the immunohistochemical evaluation, confirmed  that the inflammation in periodontitis involves elevated levels of locally produced cytokines in the periodontium, as has been previously demonstrated [28]. However, cluster analysis revealed that three of the patients (patient no. 6, 7 and 2) deviated from the clustering pattern. For example, the healthy gingival tissue collected from patient 6 clustered with the periodontitis-affected tissue, which could be due to moderate inflammatory infiltration (H&E score 2) observed in the healthy gingival tissue. The clustering pattern in tissue from patient 7, where the healthy and diseased gingival tissue also clustered together, could be partly explained by the patient's history of osteoarthritis, which is a disease associated with elevated levels of circulating proinflammatory cytokines IL-6 and TNFa [29]. The cluster pattern for patient 2 differed from the rest of the patient group, which could be related to this patient undergoing investigation for the inflammatory disease sarcoidosis, and in turn might affect the systemic inflammatory response. Previous studies report that oral manifestations of sarcoidosis include aggressive destruction of the periodontium with rapid periodontal bone loss [30,31,32]. One of these studies also emphasizes the importance of patients diagnosed with sarcoidosis to be evaluated for other systemic involvements [31]. Thus, regarding our clustering pattern, it cannot be ruled out that general health differences might have some effect on the final outcome. However, the comparison of the gene expression profiles of all individuals should minimize potential interfering signals originating from single individuals affected with other diseases. Our RNA-Seq analysis, investigating the gene expression profile in the gingival tissues showed that the genes were differentially distributed between healthy and periodontitis-affected samples. Enrichment analysis among uniquely expressed genes in the periodontitis-affected tissues showed regulated pathways indicative of inflammation, such as cytokine signaling, chemokine signaling and the JAK-STAT signaling pathway. Several cytokines such as interleukins, which are involved in periodontits, signal through the JAK-STAT signaling pathway [33]. On the other hand, in the healthy biopsies, pathways were indicative of non-inflammatory processes that may be involved in the maintenance of the healthy gingival tissue. Future studies should also include investigation of genes within these pathways, which may contribute to understanding, prevention and treatment of periodontitis.
Differential gene expression analyses of periodontitis-affected vs. healthy gingival tissues showed the majority of differentially expressed genes to be upregulated in the periodontitis-affected tissues. Furthermore, GO enrichment analysis among these differentially expressed genes demonstrated that most of these genes were involved in immune and inflammatory processes. This is in line with the increased inflammatory response in the tissue, and also in accordance with our previous microarray studies on inflammatory-stimulated cell cultures reporting that gene expression profiles of TNFa-stimulated cells show an induction of inflammatory genes [34,35].
Up to date, RNA-Seq studies aimed to identify new genes involved in the pathogenesis of periodontits have not been reported. One ab initio study by Covani et al. [36] identified genes with potential roles in periodontitis, some of which have not previously been associated with the disease. However, the protein expression of these genes in periodontitis-affected tissues has not been confirmed. In our study we aimed to identify genes involved in the pathogenesis of periodontitis. Therefore, we further searched through the differentially expressed genes, focusing on the top 50 upregulated genes. Two of these 50 upregulated genes, IRF4 and CCL18, were also detected at the protein level in periodontitis affected-tissues, supporting these genes as novel finds in the pathogenesis of periodontitis. Furthermore, these two selected genes have been reported to be involved in other chronic inflammatory diseases such as RA. The transcription factor, IRF4, has been demonstrated to be involved in T-cell-dependent chronic inflammatory diseases such as IBD [37]. Mudter et al. 2011 reported a correlation between mRNA levels of IRF4 and production of cytokines such as IL-6 and IL-17 in the inflamed colon from patients with IBD, indicating that IRF4 is involved in the regulation of chronic mucosal inflammation [37]. In addition, the gene for CCL18 was upregulated in periodontitis-affected gingival tissues. This chemokine, expressed by macrophages, monocytes and dendritic cells, has been demonstrated to be increased in synovial tissue of RA patients [38]. It has also been suggested that blockage of CCL18 expression by anti-TNF-a antibodies identifies CCL18 as an additional target for anti-TNF-a therapy in patients with RA [39,40]. Studies are currently ongoing to investigate the expression of candidate genes novel for periodontitis in a larger cohort of patients with periodontitis and healthy controls, to be able to evaluate their impact and to further explore the possible therapeutic targeting of these genes. In addition, future studies will also be performed investigating the biological significance of the down-regulated genes in periodontitis.
In conclusion, we demonstrate for the first time, using RNA-seq, profile analysis of periodontitis revealing site-specific local variation in gene expression profiles of periodontitis-affected and healthy tissues obtained from patients diagnosed with periodontitis. Furthermore, we have identified differentially expressed novel genes in gingival tissue of periodontitis. Our findings provide a first step towards a quantitative comprehensive insight into the transcriptome of gingival tissue from patients with periodontitis, to enable identification of possible diagnostic markers of periodontitis as well as potential therapeutic targets.

Supporting Information
Table S1 Uniquely expressed genes within enriched pathways in periodontitis-affected and healthy gingival tissues. (XLSX)