PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Gene. Author manuscript; available in PMC 2017 May 1.
Published in final edited form as:
PMCID: PMC4767619
NIHMSID: NIHMS754434

RNA-sequencing study of peripheral blood monocytes in chronic periodontitis

Abstract

Background

Monocytes are an important cell type in chronic periodontitis (CP) by interacting with oral bacteria and mediating host immune response. The aim of this study was to reveal new functional genes and pathways for CP at monocyte transcriptomic level.

Methods

We performed an RNA-sequencing (RNA-seq) study of peripheral blood monocytes (PBMs) in 5 non-smoking moderate to severe CP (case) individuals vs. 5 controls. We took advantage of a microarray study of periodontitis to support our findings. We also performed pathway-based analysis on the identified differentially expressed (DEx) transcripts/isoforms using DAVID (Database for Annotation, Visualization and Integrated Discovery).

Results

Through differential expression analyses at both whole gene (or whole non-coding RNA) and isoform levels, we identified 380 DEx transcripts and 5,955 DEx isoforms with a PPEE (posterior probability of equal expression) of < 0.05. Pervasive up-regulation of transcripts at isoform level in CP vs. control individuals were observed, suggesting a more functionally active monocyte transcriptome for CP. By comparing with the microarray dataset, we identified several CP-associated novel genes (e.g., FACR and CUX1) that have functions to interact with invading microorganisms or enhance TNF production on lipopolysaccharide stimulation. DAVID analysis of both the RNA-seq and the microarray datasets leads to converging evidence supporting “endocytosis”, “cytokine production” and “apoptosis” as significant biological processes in CP.

Conclusions

As the first RNA-seq study of PBMs for CP, this study provided novel findings at both gene (e.g., FCAR and CUX1) and biological process level. The findings will contribute to better understanding of CP disease mechanisms.

Keywords: periodontitis, RNA sequencing analysis, monocytes, apoptosis, endocytosis, cytokine, FCAR, CUX1

Background

Chronic periodontitis (CP) [1] is an inflammatory disease that affects the supporting structures of teeth. It is one of the most prevalent chronic diseases and a major cause of loss of permanent teeth [2], which severely affects the quality of life in old people.

While oral bacteria are necessary for initiating pathological cascades of CP [3], it is host immune response and inflammatory process that is responsible for periodontal tissue damage [4], where monocytes play a central role. Monocytes, by differentiating into dendritic cells (DCs) in periodontium [5], trigger T cell differentiation into type 1, type 2 or type 17 T-helper (Th) cells, which in turn direct cellular, humoral and innate immune responses, respectively [6, 7]. Monocytes, by interacting with oral bacteria, initiate inflammatory processes in periodontium and produce a series of inflammatory factors, such as PGE2, IL-1-β, TNF-α and MMP [8].

Although the important pathological cascades mediated by monocytes mainly occur in periodontium, peripheral blood monocytes (PBMs), as functional continuum and the sole source of the monocytes in the periodontium, play a critical role in CP pathogenesis. It was found that in CP patients, PBMs had a “specific functional profile” that favors the Th2-cell response over Th1-cell response [912]. Therefore, functional aberrance of PBMs may predispose an individual to a destructive host immune response leading to CP.

While the immune response and inflammatory processes are well known as the key mechanisms for CP, the detailed functional genes and pathways underlying the disease at the PBMs transcriptome level are still unclear. Despite two recent genome-wide gene expression studies of CP [13, 14] the studies were performed on the gingival tissue that may have a different transcriptomic landscape from PBMs, and hence the findings from those studies may not be relevant to functional genes contributing to CP at the PBMs level. Another study [15] was conducted on PBMs on CP patient using the microarray technology. However, the interesting findings from the study may need further confirmation from an independent study, such as the one presented here.

Here, to shed light on those important genes and pathways for CP, we used the latest next generation sequencing technology and performed the first RNA-sequencing (RNA-seq) analysis of PBMs for CP. Our findings contributed to further understanding of CP pathogenesis by having identified a number of novel yet functionally relevant genes and biological processes as specific key players for the disease.

Methods

Study participants characteristics

The study was approved by Institutional Review Board of Louisiana State University Health Science Center (LSUHSC). Signed informed-consent documents were obtained from all study participants before entering the study.

The study participants (cases and controls) were recruited at LSUHSC School of Dentistry (LSUHSC SOD) periodontics clinic. All case individuals (n=5) were non-smoking African American (AA) males diagnosed with generalized moderate to severe CP as per the 1999 world workshop classification [16].

Specifically, cases with >30% sites with clinical attachment loss (CAL) of ≥3–4 mm were diagnosed as generalized “moderate” chronic periodontitis and cases with >30% sites with CAL ≥5 mm were diagnosed as generalized “severe” chronic periodontitis.

Control individuals (n=5) were periodontally healthy based on clinical and radiographic examinations. None of the controls exhibited CAL ≥ 3 mm, probing depths ≥ 4 mm or radiographic evidence of bone loss. All control individuals were also non-smokers and gender and ethnically matched to the cases..

The average age of case individuals was 44.6 years with a standard deviation of 7.73 years. The average age of control individuals was 36.4 years with a standard deviation of 12.97 years. There was no statistical difference in age between case and control individuals (p = 0.26) according to the t test.

Exclusion criteria for all participants were as follows: 1) Recent treatment for periodontal disease (within a year) 2) Current and/or long-term use of NSAIDs (non-steroidal anti-inflammatory drugs) which could affect the outcome of the study; 3) Former or current smokers, since smoking is one of the most influential factors in periodontitis; 4) Diabetes mellitus; 5) Treatment with corticosteroid therapy currently or for more than 6 months duration at any time; 6) Treatment with anticonvulsant therapy or anti-depressants currently or for more than 6 months duration at any time; 7) Treatment with antibiotics currently or for more than 6 months duration at any time; 8) Evidence of other metabolic or inherited bone diseases such as hyper or hypoparathyroidism, Paget’s disease, osteomalacia, osteognesis imperfect or others; 9) Hyperthyroidism; 10) History of or current diagnosis with aggressive periodontitis.

