|Home | About | Journals | Submit | Contact Us | Français|
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.
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).
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.
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.
Chronic periodontitis (CP)  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 , which severely affects the quality of life in old people.
While oral bacteria are necessary for initiating pathological cascades of CP , it is host immune response and inflammatory process that is responsible for periodontal tissue damage , where monocytes play a central role. Monocytes, by differentiating into dendritic cells (DCs) in periodontium , 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 .
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 [9–12]. 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  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.
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 .
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 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 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 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.
We used RSEM v1.2.11  to analyze the RNA-seq data for alignment and expression calculation. Under RSEM, we used EBSeq  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 (where gi denote a gene or isoform i, C1 condition 1 and C2 condition 2) can be represented by a prior predictive distribution (which is constructed based on beta functions). Under the alternative (DE: differential expression), follows another prior predictive distribution . 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 and Zgi is Therefore, the posterior probability of DE for gi is defined by Bayes’ rule:
For more details, please refer to the supplement to the paper on EBSeq .
We took advantage of an online microarray dataset from GEO (GEO accession #: GSE6751)  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) . 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 . 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 . We then performed differential expression analysis using paired t test through the Bioconductor’s LIMMA (linear models for microarray data) package  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.
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/) .
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).
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.
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 . 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 . RNF5, down-regulated in CP patients vs. controls, was found to influence susceptibility to bacterial infection . 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 .
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  and LY6G5C , and 6 down-regulated transcripts, HLA-DRA , ZNRD1 , GABBR1 , RNF5 , TRIM26 , and PPP1R11 . 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.
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.
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).
Using Bioconductor’s LIMMA package , 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.
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.
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 . 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 . 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” . This gene’s expression in monocyte or mononuclear cells has been associated with several autoimmune diseases, such as autoimmune thyroid disease , multiple sclerosis  and ankylosing spondylitis . 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 . 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 . CLCN5 was found to play a role in immunopathogenesis of ulcerative colitis . 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 . VNN2’s up-regulation in mononuclear cells was associated with both inflammatory bowel disease and RA . HLA-DOA, as a MHC class II gene, has been associated with asthma , RA  and type I diabetes .
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.
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.
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).
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).
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).
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.)
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  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  to provide supporting evidence for our RNA-seq analyses results. We also used DAVID  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 , the term that was shared between our RNA-seq and the microarray datasets  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  analysis on the isoforms supported by the microarray dataset , 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  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  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 . 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  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.
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.
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.
The author(s) declare that they have no competing interests.
Authors’ contributionsYZL: 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.