Circulating lymphocytes and monocytes transcriptomic analysis of patients with type 2 diabetes mellitus, dyslipidemia and periodontitis

Type 2 diabetes mellitus (T2DM), dyslipidemia and periodontitis are frequently associated pathologies; however, there are no studies showing the peripheral blood transcript profile of these combined diseases. Here we identified the differentially expressed genes (DEGs) of circulating lymphocytes and monocytes to reveal potential biomarkers that may be used as molecular targets for future diagnosis of each combination of these pathologies (compared to healthy patients) and give insights into the underlying molecular mechanisms of these diseases. Study participants (n = 150) were divided into groups: (H) systemically and periodontal healthy (control group); (P) with periodontitis, but systemically healthy; (DL-P) with dyslipidemia and periodontitis; (T2DMwell-DL-P) well-controlled type 2 diabetes mellitus with dyslipidemia and periodontitis; and (T2DMpoorly-DL-P) poorly-controlled type 2 diabetes mellitus with dyslipidemia and periodontitis. We preprocessed the microarray data using the Robust Multichip Average (RMA) strategy, followed by the RankProd method to identify candidates for DEGs. Furthermore, we performed functional enrichment analysis using Ingenuity Pathway Analysis and Gene Set Enrichment Analysis. DEGs were submitted to pairwise comparisons, and selected DEGs were validated by quantitative polymerase chain reaction. Validated DEGs verified from T2DMpoorly-DL-P versus H were: TGFB1I1, VNN1, HLADRB4 and CXCL8; T2DMwell-DL-P versus H: FN1, BPTF and PDE3B; DL-P versus H: DAB2, CD47 and HLADRB4; P versus H: IGHDL-P, ITGB2 and HLADRB4. In conclusion, we identified that circulating lymphocytes and monocytes of individuals simultaneously affected by T2DM, dyslipidemia and periodontitis, showed an altered molecular profile mainly associated to inflammatory response, immune cell trafficking, and infectious disease pathways. Altogether, these results shed light on novel potential targets for future diagnosis, monitoring or development of targeted therapies for patients sharing these conditions.

index insulin resistance (Table 1). Insulin levels showed no significant differences between T2DMpoorly-DL-P and T2DMwell-DL-P groups 6 , as well as the BMI and abdominal circumference were similarly higher in the diabetic groups (T2DMpoorly-DL-P and T2DMwell-DL-P) compared with normoglycemic patients (H; Table 1). As expected, the total cholesterol, LDL, and triglyceride levels were higher in the T2DMpoorly-DL-P, T2DMwell-DL-P, and DL-P groups compared to P and H groups. Additional information can be found in our previous studies 14, 25 .
T2DM patients (T2DMpoorly-DL-P and T2DMwell-DL-P) presented P in a moderate to severe level, since they showed 4-5 mm of probing depth in ≥20% of the periodontal sites; 3-4 mm of clinical attachment loss in ≥30% of the periodontal sites, and ≥40% of the periodontal sites with bleeding on probing 14,25 . Both groups demonstrated similar periodontal tissue destruction, with bone loss and local inflammation. The T2DMwell-DL-P group showed a significant difference regarding the presence of deeper periodontal sites in comparison to groups without T2DM (DL-P, P and H) 6 . Patients from the T2DMpoorly-DL-P and T2DMwell-DL-P groups presented worse periodontal clinical condition than the P and DL-P 14,25 (Table 1).