Monocyte isolation

Monocyte isolation from whole blood was performed with a monocyte-negative isolation kit (Miltenyi Biotec Inc, Auburn, CA) following manufacturer’s recommendation. The kit is an indirect magnetic labeling system for isolating untouched monocytes from human peripheral blood mononuclear cells (PBMCs). Non-monocytes, i.e. T cells, NK cells, B cells, dendritic cells and basophils, are indirectly magnetically labeled using a cocktail of biotin-conjugated antibodies against CD3, CD7, CD16, CD19, CD56, CD123 and Glycophorin A, and Anti-Biotin MicroBeads. Isolation of highly pure unlabeled monocytes is achieved by depletion of the magnetically labeled cells. Based on our extensive experience of monocyte isolation using this kit [17, 18], the isolated monocytes can achieve a purity of at least 85%.

Total RNA extraction

Total RNA from monocytes was extracted using Qiagen RNeasy Mini kit (Qiagen, Inc., Valencia, CA). We used Agilent Bioanalyzer (Agilent, Santa Clara, CA) to control the RNA quality before RNA-seq experiment, where RNA integrity number (RIN) was confirmed at no less than 7.0 before RNA-seq experiment.

RNA-seq

RNA sequencing was performed at JP Sulzberger Columbia Genome Center. We used poly-A pull-down to enrich mRNAs from total RNA samples (200ng–1ug per sample) and proceeded on library preparation by using Illumina TruSeq RNA prep kit (Illumina Inc., San Diego, CA).

Specifically, we purified the polyA containing mRNA molecules using oligo-dT attached magnetic beads with two rounds of purification. During the second elution of the polyA RNA, the RNA was fragmented and primed for cDNA synthesis. We then performed reverse transcription on the RNA fragments to synthesize first strand cDNA using reverse transcriptase and random primers. The RNA template was then removed and a 2nd strand cDNA was synthesized. End repair was then performed to convert the overhangs into blunt ends using the End Repair Mix. After 3′ ends adenylation, multiple indexing adapters were ligated to the ends of the ds cDNA and those DNA fragments with adapter molecules were enriched with PCR reactions that constructed a library for downstream sequencing experiments. The library was validated using Agilent Bioanalyzer (Agilent, Santa Clara, CA) to ensure the size and purity of the sample by observing a band at approximately 260 bp.

Libraries were then sequenced using an Illumina HiSeq 2500 instrument (Illumina Inc., San Diego, CA). We multiplexed samples in each lane, which yielded targeted number of single-end 100bp reads for each sample, as a fraction of 180 million reads for the whole lane.

We used the RTA software (Illumina Inc., San Diego, CA) for base calling and CASAVA software (version 1.8.2) (Illumina Inc., San Diego, CA) for converting the BCL files to fastq format, coupled with adaptor trimming. The read depth for each sample was 30 million reads.

The raw RNA-seq data were deposited to GEO with an accession number of GSE61490.

Differential expression analysis of RNA-seq data

We used RSEM v1.2.11 [19] to analyze the RNA-seq data for alignment and expression calculation. Under RSEM, we used EBSeq [20] to identify differentially expressed (DEx) transcripts at whole gene (or whole non-coding RNA) level or at isoform level.

EBSeq assumes the expected count (X) for a gene or an isoform is distributed as Negative Binomial (NB), i.e., X~NB(r; p), and determines the posterior probability of differential expression of a gene or an isoform through the Bayes’ rule. To achieve that, EBSeq models a prior distribution of the data using β distribution. Under the null hypothesis (EE: equal expression), the data XgiC1,C2=XgiC1,XgiC2 (where gi denote a gene or isoform i, C1 condition 1 and C2 condition 2) can be represented by a prior predictive distribution f0Ig(XgiC1,C2) (which is constructed based on beta functions). Under the alternative (DE: differential expression), XgiC1,C2 follows another prior predictive distribution f1Ig(XgiC1,C2). Given a latent variable Zgi, where Zgi = 1 represents that the gene/isoform gi is DE and Zgi = 0 represents that gi is EE, and also Zgi~Bernoulli(p). Based on that, the marginal distribution of XgiC1,C2 and Zgi is (1-p)f0Ig(XgiC1,C2)+pf1Ig(XgiC1,C2) Therefore, the posterior probability of DE for gi is defined by Bayes’ rule:

pf1Ig(XgiC1,C2)(1-p)f0Ig(XgiC1,C2)+pf1Ig(XgiC1,C2)

For more details, please refer to the supplement to the paper on EBSeq [20].

An online microarray dataset to provide supportive evidence

