Distinct gene expression characteristics in epithelial cell-Porphyromonas gingivalis interactions by integrating transcriptome analyses

Porphyromonas gingivalis is a pivotal periodontal pathogen, and the epithelial cells serve as the first physical barrier to defend the host from bacterial attack. Within this host-bacteria interaction, P. gingivalis can modify the host immune reaction and adjust the gene expression, which is associated with periodontitis pathogenesis and developing strategies. Herein, a meta-analysis was made to get the differential gene expression profiles in epithelial cells with or without P. gingivalis infection. The network-based meta-analysis program for gene expression profiling was used. Both the gene ontology analysis and the pathway enrichment analysis of the differentially expressed genes were conducted. Our results determined that 290 genes were consistently up-regulated in P. gingivalis infected epithelial cells. 229 gene ontology biological process terms of up-regulated genes were discovered, including “negative regulation of apoptotic process” and “positive regulation of cell proliferation/migration/angiogenesis”. In addition to the well-known inflammatory signaling pathways, the pathway associated with a transcriptional misregulation in cancer has also been increased. Our findings indicated that P. gingivalis benefited from the survival of epithelial cells, and got its success as a colonizer in oral epithelium. The results also suggested that infection of P. gingivalis might contribute to oral cancer through chronic inflammation. Negative regulation of the apoptotic process and transcriptional misregulation in cancer pathway are important contributors to the cellular physiology changes during infection development, which have particular relevance to the pathogenesis and progressions of periodontitis, even to the occurrence of oral cancer.


Introduction
Periodontitis is a ubiquitous inflammatory disease and the primary cause of tooth extraction in adults resulted from infection of periodontopathic bacteria. The host manifests an inflammatory immune response to this bacterial infection. Porphyromonas gingivalis is an opportunistic pathogen and mainly colonized in periodontal tissues [1]. More and more evidence shows that it is also a pathogen of systemic diseases. P. gingivalis may modify gene expression, even host immune response by degrading host cell surface proteins or receptors [2,3]. Within this host-bacteria balance, the gene expression is directly correlated with individual susceptibility and influences the pathogenesis of periodontitis and disease progression. Several studies have suggested that plenty of genes were regulated by a host-protective response in gingival tissue of periodontitis. These related genes are critical components of immunological response, biological behavior, and metabolic signal pathways. Recently Kebschull M et al. have employed microarray analyses to speculate the gene expression signatures in periodontitis [4][5][6][7][8]. It is proved that microarray could supply more insight into the etiology of periodontitis.

Ivyspring International Publisher
The growing microarray data offer some valuable clue for analyzing the host immune response of periodontitis. However, the discrimination against critical genes and relevant pathways based on these reports were limited because of the sample size, differences in research design and reporting methods in the separate studies.
A meta-analysis deals with the transcriptome data by combining the individual microarray study, and counters the above-mentioned disadvantage [9]. In the current study, a meta-analysis consisting of three separate microarray datasets was made to distinguish the gene expression profile in the epithelial cells with or without P. gingivalis infection. The present study gives us a synthesized assessment of the gene expression signatures in the epithelial cells with P. gingivalis infection. The differentially expressed genes (DEGs), the gene ontology (GO) terms and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways concerned in the transcription signatures were identified [10][11][12][13]. Our analysis would offer new understanding in the comprehension of periodontitis pathobiology, more information on the design of future research, and more tactics to block the progression of periodontitis.