Gene expression analysis and functional enrichment analyses by Ingenuity Pathway Analysis (IPA) and Gene Set Enrichment Analysis (GSEA).
In the circulating lymphocytes and monocytes gene expression analysis, when considering the T2DMpoorly-DL-P versus H comparison, we identified 1374 up-and down-regulated DEGs, while 869 were identified in the T2DMwell-DL-P versus H, 521 up-and down-regulated DEGs (DL-P versus H) and 564 up-and down-regulated DEGs (H versus P). The segregation of these pairwise comparisons of subjects by hierarchical cluster analysis is shown in Supplementary Figure 1.
To get a better insight into the biological function of the genes regulated by these inflammatory diseases, we performed Ingenuity Pathway Analysis (IPA) on all differentially expressed genes. Table 2 summarizes the networks identified in each pairwise group comparison, while the Top Canonical Pathways are found in the Supplementary Table 1. Noteworthy, Table 2 shows repeatedly networks associated with immunological pathways. Some of these enriched networks are illustrated in Figs. 1 and 2.
We also performed the gene set enrichment analysis (GSEA) in the same group comparisons as in the IPA analysis. The top 20 enriched (curated) gene sets in the circulating lymphocytes and monocytes of each group comparison by GSEA are shown in the Table 3. In both comparisons, T2DMwell-DL-P vs. H and DL-P vs. H) there were fewer than 20 statistically significant enriched gene sets. In Table 3 we highlighted in gray the gene sets that are unique in each group comparison. From these unique gene sets, we selected in bold the gene sets that  www.nature.com/scientificreports www.nature.com/scientificreports/ included genes enriched in the disease groups that were also found in the top IPA networks (Table 2). For example, in Table 3 the two selected gene sets (bold and gray highlighted) that are related to interferon responsive genes and interferon alpha and beta signaling genes where also found in the lists of the top IPA networks. The plots of those selected gene sets are shown in Fig. 3.
The Venn diagram in Fig. 4 exhibits the up-regulated differentially expressed genes that were unique in circulating lymphocytes and monocytes of each disease group compared to the healthy control group. Analyzing the Fig. 4 in combination with the Table 4, we observed that the comparison of T2DMwell-DL-P vs. H revealed 177 up-regulated genes. From those 177 up-regulated genes, 17 were related to immune response. Table 4 shows that in the first three group comparisons, as the metabolic impairment of the patient decreases, the number of up-regulated genes related with the immune response also decreases.
Validation by RT-qPCR of selected genes from the network analysis. Among the circulating lymphocytes and monocytes top up-and down-regulated genes present in the enriched biological pathways, we selected three to five genes from each disease group to validate by RT-qPCR. The results are shown in Figs. 5 and 6. The results of RT-qPCR comparing T2DMpoorly-DL-P vs. H, which assessed the influence of poor glycemic control of T2DM added to the presence of dyslipidemia and periodontitis, validated the importance of the TGFB1I1 (transforming growth factor beta 1 induced transcript 1) and the VNN1 (vanin 1) genes, since they were upregulated on T2DMpoorly-DL-P subjects, while the HLADRB4 (major histocompatibility complex, class II, DR beta 4) and the CXCL8 (C-X-C motif chemokine ligand 8) genes were downregulated in the same group (Figs. 5A and 6C).
On the T2DMwell-DL-P vs. H comparison, which evaluated the stimulus of T2DM with good glycemic control and the presence of dyslipidemia and periodontitis, the RT-qPCR results from circulating lymphocytes and monocytes validated that the BPTF (bromodomain PHD finger transcription factor) and PDE3B (phosphodiesterase 3B) genes were upregulated in H, and the FN1 (Fibronectin 1) gene was downregulated in the same group (Fig. 5B).
On the DL-P vs. H comparison, which evaluated the effect of the dyslipidemia and periodontitis, the qPCR confirmed that the HLADRB4 and CD47 (CD47 molecule) genes were upregulated in H, while the DAB2 (DAB adaptor protein 2) gene was downregulated in H (Fig. 6A,D). Lastly, the P vs. H comparison, which evaluated the association with periodontitis, the RT-qPCR from circulating lymphocytes and monocytes confirmed that the HLADRB4 and ITGB2 (integrin subunit beta 2) genes were upregulated in H, while the IGHDL-P (immunoglobulin heavy constant gamma 3) gene was downregulated in H (Fig. 6B,E).