We took advantage of an online microarray dataset from GEO (GEO accession #: GSE6751) [15] as an independent dataset to support the findings from our RNA-seq data. The dataset contains 15 CP patients, whose PBMs expression was analyzed 1 week before periodontal treatment (time point #1), at treatment initiation (time point #2), 6 weeks after the treatment (time point #3) and 10 weeks after the treatment (time point #4) [15]. At 10 weeks after treatment, considerable improvement in periodontal status, including bleeding sites per individual (decreased from an average of 81% before treatment to an average of 25% after treatment) and deep pockets per individual (decreased from an average of 69.2 before treatment to an average of 9.2 after treatment), was achieved [15]. The treatment involved only periodontal surgery and tooth extractions, if necessary, without local or systemic antibiotics.

We took advantage of the 15 individuals’ PBMs samples before treatment and the samples 10 weeks after treatment as a “paired” cohort with difference in periodontal status. Given the considerable improvement of periodontal status 10 weeks after treatment vs. before treatment, we treat the individuals before treatment as “periodontally unhealthy status” and the individuals 10 weeks after treatment as “periodontally healthy status”. Comparing the two statuses in gene expression may shed lights on those genes differentially expressed in terms of “difference in periodontal health”. Such differentially expressed genes may provide supporting evidence for those differentially expressed genes as identified in our RNA-seq data, where two groups of individuals are also different in terms of periodontal health.

The PBMs gene expression was assayed using Affymetrix U133 Plus 2.0 arrays. We analyzed the raw cel files for each individual before the periodontal treatment (time point #1) and 10 weeks after the treatment (time point #4). The cel files were normalized using RMA [21]. We then performed differential expression analysis using paired t test through the Bioconductor’s LIMMA (linear models for microarray data) package [22] to compare the PBMs gene expression before vs. after periodontal treatment. Genes and their differential expression p values as assessed above were used as the basis for supporting those DEx transcripts identified in our RNA-seq data.

Enrichment analysis using DAVID

GO (gene ontology) and other pathway/functional term enrichment analysis was performed on a number of group of genes (e.g., DEx transcripts or isoforms identified in CP vs. control individuals) using DAVID (Database for Annotation, Visualization and Integrated Discovery) v6.7 software package (http://david.abcc.ncifcrf.gov/) [23].

Results

Quality control (QC) analysis

The mean number of reads over the 10 samples is 26,549,771, with a range from 21,697,093 reads to 31,754,419 reads. The average mapping rate (ratio of number of uniquely mapped reads over total number of reads) across the 10 samples is 75.93%, with a range from 69% to 84%.

The sequence QC results on base quality are quite excellent. According to the FastQC software (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/), all of the 10 samples pass the QC (with a green light tick) for both the “Per base sequence quality” and “Per sequence quality scores”, two important metrics for base quality QC. Specifically, the average quality score per read is either 37 or 38 for all the samples, which is very good as the maximum quality score is 40.

Our sequence read depth turns out to be sufficient for quantifying lowly expressed transcripts as we have detected a large number of transcripts with a low count. As shown in our data, we have observed >4,400 transcripts with a count of 10–50 for at least one sample. (To preclude false expression signals, these transcripts do not include those whose counts are either 0 or lower than 10 in any sample.) Among these, there are >600 transcripts with an average count of <50 across the 10 samples, representing transcripts with an overall low expression for all the 10 samples. Note that although the level of expression is very low these >600 transcripts are consistently detected in all the samples. Hence it is unlikely that these detected transcripts are false signals (e.g., due to mistaken alignment).

RNA-seq analysis at the whole gene level

At the whole gene (or whole non-coding RNA) level, we identified a total of 380 DEx transcripts in CP vs. control individuals, which achieved a PPEE (posterior probability of equal expression) of < 0.05 and a PPDE (posterior probability of differential expression) of > 0.95. Among these transcripts, 228 were up-regulated and 152 were down-regulated in CP vs. control individuals.

As a large number of DEx transcripts achieved a PPEE of 0 or near 0 (or PPDE of 1 or near 1), to rank the DEx transcripts, we used the measure of PostFC (posterior fold change). Listed in Table 1 are the top 10 up-regulated transcripts with the largest PostFC in CP vs. control individuals and the top 10 down-regulated transcripts with the smallest PostFC in CP vs. controls.

Table 1
The DEx transcripts identified at the whole gene level with the most extreme fold changes

Among the transcripts with the most extreme PostFC in CP vs. control individuals (Table 1), representing those transcripts that may be strongly regulated in terms of CP status, LST1 was previously found to be up-regulated in response to LPS and bacterial infection and in blood of rheumatoid arthritis (RA) patients as compared to controls [24]. This gene’s known function is consistent with the context of CP pathogenesis, where LPS and bacterial infection initiates the pathologic cascade of the disease. Another up-regulated gene, LILRA5 was found to be expressed by synovial tissue macrophages in RA patients and induce pro-inflammatory cytokines [25]. RNF5, down-regulated in CP patients vs. controls, was found to influence susceptibility to bacterial infection [26]. Another down-regulated transcript, IER3, was reported in a study, where deficiency of IER3 expression in mice may result in aberrant immune regulation and enhanced inflammation [27].

As another interesting finding, almost half of the transcripts (n=8) presented in Table 1 are located in the MHC (major histocompatibility complex) regions, which include 2 up-regulated transcripts, LST1 [28] and LY6G5C [29], and 6 down-regulated transcripts, HLA-DRA [30], ZNRD1 [31], GABBR1 [32], RNF5 [33], TRIM26 [33], and PPP1R11 [34]. The co-localization to the MHC regions of these strongly regulated transcripts in terms of CP status once again suggests immune aberration as an important component of CP pathogenesis.

Pathway-based analysis using DAVID for those transcripts identified at the whole gene level

Using DAVID, we analyzed the 380 DEx transcripts that achieved a PPEE of <0.05. The top ten enriched pathways/functional terms are listed in Table 2. These terms achieved significant enrichment p values even after adjusted with the Bonferroni method.

Table 2
Top 10 enriched functional terms identified by DAVID for the DEx transcripts identified at the whole gene level

Among these ten terms, two are related to immune response (GO:0006955 and SP_PIR_keyword immune response), three are related to immunoglobulin/major histocompatibility complex (MHC) (IPR003006, IPR003597 and GO:0042611) and one is related to antigen processing and presentation (GO:0019822).

Differential expression analysis of the microarray dataset

Using Bioconductor’s LIMMA package [35], we identified 2,354 probe sets, corresponding to 2,092 annotated genes, which were DEx in periodontally unhealthy (one week before treatment) vs. periodontally healthy statuses (10 weeks after periodontal treatment), among which, 1,310 genes were found to be up-regulated and 782 genes down-regulated before vs. after treatment at a p value of <0.05.

Supporting evidence from the microarray data on those transcripts identified at the whole gene level

At the whole gene level, we identified in the RNA-seq data 380 DEx transcripts, which achieved a PPEE (posterior probability of equal expression) of < 0.05 and a PPDE (posterior probability of differential expression) of > 0.95. We selected a subset with a higher significance, which were 210 DEx transcripts that achieved a PPDE =1 (and a PPEE < 4.70E-10), for finding supporting evidence from the microarray dataset. Those transcripts in the list of 210 DEx transcripts were checked for differential expression signals in the microarray data. Through this analysis, we identified 11 genes that achieved a PPEE of 1.78E-15 or lower in the RNA-seq dataset and a p value of <0.05 in the microarray dataset. The 11 genes also had the same direction of regulation across the two datasets. The direction of regulation in the RNA-seq dataset is defined based on the ratio of posterior mean expression in CP patients over controls, with a ratio of >1 as up-regulation and a ratio of <1 as down-regulation. The direction of regulation in the microarray dataset is defined based on the ratio of mean expression in the periodontally unhealthy status (before treatment) over periodontally healthy status (after treatment), with a ratio of >1 as up-regulation and a ratio of <1 as down-regulation. Table 3 shows the information of the 11 genes supported by the microarray dataset.

Table 3
The DEx transcripts identified at the whole gene level and supported by the PBMs microarray dataset

Among the 11 DEx transcripts with supporting evidence from the microarray data (Table 3), FCAR achieved the most significant p value (7.78E-4) in the microarray study [15]. In addition, multiple probe sets (207674_at, 211307_s_at, 211816_x_at) from this gene were associated with periodontal health status in the microarray study [15]. According to OMIM description (entry #: 147045), FCAR “interacts with aggregated IgAs, such as IgA coated on the surface of an invading microorganism, and mediates several immunologic defense processes such as phagocytosis, antibody-dependent cell-mediated cytotoxicity, and stimulation of the release of inflammatory mediators” [36]. This gene’s expression in monocyte or mononuclear cells has been associated with several autoimmune diseases, such as autoimmune thyroid disease [37], multiple sclerosis [38] and ankylosing spondylitis [39]. Given FCAR’s roles in immunologic defense of microorganism and the gene’s importance in autoimmune diseases, FCAR may also play a key role in CP pathogenesis.

In addition to FCAR, other 6 DEx genes supported by the microarray dataset (Table 3) may also be involved in CP pathogenesis as they have a close link with immune response or immune-related diseases. CUX1 was found to enhance TNF production on LPS stimulation [40]. RNASE3 is a mediator in host immune response to parasites, bacteria and viruses and may also cause side-effects on the host’s own tissues [41]. CLCN5 was found to play a role in immunopathogenesis of ulcerative colitis [42]. REL is a transcription factor that regulates innate and adoptive immunity. In particular, REL is essential for IL-12 and IL-23 transcription by macrophages and dendritic cells, which is crucial for T-cell differentiation and effector functions [43]. VNN2’s up-regulation in mononuclear cells was associated with both inflammatory bowel disease and RA [44]. HLA-DOA, as a MHC class II gene, has been associated with asthma [45], RA [46] and type I diabetes [47].

Supporting evidence at the pathway level from the microarray data for those transcripts identified at the whole gene level

Using DAVID, we analyzed the 380 DEx transcripts (PPEE < 0.05) identified at the whole gene level from the RNA-seq data. We did the same analysis using DAVID on those 2,092 DEx genes (p < 0.05) as detected in the microarray dataset. The identified pathways and functional terms that achieved a Bonferroni enrichment p value of < 0.05 in the RNA-seq data were checked for enrichment signals in the DAVID analytical results of the microarray data. Through this analysis, we identified one functional term “SP_PIR_Keyword: immune response”, which achieved an enrichment p value of 7.27E-8 (Bonferroni p value = 2.58E-5) in the RNA-seq data and an enrichment p value of 2.24E-2 in the microarray data.

In the RNA-seq data, the component genes that fall into this term include HLA-DRA, HLA-H, HLA-DMB, HLA-DPB1, MICB, HLA-DOB, HLA-DOA, C1QB, HLA-F, TMED7, LILRB4, TLR9, TAP2, HLA-DQA2, LST1 and LILRA5. For the microarray data, the component genes that fall into this term include CD7, IL23A, LIME1, DCLRE1C, POLR3G, TLR3, POLR3B, IL2, HLA-DOA, NLRX1, POLR3E, IFIH1, TLR10, HFE, CD300LB, INPPL1, CNPY3, TLR7, MYD88, TLR6, TLR5, LILRA1, TLR4, CR1, HLA-DRB1, CLU, CLEC7A, TLR8, CLEC5A, LY86, CLEC4E, CLEC4A and TLR2.

Differential expression analysis at the isoform level using the RNA-seq data

We identified 5,955 DEx at the isoform level, which achieved a PPEE of <0.05 and a PPDE of >0.95. Due to the large number of DEx isoforms identified, to control for false positives, we focused on a subset of isoforms that achieved a higher significance for differential expression, which are 2,344 DEx transcripts that achieved a PPDE of 1 (and a PPEE from 0 to 4.89E-10).

The majority of the 2,344 isoforms (~85%) were up-regulated in CP vs. control individuals. These isoforms are from 1,816 distinct genes or non-coding RNAs. We used DAVID to analyze these isoforms and the top 10 most enriched functional terms are listed in Table 4.

Table 4
Top 10 enriched functional terms identified by DAVID for the DEx isoforms

Supporting evidence from the microarray data on the DEx isoforms identified from the RNA-seq data

Due to the large number (n = 2,344) of DEx isoforms that achieved a PPDE of 1, we focused on the subset of transcripts (with a higher significance) that achieved both a PPDE of 1 and PPEE of 0 (1,843 isoforms) to find supporting evidence from the microarray data. For these isoforms, we checked for the DEx signals at the gene level from the microarray data. The analysis identified supporting evidence for 213 isoforms from 177 distinct genes or non-coding RNAs, which achieved a p value of less than 0.05 in the microarray data. Again they all had the same direction of regulation in terms of CP status across the two datasets. Almost all of the 213 isoforms/genes were up-regulated in CP cases vs. controls (or periodontally unhealthy status vs periodontally healthy status); among these isoforms, only those from 7 genes were down-regulated in CP cases vs. controls, which are URGCP, RPS20, FAM98A, XBP1, G3BP1, NFAT5 and ZNF207.

For these 213 isoforms/genes supported by the microarray dataset, we used DAVID to map them to functional terms and pathways. The top 10 most enriched functional terms/pathways were all GO terms (listed in Table 5). Almost all of these GO terms were significant even after Bonferroni correction (Bonferroni enrichment p values from 5.83E-3 to 5.06E-2). Among these terms are endocytosis (GO:0006897), cytokine production (GO:0001816) and apoptosis (GO:0006915).

Table 5
Top 10 enriched functional terms identified by DAVID for the DEx isoforms supported by the microarray dataset

Thirteen genes (DEx in both RNA-seq and microarray datasets) belong to the GO term “GO:0006897~endocytosis”, which are DENND1A, RUFY1, CORO1C, ASGR2, APP, DAB2, PICALM, CD36, AP1S2, CLEC7A, THBS1, CLCN5 and RIN3. All of these genes/isoforms were up-regulated in CP patients vs. controls (or periodontally unhealthy vs. periodontally healthy statuses).

Seven genes (DEx in both RNA-seq and microarray datasets) belong to the GO term “GO:0001816~cytokine production”, which are NLRC4, G6PD, MYD88, NFAT5, TLR4, NLRP3 and PTAFR. Except for NFAT5, all of the genes were up-regulated in CP patients vs. controls (or periodontally unhealthy vs. periodontally healthy statuses).

Twenty genes (DEx in both RNA-seq and microarray datasets) belong to the GO term “GO:0006915~apoptosis”, which are ARHGEF2, SGK1, DNM1L, XIAP, UBE4B, CIDEB, STK17B, TRIO, NLRP3, BCL2L13, NCSTN, TNFRSF1A, PEA15, NLRC4, APP, GSN, HIPK3, BNIP3L, NLRP12, THBS1. All of these genes were up-regulated in CP patients vs. controls (or periodontally unhealthy vs. periodontally healthy statuses). Among these genes, ten are categorized by DAVID as apoptosis induction genes (BCL2L13, BNIP3L, NLRC4, NLRP3, ARHGEF2, ARP, CIDEB, NCSTN, STK17B, TRIO), whereas only 5 were categorized as anti-apoptosis genes (BNIP3L, XIAP, HIPK, PEA15 and THBS1).

Supporting evidence from the microarray data at the pathway level for the DEx isoforms identified in the RNA-seq data

We used DAVID to map those 2,344 DEx isoforms, which achieved a PPDE of 1 (and a PPEE from 0 to 4.89E-10) in the RNA-seq dataset. For those functional terms/pathways that achieved a Bonferroni-corrected enrichment p value of < 0.05, we cross-checked the enrichment signals against the DAVID analytical results from the microarray data.

Through this analysis, we identified 7 functional terms that achieved a Bonferroni-corrected enrichment p value of <0.05 in both the RNA-seq data and the microarray data. These functional terms include 6 GO terms and 1 SP_PIR_Keywords term (listed in Table 6). Among these terms, 3 are related to apoptosis: GO:0042981~regulation of apoptosis, GO:0043067~regulation of programmed cell death and GO:0010941~regulation of cell death, and other 4 are GO:0031982~vesicle, GO:0031410~cytoplasmic vesicle, GO:0010627~regulation of protein kinase cascade and actin-binding (SP_PIR_KEYWORDS).

Table 6
The functional terms identified by DAVID on the DEx isoforms and supported by the microarray dataset

In the RNA-seq data, 121 genes belong to the GO apoptosis-related terms. In the microarray data, 137 genes belong to the GO apoptosis-related terms. By comparing the 121 vs. 137 component genes across the two datasets, 22 genes (29 isoforms from RNA-seq data vs. 27 probe sets from microarray data) were shared between the two datasets (Table 7). These isoforms/probe sets were all up-regulated in CP vs. control individuals (or periodontally healthy vs. unhealthy statuses). According to DAVID analysis, among the 22 genes, 15 are classified as those for apoptosis induction (TRIO, ARHGEF2, BCL6, TLR4, STK17B, NLRP3, BCL2L13, IKBKG, BNIP3L, APP, NLRC4, RXRA, CIDEB, NCSTN, BID), whereas only 8 as those for anti-apoptosis (HIPK3, BCL6, BNIP3L, PEA15, XIAP, NUP62, THBS1, MYD88) (Table 7). (Note that BCL6 and BNIP3L are classified as both apoptosis induction and anti-apoptosis genes.)

Table 7
The component genes of the apoptosis GO term shared by the RNA-seq and the microarray datasets

Discussion

We performed the first RNA-seq study of PBMs for CP. Although we had a small sample size, the DEx transcripts as detected by the RNA-seq analysis were biologically relevant; the top 10 enriched functional terms identified by DAVID [23] for the DEx transcripts identified at the whole gene level (Table 2) were almost all related to immune response, antigen processing and presentation and immunoglobulin/MHC. As host immune response is well recognized as the key basis for CP pathogenesis, this result (Table 2) lent strong support to the biological validity of DEx transcripts identified in our RNA-seq study.

We performed differential expression analysis of the RNA-seq data at both whole gene (or whole non-coding RNA) level and isoform level. We used a PBMs microarray dataset [15] to provide supporting evidence for our RNA-seq analyses results. We also used DAVID [23] to map the DEx transcripts/isoforms identified in the RNA-seq data and supported by the microarray datasets to major functional terms (e.g., GO terms) in order to gain an overall picture of the functions of the identified genes.

Overall, the majority of the transcripts DEx at the whole gene level and supported by the microarray dataset are immune system or immune disease related genes, which were all up-regulated in CP vs. control individuals. Through analysis of the DEx transcripts at the functional term (or pathway) level using DAVID [23], the term that was shared between our RNA-seq and the microarray datasets [15] is “immune response (a SP_PIR_Keyword)”. Taking the above evidence, our findings confirmed immune response as the key mechanism for CP pathogenesis. Many genes identified in this study, in spite of their functional relevance to immune response, inflammation and other aspects of periodontitis pathogenesis, have not been linked to the disease in previous studies. Hence our findings also furnished novel genetic clues to the field and contributed to better delineation of the disease mechanism at the genetic/genomic level.

Additional new information was revealed by isoform level analysis of the RNA-seq data. First, the majority (~85%) of the DEx isoforms were up-regulated in CP patients vs. controls. This observation suggests that PBMs may be more functionally active at transcriptomic level in CP status than in healthy status, consistent with the key functional roles that PBMs may play in the pathologic process of CP.

Second, through DAVID [23] analysis on the isoforms supported by the microarray dataset [15], we identified several GO terms, including endocytosis (GO:0006897), cytokine production (GO: 0001816) and apoptosis (GO: 0006915). Essentially all of the 40 component genes of the three GO terms were up-regulated in CP patients vs. controls, which may suggest enhanced monocytic endocytosis, cytokine production and apoptosis as the key potential mechanisms for CP.

The microarray dataset [15] is not a perfect dataset for supporting our RNA-seq study. However, to our knowledge, this is the only dataset that used PBMs to study gene expression changes related to periodontitis. As findings from gene expression studies are tissue/cell type dependent, it is critical that our RNA-seq study on PBMs is compared with a study using PBMs also. In terms of study design, the microarray study [15] did not compare periodontitis patients with healthy controls as in our RNA-seq study. However, the microarray study data do allow us to compare a periodontally unhealthy status (at the baseline) with a periodontally healthy status (at 10 weeks when periodontal tissue is much healthier compared with the baseline) through a paired study design [15]. From the perspective of periodontal health, our RNA-seq study had a similar study goal, which is also to compare periodontally unhealthy status (in diseased individuals) with periodontally healthy status (in healthy individuals), although this comparison was done cross-sectionally. Therefore, although we acknowledge the limitation of the microarray study [15] as a supportive dataset for our RNA-seq study, due to the above reasons, we are confident that the comparison between the two datasets is necessary and valuable and can capture “converging” evidence for transcripts associated with periodontitis, which may represent much solid scientific findings than those from the RNA-seq study alone. Indeed, the value of using the microarray dataset is evidenced by the functional relevance of the identified genes shared by the two datasets (Table 3).

This study has two major limitations. First, it has a small sample size, which leads to a limited statistical power and more genes/transcripts associated with CP may be identified with a larger sample size. Second, the sequencing depth (100bp single ended reads) is also limited and more isoforms, especially those isoforms with a low expression, may be detected with a higher sequencing depth. Therefore, the findings from this study may need to be further validated with independent studies of larger sample sizes and higher sequencing depths. Such studies are also needed to identify more genes and isoforms important to CP.

Conclusions

In summary, through RNA-seq analysis, we identified a large number of transcripts DEx in CP vs. control individuals at both whole gene and isoform levels. At the whole gene level, the identified transcripts are enriched in immune response related functional terms (Table 2). Comparison analysis using the microarray dataset for these transcripts also identified a number of novel genes related to immune response (Table 3), such as FCAR and CUX1, which were not reported for CP before. New, interesting information emerged from the isoform level analyses, including the pervasive up-regulation of PBMs gene expression, and enhanced endocytosis, cytokine production and apoptosis in CP vs. control status (Tables 5 and and6).6). Overall, our findings highlighted the significant role of PBMs for CP pathogenesis and identified a number of novel genes and potentially new mechanisms for the disease. These findings may represent novel biomarkers/signatures for clinical diagnosis and treatment of CP and provide the basis for further biological and functional investigations.

Highlights

  • The functions of monocytes in periodontitis were first explored with RNA-seq.
  • There was pervasive up-regulation of transcripts at isoform level in periodontitis.
  • Novel signatures: FACR & CUX1 genes, endocytosis, cytokine production & apoptosis

Acknowledgments

Yao-Zhong Liu, Yu Zhou and Hong-Wen Deng were partially supported by National Institute of Dental and Craniofacial Research grant RC2DE020756 and Fogarty International Center grant R03TW008221.

Abbreviation list

CP
chronic periodontitis
RNA-seq
RNA-sequencing
PBMs
Peripheral blood monocytes
DEx
Differentially expressed
DAVID
Database for Annotation, Visualization and Integrated Discovery
PPEE
posterior probability of equal expression
PPDE
posterior probability of differential expression
DCs
dendritic cells
Th
T-helper
LSUHSC
Louisiana State University Health Science Center
LSUHSC SOD
LSUHSC School of Dentistry
AA
African American
CAL
clinical attachment loss
PBMCs
peripheral blood mononuclear cells
RIN
RNA integrity number
NB
Negative Binomial
GO
Gene ontology
PostFC
posterior fold change
RA
rheumatoid arthritis
MHC
major histocompatibility complex

Footnotes

Competing interests

The author(s) declare that they have no competing interests.

Authors’ contributions

YZL: data analysis and manuscript writing

PM: coordination of subject recruitment, manuscript writing

JP: subject recruitment

YZ: monocyte isolation, RNA extraction

MB: RNA-seq data analysis

MS: RNA-seq data analysis

YPW: manuscript writing

EF: coordination of RNA-seq data analysis, manuscript writing

HWD: manuscript writing

.

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Reference List

1. Armitage GC. Development of a classification system for periodontal diseases and conditions. Ann Periodontol. 1999;4:1–6. [PubMed]
2. 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:914–920. [PubMed]
3. Ezzo PJ, Cutler CW. Microorganisms as risk indicators for periodontal disease. Periodontol 2000. 2003;32:24–35. [PubMed]
4. Taubman MA, Valverde P, Han X, Kawai T. Immune response: the key to bone resorption in periodontal disease. J Periodontol. 2005;76:2033–2041. [PubMed]
5. Pulendran B, Palucka K, Banchereau J. Sensing pathogens and tuning immune responses. Science. 2001;293:253–256. [PubMed]
6. Mosmann TR, Coffman RL. TH1 and TH2 cells: different patterns of lymphokine secretion lead to different functional properties. Annu Rev Immunol. 1989;7:145–173. [PubMed]
7. Korn T, Bettelli E, Oukka M, Kuchroo VK. IL-17 and Th17 Cells. Annu Rev Immunol. 2009;27:485–517. [PubMed]
8. Offenbacher S, Heasman PA, Collins JG. Modulation of host PGE2 secretion as a determinant of periodontal disease expression. J Periodontol. 1993;64:432–444. [PubMed]
9. Sigusch B, Klinger G, Glockmann E, Simon HU. Early-onset and adult periodontitis associated with abnormal cytokine production by activated T lymphocytes. J Periodontol. 1998;69:1098–1104. [PubMed]
10. Petit MD, Wassenaar A, van d V, van Eden W, Loos BG. Depressed responsiveness of peripheral blood mononuclear cells to heat-shock proteins in periodontitis patients. J Dent Res. 1999;78:1393–1400. [PubMed]
11. Tew JG, Burmeister JA, Palcanis KG, Ranney RR. Spontaneous lymphocyte proliferation and the periodontal status of young adults. J Periodontal Res. 1983;18:534–540. [PubMed]
12. Boyatzis S, Seymour GJ. Effect of age and periodontal disease status in man on the spontaneous proliferation of peripheral blood lymphocytes. Arch Oral Biol. 1986;31:749–755. [PubMed]
13. Davanian H, Stranneheim H, Bage T, Lagervall M, Jansson L, Lundeberg J, Yucel-Lindberg T. Gene expression profiles in paired gingival biopsies from periodontitis-affected and healthy tissues revealed by massively parallel sequencing. PLoS ONE. 2012;7:e46440. [PMC free article] [PubMed]
14. Kebschull M, Demmer RT, Grun B, Guarnieri P, Pavlidis P, Papapanou PN. Gingival tissue transcriptomes identify distinct periodontitis phenotypes. J Dent Res. 2014;93:459–468. [PMC free article] [PubMed]
15. Papapanou PN, Sedaghatfar MH, Demmer RT, Wolf DL, Yang J, Roth GA, Celenti R, Belusko PB, Lalla E, Pavlidis P. Periodontal therapy alters gene expression of peripheral blood monocytes. J Clin Periodontol. 2007;34:736–747. [PMC free article] [PubMed]
16. Armitage GC. Development of a classification system for periodontal diseases and conditions. Ann Periodontol. 1999;4:1–6. [PubMed]
17. Deng FY, Lei SF, Zhang Y, Zhang YL, Zheng YP, Zhang LS, Pan R, Wang L, Tian Q, Shen H, Zhao M, Lundberg YW, Liu YZ, Papasian CJ, Deng HW. Peripheral blood monocyte-expressed ANXA2 geneis involved in pathogenesis of osteoporosis in humans. Mol Cell Proteomics. 2011 [PMC free article] [PubMed]
18. Shen H, Qiu C, Li J, Tian Q, Deng HW. Characterization of the DNA methylome and its interindividual variation in human peripheral blood monocytes. Epigenomics. 2013;5:255–269. [PMC free article] [PubMed]
19. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323. [PMC free article] [PubMed]
20. Leng N, Dawson JA, Thomson JA, Ruotti V, Rissman AI, Smits BM, Haag JD, Gould MN, Stewart RM, Kendziorski C. EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments. Bioinformatics. 2013;29:1035–1043. [PMC free article] [PubMed]
21. Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, Speed TP. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics. 2003;4:249–264. [PubMed]
22. Smyth GK. Limma: linear models for microarray data. In: Gentleman R, Carey V, Dudoit S, Irizarry R, Huber W, editors. Bioinformatics and Computational Biology Solutions using R and Bioconductor. Springer; 2005. pp. 397–420.
23. Huang dW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44–57. [PubMed]
24. Mulcahy H, O’Rourke KP, Adams C, Molloy MG, O’Gara F. LST1 and NCR3 expression in autoimmune inflammation and in response to IFN-gamma, LPS and microbial infection. Immunogenetics. 2006;57:893–903. [PubMed]
25. Mitchell A, Rentero C, Endoh Y, Hsu K, Gaus K, Geczy C, McNeil HP, Borges L, Tedla N. LILRA5 is expressed by synovial tissue macrophages in rheumatoid arthritis, selectively induces pro-inflammatory cytokines and IL-10 and is regulated by TNF-alpha, IL-10 and IFN-gamma. Eur J Immunol. 2008;38:3459–3473. [PubMed]
26. Kuang E, Okumura CY, Sheffy-Levin S, Varsano T, Shu VC, Qi J, Niesman IR, Yang HJ, Lopez-Otin C, Yang WY, Reed JC, Broday L, Nizet V, Ronai ZA. Regulation of ATG4B stability by RNF5 limits basal levels of autophagy and influences susceptibility to bacterial infection. PLoS Genet. 2012;8:e1003007. [PMC free article] [PubMed]
27. Arlt A, Schafer H. Role of the immediate early response 3 (IER3) gene in cellular stress response, inflammation and tumorigenesis. Eur J Cell Biol. 2011;90:545–552. [PubMed]
28. Rollinger-Holzinger I, Eibl B, Pauly M, Griesser U, Hentges F, Auer B, Pall G, Schratzberger P, Niederwieser D, Weiss EH, Zwierzina H. LST1: a gene with extensive alternative splicing and immunomodulatory function. J Immunol. 2000;164:3169–3176. [PubMed]
29. Mallya M, Campbell RD, Aguado B. Characterization of the five novel Ly-6 superfamily members encoded in the MHC, and detection of cells expressing their potential ligands. Protein Sci. 2006;15:2244–2256. [PubMed]
30. Lee JS, Trowsdale J, Bodmer WF. cDNA clones coding for the heavy chain of human HLA-DR antigen. Proc Natl Acad Sci U S A. 1982;79:545–549. [PubMed]
31. Fan W, Wang Z, Kyzysztof F, Prange C, Lennon G. A new zinc ribbon gene (ZNRD1) is cloned from the human MHC class I region. Genomics. 2000;63:139–141. [PubMed]
32. Hsu WL, Tse KP, Liang S, Chien YC, Su WH, Yu KJ, Cheng YJ, Tsang NM, Hsu MM, Chang KP, Chen IH, Chen TI, Yang CS, Goldstein AM, Chen CJ, Chang YS, Hildesheim A. Evaluation of human leukocyte antigen-A (HLA-A), other non-HLA markers on chromosome 6p21 and risk of nasopharyngeal carcinoma. PLoS ONE. 2012;7:e42767. [PMC free article] [PubMed]
33. de Jong S, van Eijk KR, Zeegers DW, Strengman E, Janson E, Veldink JH, van den Berg LH, Cahn W, Kahn RS, Boks MP, Ophoff RA. Expression QTL analysis of top loci from GWAS meta-analysis highlights additional schizophrenia candidate genes. Eur J Hum Genet. 2012;20:1004–1008. [PMC free article] [PubMed]
34. Amadou C, Ribouchon MT, Mattei MG, Jenkins NA, Gilbert DJ, Copeland NG, Avoustin P, Pontarotti P. Localization of new genes and markers to the distal part of the human major histocompatibility complex (MHC) region and comparison with the mouse: new insights into the evolution of mammalian genomes. Genomics. 1995;26:9–20. [PubMed]
35. Smyth GK. Limma: linear models for microarray data. In: Gentleman R, Carey V, Dudoit S, Irizarry RA, Huber W, editors. Bioinformatics and Computational Biology Solutions Using R and Bioconductor. New York: Springer; 2014. pp. 397–420.
36. Herr AB, Ballister ER, Bjorkman PJ. Insights into IgA-mediated immune responses from the crystal structures of human FcalphaRI and its complex with IgA1-Fc. Nature. 2003;423:614–620. [PubMed]
37. Heul-Nieuwenhuijsen L, Padmos RC, Drexhage RC, de Wit H, Berghout A, Drexhage HA. An inflammatory gene-expression fingerprint in monocytes of autoimmune thyroid disease patients. J Clin Endocrinol Metab. 2010;95:1962–1971. [PubMed]
38. Achiron A, Gurevich M, Magalashvili D, Kishner I, Dolev M, Mandel M. Understanding autoimmune mechanisms in multiple sclerosis using gene expression microarrays: treatment effect and cytokine-related pathways. Clin Dev Immunol. 2004;11:299–305. [PMC free article] [PubMed]
39. Montenegro V, Chiamolera M, Launay P, Goncalves CR, Monteiro RC. Impaired expression of IgA Fc receptors (CD89) by blood phagocytic cells in ankylosing spondylitis. J Rheumatol. 2000;27:411–417. [PubMed]
40. Ikeda T, Kumagai E, Iwata S, Yamakawa A. Soluble CD26/Dipeptidyl Peptidase IV Enhances the Transcription of IL-6 and TNF-alpha in THP-1 Cells and Monocytes. PLoS ONE. 2013;8:e66520. [PMC free article] [PubMed]
41. Topic RZ, Dodig S. Eosinophil cationic protein--current concepts and controversies. Biochem Med (Zagreb) 2011;21:111–121. [PubMed]
42. Alex P, Ye M, Zachos NC, Sipes J, Nguyen T, Suhodrev M, Gonzales L, Arora Z, Zhang T, Centola M, Guggino SE, Li X. Clcn5 knockout mice exhibit novel immunomodulatory effects and are more susceptible to dextran sulfate sodium-induced colitis. J Immunol. 2010;184:3988–3996. [PMC free article] [PubMed]
43. Visekruna A, Volkov A, Steinhoff U. A key role for NF-kappaB transcription factor c-Rel in T-lymphocyte-differentiation and effector functions. Clin Dev Immunol. 2012;2012:239368. [PMC free article] [PubMed]
44. Bovin LF, Brynskov J, Hegedus L, Jess T, Nielsen CH, Bendtzen K. Gene expression profiling in autoimmune diseases: chronic inflammation or disease specific patterns? Autoimmunity. 2007;40:191–201. [PubMed]
45. Yucesoy B, Johnson VJ, Lummus ZL, Kashon ML, Rao M, Bannerman-Thompson H, Frye B, Wang W, Gautrin D, Cartier A, Boulet LP, Sastre J, Quirce S, Tarlo SM, Germolec DR, Luster MI, Bernstein DI. Genetic variants in the major histocompatibility complex class I and class II genes are associated with diisocyanate-induced Asthma. J Occup Environ Med. 2014;56:382–387. [PMC free article] [PubMed]
46. Reynolds RJ, Kelley JM, Hughes LB, Yi N, Bridges SL., Jr Genetic association of htSNPs across the major histocompatibility complex with rheumatoid arthritis in an African-American population. Genes Immun. 2010;11:94–97. [PMC free article] [PubMed]
47. Santin I, Castellanos-Rubio A, Aransay AM, Gutierrez G, Gaztambide S, Rica I, Vicario JL, Noble JA, Castano L, Bilbao JR. Exploring the diabetogenicity of the HLA-B18-DR3 CEH: independent association with T1D genetic risk close to HLA-DOA. Genes Immun. 2009;10:596–600. [PMC free article] [PubMed]