Gene expression array data collection
An electronic search was performed in Gene Expression Omnibus (GEO, NCBI) database. Microarray data of the gene expression with the keywords "Porphyromonas gingivalis" or "gingival epithelial cells" or "oral epithelial cells" were downloaded. Epithelial cells infected with P.gingivalis were considered as "case group", while non-infected epithelial cells as "control group". Sample-sourced datasets from epithelial cells infected with other bacteria were excluded. Three separate microarray datasets were enrolled with raw data. The details of these datasets were listed in Table 1.
The basic data were obtained from the 3 individual studies, including GEO accession; sample source; platform; numbers of cases or controls in the array data. Two of the datasets were performed in Affymetrix Human Genome U133 Array. The third dataset was conducted in Affymetrix Human OElncRNAs520855F Array. This analysis contains totally 22 GEO transcriptome datasets for epithelial cells infected with P.gingivalis or non-infected control. Microarray data for epithelial cells with P.gingivalis infection (n = 11) and non-infected control (n = 11) from GEO transcriptome database (http://www.ncbi. nlm.nih.gov/geo/) were compared for their transcriptome profiles.
The research was authorized by the Ethics Committee of School of Stomatology, China Medical University (Shenyang). They determined the current analysis of the public data set did not require the consents of the patients.

Differential gene expression analysis
DEGs were identified in epithelial cells with P.gingivalis infection or not, the datasets extracted from the individual microarray assays were submitted to the network-based meta-analysis program for gene expression profiling (http://www.inmex.ca/INMEX/), starting with multiple gene expression tables for meta-analysis [13][14][15][16]. First, the gene ID or probe ID was converted to its identical Entrez ID using the gene/probe conversion tool in INMEX. Second, the data were enrolled, processed, annotated and the intensity measurement of every ID was log2 transformed and normalized to zero mean and unit variance. While performing differential expression analysis on individual data set, the False Discovery Rate (FDR) of Benjamini-Hochberg's was set 0.05 in order to adjust the cut-off P-value. Third, data integrity was assessed as described previously [17][18]. Here we used the batch effect correction option so as to minimize the potential batch effect [19]. Then the meta-analysis for DGEs was made exploiting the INMEX web-based computational tool. Statistical analysis was performed by the INMEX program. The combined P values were detected according to Fisher's method. A P value of less than 0.05 was considered as the criterion of statistically significant in our study. The genes selected were classified based on the grade products of the combined P-value. The statistics were downloaded and visual exploration was performed. A Venn diagram was made to compare the difference between the findings of the meta-analysis and the individual study (http://bioinformatics.psb.ugent .be/webtools/Venn/).

Gene Ontology Biological Process terms and KEGG pathway analysis
Subsequently, we exploited a web-based tool (named Database for Annotation, Visualization, and Integrated Discovery) to conduct the GO Biological Process (GO_BP) terms and KEGG pathway analysis [9,18,20]. DAVID Version 6.8 was used here. In this database, genes were divided into different classes according to their biological processes or molecular functions. Using this GO analysis we compared the DEGs and distributed them into a functional systematization. Gene Entrez ID was uploaded and analyzed for its GO Biological Process Annotation with functional annotation chart after the identifier was selected. Pathway annotations of the DEGs were acquired from the KEGG database [12,19]. Pathway categories with a Benjamini-corrected P-value < 0.05 were utilized to determine a critical analysis. Our figures listed the typical GO biological process terms, which were chosen from the functional annotation charts that were significantly enriched at the top. Representative KEGG pathways chosen from the most significantly enriched charts were shown in our figures.

DEGs in infected and non-infected epithelial cells
First of all, we compared the genome-wide gene expressions of infected and non-infected epithelial cells across microarray datasets. There were 362 genes differentially expressed in the two groups of epithelial cells. Among the 362 DEGs, 290 genes were significantly up-regulated (Table S1), and 72 genes were significantly down-regulated (Table S2). The up-regulated DEGs consisted of PLK2, CYP24A1, TNFAIP3, PTGS2, EDN1, IL6, GADD45A, IL1B, etc. The top 50 up-regulated DEGs were shown in the heatmap (Fig. 1).
The most consistently up-and down-regulated DEGs of top 20 in epithelial cells, in the comparison of infected and non-infected epithelial cells, were shown in Table 2.
We compared the difference between the findings of the meta-analysis and the separate study, and a Venn diagram was made. 285 DEGs were distinguished in the meta-analysis only, which called gained DEGs. While only 16 DEGs were found in the separate study only, which called lost DEGs (Fig. 2).

Functional classification and pathway assignment of DEG
GO analysis was carried out, and the DEGs were classified into different hierarchical categories based on the GO database. 229 GO Biological Process terms of up-regulated genes were discovered (Table S3). It has been shown that the significantly enriched GO terms for the DEGs included "negative regulation of apoptotic process", "aging", "positive regulation of cell proliferation", "positive regulation of cell migration", and "positive regulation of angiogenesis". The top 20 enriched GO terms for up-regulated DEG in epithelial cells were shown in Fig. 3.
Here, we paid more attention to the first GO terms, namely negative regulation of apoptotic process. And the involved DEGs associated with up-regulated biological process categories of negative regulation of apoptotic process were listed in Table 3.
KEGG pathway analysis was made in order to annotate the functional categories of DEGs as well.
Here, 38 up-regulated pathways were found (Table  S4). The 20 most significantly up-regulated pathways for DEGs were demonstrated in Fig. 4.
Besides the well-known inflammatory signaling pathways, such as nuclear factor kappa-light-chainenhancer of activated B cells (NF-κB) signaling pathway, tumor necrosis factor (TNF) signaling pathway, mitogen-activated protein kinase (MAPK) signaling pathway, we found that the pathway associated with a transcriptional misregulation in cancer, Janus kinases/signal transducer and activator of transcription (JAK/Stat) pathway, and transforming growth factor-β (TGF-β) signaling pathway, had also been increased. Representatively, the involved DEGs of the KEGG pathways of transcriptional misregulation in cancer were listed in Table 4.

Discussion
As a Gram-negative anaerobic bacterium, P. gingivalis is one of the most prominent pathogens in the etiology of chronic periodontitis. Epithelial cells are now considered to be the first immune barrier in the pathogenesis of periodontal disease confronted by the microorganism that colonizes or even intrudes mucosal surfaces. The intercommunication of P. gingivalis with epithelial cells involves several intricate signaling pathways [21][22][23]. Recognition of the pertinent genes and signal pathways concerned is pivotal, so we could comprehend the cellular and molecular changes during the periodontitis onset and disease progression. Considering P. gingivalis can not only invade the local periodontal tissue but also persist in other systems of the human being (for example esophagus tissue, respiratory and vascular), we should be more cautious in the active performance of P. gingivalis with epithelial cells. Now we have obtained the transcriptomic profile with our meta-analysis using several public microarray data. Taking advantage of this enhanced statistical power, our results were more reliable [24][25]. We identified the DEGs, the corresponding Biological Process and pathways in P. gingivalis-infected epithelial cells based on comparisons to the GO [8,9] and KEGG [12][13] databases, respectively. The present results indicated the changes in gene expression feature modulated by some signal pathways were related to host-bacterial balance.  In our meta-analysis we showed several genes about negative regulation of the apoptotic process, including BCL2A1, BCL2, BARD1, CFLAR, were significantly up-regulated in epithelial cells infected with P. gingivalis compared to non-infected ones. It is conceived that P. gingivalis improves the survival of host cells by directly inhibiting different pro-apoptotic pathways [3,[26][27][28]. Inducing the durability of epithelial cells is crucial for the survival of P. gingivalis in these cells. It is reported that this resistance to apoptosis may be related to the BCL2 protein family, including BCL-2 and BAX [29]. This suppressed apoptosis during P. gingivalis infection was phosphatidylinositol 3 kinase/Protein Kinase B (PI3K/Akt) -dependent, furthermore, the PI3K signal pathway could regulate transcription or post-translational modification of certain BCL2 family member [3].
Our results verified that P. gingivalis infection has numerous anti-apoptotic effects on epithelial cells (such as activation of JAK/Stat pathways, NF-κB signal pathway, etc). The diversity of these multiple genes and pathways demonstrates that P. gingivalis may affect apoptotic pathways at different levels. Some of the up-regulated genes and pathways surround the mitochondria, which is a critical component of the intrinsic apoptosis process. The analysis of our study, along with others reports, showed that P. gingivalis suppressed the apoptosis of epithelial cell partially through the JAK/Stat signal pathway, which regulates the inherent mitochondria cell death mediated by the caspase-dependent apoptotic pathway.
At the meantime, P. gingivalis participates in the regulation of the extrinsic apoptotic pathway. Our KEGG pathway analysis demonstrated NF-κB signal pathway was up-regulated. It is reported that NF-κB, together with TNF and FasL, takes its role in the composition of the extrinsic apoptotic process [30]. NF-κB activation serves as a primary mechanism to protect cells against apoptotic stimulus such as TNF [31]. Current findings and early reports verify that the ability of P. gingivalis to benefit the survival of epithelial cells may assist its successful colonization in the oral mucosal epithelium.
IL6 is regarded as a tumor-promoting cytokine in diverse kinds of cancer including oral squamous cell carcinoma (OSCC) [32][33][34][35]. It is considered as one of the most promising markers for early diagnosis of tongue squamous cell carcinoma, which was the most usual kind of OSCC. Geng et al. revealed that both CEBPB and IL6 were up-regulated in late OSCC and CEBPB was activated as an upstream regulator of IL6. Moreover, CEBPB was recently identified as one of the master regulators in cancer biology [35,36].
Nowadays, more researches have reported a significant increase of CXCL8 and its mRNA in P. gingivalis-infected epithelial cells or oral squamous cells [37][38][39]. Ha's study found that infection by P. gingivalis stimulated the tumorigenic features and invasiveness of oral squamous cells, which was correlated with enhanced expression of CXCL8 [37]. PLAU, the serine protease, has been indicated in promoting pericellular proteolysis and oncogenic signal transduction in cancer cells. Yoshizawa implied that cases with PLAU/PLAU receptor/maspin expression should be considered as the high risk of poor prognosis of OSCC, which may lead to severe invasiveness, as well as the increased risk for cervical lymph node metastasis [40].
In addition, over-expression of MEIS1 in acute leukemia has been observed [41], which is related to the shorted latency and accelerated progression of acute leukemia [42]. The enhancer E9 mediates an auto-regulatory loop which contributes to the consistent expression of MEIS1 in leukemia. It was reported that MPO levels were dramatically enhanced in the bronchoalveolar lavage fluid and positively related to the carcinogenicity of inflammatory in lung cancer formation [43][44]. MPO is conceived as a critical mediator of inflammatory-promoted tumorigenesis, and its pro-tumor effects may take place in the window just after tumor onset in the lung [45].
So far, several studies have attempted to elucidate the role of P. gingivalis in contribution to oral cancer. Previous research by our group showed P. gingivalis enhanced the cell proliferation, migration, and invasive ability of human oral epithelial cells. Thus it promoted the cells tumorigenic properties [2]. Gonda et al. reported that chronic inflammation induced by bacterial infection drove carcinogenesis, involving induction, progression, invasion, and metastasis [33]. It is believed that the inflammatory microenvironment was an integral part of oral cancer [2,[32][33][34][35]46,47].
Our KEGG pathway analysis further confirmed P. gingivalis played a role in promoting tumor transformation in some inflammatory conditions. Here pathways enrichment investigation showed that several pathways were enhanced significantly, including the NF-κB pathway, TNF pathway, TGF-beta pathway, MAPK signal pathway, and cytokine-cytokine receptor interaction. These inflammation-related signaling pathways were listed as the top 5 of our KEGG pathway analysis results.
Here we believe that this chronic bacteria infection should be identified as probable effector which may result in oral cancer.
Generally speaking, the current meta-analysis using microarray datasets offered us an outlook about the gene expression characteristics in epithelial cells infected with P.gingivalis. The bioinformatics analyses of the up-regulated genes showed us that the GO terms of the negative regulation of apoptotic process, several inflammation-related signaling pathways, and the transcriptional misregulation in cancer were significantly enriched. It provided us some new sights into the alteration in cytopathogenicity that took place in the process of infection, which had particular relevance to periodontitis pathogenesis and developing strategies. Meanwhile, it provided more evidence to support the view that P.gingivalis might be related to the occurrence and development of OSCC. Genomes; NF-κB: nuclear factor kappa-light-chain-enhancer of activated B cells; OSCC: oral squamous cell carcinoma; TGF-β: transforming growth factor-beta; TNF: tumor necrosis factor.

Supplementary Material
Supplementary