Discussion
To the best of our knowledge, this is the first study describing functional information related to the circulating lymphocytes and monocytes transcriptional profile of individuals affected (simultaneously or not) by three of the currently most prevalent diseases: T2DM, dyslipidemia and periodontitis. Using a gene expression dataset generated from PBMC followed by functional enrichment analyses, we showed that validated DEGs were implicated  www.nature.com/scientificreports www.nature.com/scientificreports/ mainly in immune cell trafficking, antimicrobial and inflammatory response, cell-to-cell signaling and interaction, infectious and cardiovascular diseases, cellular growth and proliferation.
When we evaluated the circulating lymphocytes and monocytes expression profiling of T2DMpoorly-DL-P vs H, we sought to find DEGs responsible for the influence of poor glycemic control of T2DM in the presence of dyslipidemia and periodontitis. The RT-qPCR results demonstrated that the TGFB1I1 and VNN1 genes were upregulated in T2DMpoorly-DL-P in comparison to H. Functional information of the VNN1 gene can be found in the Supplementary Discussion section. Based on the functional network analysis of the microarray data, these genes are associated with cellular signaling and movement, cancer, molecular transport, vitamin/mineral metabolism, organism injury and abnormalities. Considering the upregulated genes in this pairwise comparison, the TGFB1I1 gene encodes a secreted ligand of the TGF-beta (transforming growth factor-beta) superfamily of proteins. Ligands of this family bind various TGF-beta receptors leading to recruitment and activation of the SMAD family, which are signal transducers, and transcriptional modulators that mediate the TGF-beta signal and thus regulate multiple cellular processes 26 . The protein encoded by the TGFB1I1 gene regulates cellular proliferation, differentiation and growth, and modulate the expression and activation of other growth factors including interferon gamma (IFN-γ) and tumor necrosis factor alpha (TNF-α) 26,27 . Our present results lead to the hypothesis that the circulating TGFB1I1 overexpression in T2DMpoorly-DL-P subjects could contribute to an exaggerated cellular activation and immune response compared to individuals in H. To this date, we were unable to find previously published data showing the circulating lymphocytes and monocytes expression of this gene in the three pathologies studied here. Based on these findings, further studies should be designed and performed to  Table 2. Straight lines represent direct interaction, and dotted lines represent indirect interaction.  Table 2. Straight lines represent direct interaction, and dotted lines represent indirect interaction. (2020) 10:8145 | https://doi.org/10.1038/s41598-020-65042-9 www.nature.com/scientificreports www.nature.com/scientificreports/  www.nature.com/scientificreports www.nature.com/scientificreports/ validate and investigate the effects on peripheral cells and organs affected by these diseases. Otherwise, increased TGFB1I1 mRNA expression level was identified in individuals exhibiting rheumatoid arthritis in comparison with individuals with osteoarthritis, both conditions presenting the fibroblast like synoviocytes 28 .
Interestingly, the CXCL8 (C-X-C Motif Chemokine Ligand 8) gene was downregulated in T2DMpoorly-DL-P compared to H subjects. The CXCL8 gene, previously named as interleukin 8 (IL-8), encodes the protein that is a key chemokine in the initiation and amplification of acute inflammatory reactions and in chronic inflammatory processes because it attracts and activates neutrophils in inflammatory regions 29,30 . Regarding periodontitis, there is an interplay between microbial species within subgingival biofilms and the adjacent periodontal tissues 6 , which elicit the innate immunity by releasing chemokines, such as the CXCL-8/IL-8 into the gingival crevicular fluid (GCF) 31 necessary to recruiting neutrophils. We hypothesize that there should be an equilibrated CXCL8/IL8 mRNA and protein production to obtain adequate host immune response against periodontal pathogens driving to a healthy periodontal status. In spite to contradictory studies regarding the salivary IL-8 levels in periodontitis patients Khalaf et al. 32 reported significantly higher IL-8 levels in saliva from periodontally healthy individuals in comparison with those affected by periodontitis. This result is in agreement with a meta-analysis that showed higher levels of IL-8 (pg/µL) in the GCF of periodontally healthy control subjects compared with periodontitis patients 33 . We also found that the circulating lymphocytes and monocytes CXCL8/IL8 mRNA levels of periodontally healthy subjects were significantly lower when compared to patients affected concomitantly by T2DM, dyslipidemia and periodontitis. In addition, to the evidences of the systemic higher expression of pro-inflammatory cytokines, such as IL-1β, IL-6 and TNF-α in PBMC of T2DM patients [34][35][36][37] , our results suggest that the circulating lymphocytes and monocytes present a hyper-inflammatory state, together with lower levels of the CXCL8/IL8 expression, might contribute to the disturbance in the orchestration of the immune response in patients affected by these diseases.
Interestingly, the IPA revealed in a unique network (Fig. 1A) the downregulation of CXCL8/IL8 gene in T2DMpoorly-DL-P compared to H subjects contrasting with the upregulation of the TGFB1I1 gene in the circulating lymphocytes and monocytes of the same subjects. To our knowledge, it is the first time that both CXCL8/ IL8 and TGFB1I1 genes are found in blood peripheral cell of patients with this combination of diseases and significantly contrasting to each other. Obviously, since the mechanism by which this network might affect the characteristics of the investigated three diseases is largely unknown, further studies are needed to explain these    Table 3.  Table 4. www.nature.com/scientificreports www.nature.com/scientificreports/ The findings of Fig. 6A that showed the validated downregulation of the CD47 gene and upregulation of the DAB2 gene in circulating lymphocytes and monocytes of DL-P patients are in agreement with the IPA analysis ( Fig. 2A). These genes are correlated on the same network associated with cell-to-cell signaling and interaction, cell morphology and cellular movement. The IPA analysis also showed the validated HLADRB4 and ITGB2 genes both downregulated in Periodontitis patients (group P), interacting via TNF on the same network (Fig. 2B) of cellular compromise, inflammatory response, and infectious diseases. Finally, the HLADRB4 gene was validated in peripheral blood of T2DMpoorly-DL-P, DL-P and P versus H comparisons, and was always found to be upregulated in H (systemic and oral healthy subjects) (Figs. 4E and 6C-D)   www.nature.com/scientificreports www.nature.com/scientificreports/ this class II molecule is a heterodimer consisting of an alpha (DRA) and a beta (DRB) chain, both anchored in the membrane 38 . Class II molecules are expressed in antigen presenting cells (B lymphocytes, dendritic cells, macrophages) 39 , playing a central role in the immune system by presenting peptides derived from extracellular proteins 40 . Our present results lead to the hypothesis that the HLADRB4 expression levels indicate the healthy status of the patients, since it was consistently upregulated in circulating lymphocytes and monocytes of systemically healthy individuals without periodontitis (H). Further studies will need to investigate a potential correlation between circulating high levels of HLADRB4 and low levels of pro-inflammatory expression in individuals with some chronic inflammatory disease, such as poor metabolic controlled T2DM, dyslipidemia and periodontitis.
The Venn diagram (Fig. 4) shows the number of genes specifically upregulated in circulating lymphocytes and monocytes of each group comparisons. Deeper analyses of those findings, followed by further investigations to assess their validity and functionality, might be useful to indicate blood biomarkers for diagnosis of the combined pathologies, for monitoring those patients in treatment, or to provide insights into novel therapeutic strategies for disease prevention and treatment.
Enriched gene sets detected by the GSEA included genes enriched in the circulating lymphocytes and monocytes of diseased groups that were also find in the top IPA networks. Therefore, taking together the functional enrichment analyses by IPA and GSEA, we observed many similarities regarding biological functions identified for each group comparison, but the peripheral blood molecular signaling specificities demand deeper analyses, as well as additional in vitro and in vivo experiments. Comparing the Venn diagram and the gene ontology findings showed in Table 4, the comparison T2DMpoorly-DL-P vs. H group showed twice the number of up-regulated genes related with the immune response (37) in comparison to the 17 up-regulated genes in the same category in the T2DMwell-DL-P vs. H group. Even though both groups are affected by T2DM and it is recognized that these patients reflect a hyperinflammatory immune response state 41 , our results suggest that the patients with increased metabolic impairment can elicit the immune response in a broader and different way than the patients presenting adequate metabolic control. Nonetheless, the up-regulated group-specific genes expressed in circulating lymphocytes and monocytes as showed in the Venn diagram (Fig. 4) should be further investigated.
Despite the significant amount of novel data showed here, the main limitations of this study are (i) the transcriptome signature presented here is specific to circulating lymphocytes and monocytes; therefore the transcriptome signature of any organ or any other cell type (e.g. liver, pancreas, kidney, periodontium, neutrophils, etc) was not included. Consequently, the results showed here should be placed in the context of the circulating lymphocytes and monocytes, which in fact, was a main goal of our study. We acknowledge the limitations to discuss how the circulating lymphocytes and monocytes, and their molecular signature, may impact how these cells home and act in target www.nature.com/scientificreports www.nature.com/scientificreports/ tissues of T2DM, Dyslipidemia and Periodontitis; (ii) the lack of longitudinal data tracking the disease progression in our cohort of patients, and (iii) the absence of genetic analyses at the translational level. Further studies enrolling additional larger and ethnically diverse population of patients are important to search for external validation of the present findings. Moreover, the present study was not designed to assess whether the disease phenotype of each group of patients was the cause or consequence of the gene expression signature. Of importance, even though we believe that the PBMC transcriptome data herein presented is novel and relevant, we clearly acknowledge its limitations when compared to big-data approaches, such as the Schüssler -Fiorenza Rose et al. 42 , multi-omics studies including the microbiome such as the Zhou et al. 43 , or other large-scale dataset methodologies developed to integrate the transcriptome within the context of the interactome, similarly to Li et al. 44 , which developed a PBMC transcriptome analysis, and constructed an interactome to investigate the multi-layered regulatory pathways in T2DM.
In addition to big-data and host-microbiome interaction studies, methods such as transcriptome, proteome and interactome of other cell types could contribute to further elucidate and/or corroborate our findings. Therefore, since T2DM and dyslipidemia are complex multi-system diseases, our present PBMC transcriptome is able to demonstrate one layer of the multi-omics signatures that could be investigated. Certainly, further studies are important to reach the biological mechanisms mediating the association of the DEGs in PBMC with each specific pathological combination, preferably in a broader and deeper way, such as by using multi-omics strategies to understand the interactions of the different biological aspects of each combination of complex diseases.
Here we identified the differences in the circulating lymphocytes and monocytes expressed genes in patients with different combinations of T2DM, DL and P. We believe that our results represent an important advance in the field, and bring to light a diverse range of molecules in the PBMC that can be further investigated for their validity as biomarkers for a combination of pathologies as well as their potential functional role in the physiological and pathological environment. Strong aspects of this study include (i) the selected patients were subjected to strict clinical inclusion criteria accurately distinguishing the pathological condition between patients, permitting to identify specific PBMC DEGs related to each of these clinical conditions; (ii) the data was consistently validated for the PBMC expression levels by RT-qPCR method, demonstrating the internal validity of the study.
In conclusion, we identified and validated a gene expression signature for circulating lymphocytes and monocytes from T2DM (poorly or well-controlled), dyslipidemia and periodontitis patients and showed the most relevant functional pathways of each pathological combination. Additional studies are needed to externally validate our findings and evaluate whether these genes can in the future be used as peripheral blood markers for risk assessment of each complex disease combination, as well as can provide novel insights into targeted therapeutic strategies for disease prevention and/or treatment. www.nature.com/scientificreports www.nature.com/scientificreports/

Materials and Methods
Study population. All experimental protocols were approved by the Ethics in Human Research Committee of School of Dentistry at Araraquara (UNESP -São Paulo State University, Araraquara, Brazil; Protocol number 50/06), and all the volunteers signed an informed consent. All methods were carried out in accordance with the principals of the Declaration of Helsinki.
From 2009 to 2011, we evaluated 1788 patients; the study's inclusion criteria included age range from 35 to 60 years, presence of at least 15 natural teeth and similar socio-economic level 6 . We analyzed 150 patients who were divided into five groups of 30 patients each, based upon diabetic, dyslipidemic and periodontal status: poorly controlled T2DM with dyslipidemia and periodontitis (T2DMpoorly-DL-P); well-controlled T2DM with dyslipidemia and periodontitis (T2DMwell-DL-P); normoglycemic individuals with dyslipidemia and periodontitis (DL-P); systemically healthy individuals with periodontitis (P); and systemically healthy individuals without periodontitis (H) 14 . The exclusion criteria were previously described in Corbi et al. 25 .
Physical, biochemical and periodontal evaluations. Each subject was submitted to physical examination including anthropometric data such as abdominal circumference (cm), hip (cm), waist (cm) and height (m), weight (kg), and body mass index (BMI) 6 .
Blood samples were collected after a 12-hour overnight fast for the evaluation of fasting plasma glucose (mg/ dL) by modified Bondar & Mead method, glycated haemoglobin (HbA1c) by enzymatic immunoturbidimetry, insulin levels by the chemiluminescence method (U/L) and high-sensitivity C-reactive protein by the nephelomeric method 6 . Insulin resistance was evaluated by calculation of the homeostasis model assessment (HOMA). The same laboratory performed all analyses. Patients were considered as nondiabetics (normoglycemic individuals) if they presented fasting glucose levels <100 mg/dL and HbA1c < 6.5%. T2DM patients were diagnosed by an endocrinologist who monitored their glycemic levels by evaluation of HbA1c; patients were divided into poorly controlled (HbA1c ≥ 8.5%, 64 mmol/mol) or well-controlled patients (HbA1c < 7.0%) 14 . Diabetic individuals with HbA1c between 7.0-8.5% were not included.
Lipid profile [total cholesterol (TC), triglycerides (TGs), and HDL] was performed by enzymatic methods. LDL was determined by the Friedewald formula. To avoid the inclusion of individuals with transitory dyslipidemia, the cutoff points used here were the highest values according to the National Cholesterol Educational Program (NCEP) Adult Treatment III (ATP III) 45 : TC ≥ 240 mg/dL, LDL ≥ 160 mg/dL, HDL ≤ 40 mg/dL, and TGs ≥200 mg/dL 6 .
Diagnosis of periodontal disease as defined by the American Academy of Periodontology 46 includes local signs of inflammation and tissue destruction (presence of deep periodontal pockets ≥6 mm) and loss of the connective tissue attachment of gingiva to teeth (clinical attachment loss ≥4 mm) in at least 4 non-adjacent teeth. All patients were subjected to a periodontal clinical examination performed in six sites per tooth by a single trained calibrated examiner (A.B.S., Kappa = 0.89) as described previously 6,14,25 . Periodontal pocket depth, clinical attachment loss, and bleeding on probing were evaluated with a periodontal probe PCP UNC 15 (Hu-Friedy ® ). Severe periodontitis criteria was defined as the presence of deep periodontal pockets ≥6 mm with clinical attachment loss ≥5 mm and bleeding on probing in at least 8 sites distributed in different quadrants of the dentition 40,25 . Microarray, RT-qPCR and statistical analysis. Total RNA was extracted from peripheral blood mononuclear cells (PBMCs) since we focused on the circulating lymphocytes and monocytes gene expression signature. After extraction and purification of RNA samples, only those in the λ(260/280) and λ(260/230) ratios between 1.8 to 2.2 were used. Details regarding the aforementioned methodologies can be found in the Supplementary Materials and Methods. Microarray data (U133 Plus 2.0, Affymetrix Inc., Santa Clara, CA, USA) was generated from patients T2DMpoorly-DL-P (n = 5), T2DMwell-DL-P (n = 7), DL-P (n = 6), P (n = 6) and H (n = 6), after considering greater homogeneity regarding biochemical, lipid and clinical periodontal parameters 6 .
The raw.CEL files were preprocessed using the Robust Multichip Average (RMA) strategy, as implemented in the Affy package 47 , which performs background correction through a normal-exponential convolution model, quantile normalizes the probe intensities and summarizes them into probeset-level quantities using an additive model fit through the median-polish strategy. The log2-expression quantities resulted from the RMA method were further processed by the RankProd package 48 . Probesets that presented > FC log 1 2 and percentage of false prediction < . pfp ( ) 0 01 were called up-or downregulated, depending on the sign of the log-FC. Pairwise comparisons were made between the H (systemically and oral healthy subjects) and each of the T2DMpoorly-DL-P, T2DMwell-DL-P, DL-P and P groups (T2DMpoorly-DL-P versus H; T2DMwell-DL-P versus H; DL-P versus H and P versus H). The significant probesets from the aforementioned workflow were analysed with the Ingenuity Pathway Analysis software (IPA; Build 470319 M, Version 43605602, Qiagen, Redwood City, CA) as well as the Gene Set Enrichment Analysis software (GSEA; version 4.0.3 http://www.broad.mit.edu/gsea). Additional methods details from processing and analyses can be found at the Supplementary Materials and Methods section.
Independent validation of selected DEGs was performed by RT-qPCR analysis in a total of 150 patients (n = 30 each group, including the patients chosen for microarray analysis). Statistical analyses were performed in the GraphPad Prism software (version 5.0) using a significance level of 0.05 6 .

Data availability
The data that support the findings of this study are available from the corresponding author, RMSC, upon reasonable request.