Systems Immunology of Diabetes-Tuberculosis Comorbidity Reveals Signatures of Disease Complications

Abstract

Comorbid diabetes mellitus (DM) increases tuberculosis (TB) risk and adverse outcomes but the pathological interactions between DM and TB remain incompletely understood. We performed an integrative analysis of whole blood gene expression and plasma analytes, comparing South Indian TB patients with and without DM to diabetic and non-diabetic controls without TB. Luminex assay of plasma cytokines and growth factors delineated a distinct biosignature in comorbid TBDM in this cohort. Transcriptional profiling revealed elements in common with published TB signatures from cohorts that excluded DM. Neutrophil count correlated with the molecular degree of perturbation, especially in TBDM patients. Body mass index and HDL cholesterol were negatively correlated with molecular degree of perturbation. Diabetic complication pathways including several pathways linked to epigenetic reprogramming were activated in TBDM above levels observed with DM alone. Our data provide a rationale for trials of host-directed therapies in TBDM, targeting neutrophilic inflammation and diabetic complication pathways to address the greater morbidity and mortality associated with this increasingly prevalent dual burden of communicable and non-communicable diseases.

Conflict of interest statement

The authors declare that they have no competing interests.

Figures

Figure 1

Factors associated with TB and/or DM in the cohort. (A) Data represent median and interquartile ranges. The Kurskal Wallis test with Dunn’s multiple comparisons ad hoc test were employed to compare the values detected between the study subgroups. P-values for the Kruskal Wallis tests are shown in Appendix File S1. P-values were adjusted for multiple comparisons using Holm-Bonferroni’s approach as described in methods. Only comparisons with significant P-values after adjustments are displayed (*P < 0.05, **P < 0.01, ***P < 0.001). (B) Bayesian network with bootstrap (100x) was used to illustrate the statistically significant associations between the parameters and the presence of TB and/or DM in the study population. Lines represent direct associations. Associations that remained statistically significant on >30 times out of 100 bootstraps are plotted. Numbers of times each association persisted during bootstrap are shown. Bold lines highlight the strongest associations, which persisted more than 60 times in the bootstrap. Vitamin D is not presented in the picture, as it would be solitary node with no arcs attached.

Figure 2

Cytokine profiling in plasma of patients with TB and/or DM. (A) Cytokine concentrations were z-score normalized across all subjects. Each column represents one patient. The cytokine names are shown to the right of the heat map. Cytokine profiles were ordered by hierarchical clustering (Euclidean distance and clustered with ward method). The condition tree at the top shows 3 main clusters. The number of subjects from each clinical phenotype (class) on each cluster is indicated. (B) Differences in cytokine levels for each disease compared to healthy subjects. Differences which did not reach statistical significance (Adjusted P < 0.05, fold-change >1.4) are represented as grey bars.

Figure 3

Transcriptomic changes with TB and/or diabetes compared to healthy subjects. (A) Differentially expressed genes (DEGs) in patients with TB and/or DM compared to healthy subjects (adjusted P < 0.05, fold-change >1.4). The Venn diagram shows the number of DEGs in common to two or more clinical phenotypes or unique to only one disease. (B) Expression patterns of the DEGs from (A). Cohort subgroups are shown by the colored bars. Each column represents one subject. The genes (rows) were normalized by z-score across all samples.

Figure 4

Blood transcriptome of patients with dual TB and DM burden compared to TB or DM alone. (A) Differentially expressed genes between patients with distinct clinical phenotypes (class). The comparison between two classes is indicated at the top of each bar. Red and blue bars represent the number of DEGs which are up- or down-regulated, respectively in the first class of the comparison. The purple bar shows the number of DEGs in the comparison TBDM vs. DM (“TBDM total”). The number of DEGs which are in common to two comparisons are shown inside green areas connecting two horizontal bars (“TB shared”). Genes which are uniquely differentially expressed in TBDM vs. DM is shown by the brown bar (“TBDM unique”). (B) Top transcriptional regulators in “TB shared” DEGs. Ingenuity’s Upstream Regulator Analysis was performed using “TB shared” (green bars), “TBDM unique” (brown bars), and “TBDM total” DEGs (purple bars). The x-axis indicates the significance of the enrichment for the upstream regulators on the left. (C) Top transcriptional regulators in “TBDM unique” DEGs identified by the same approach as in (B).

Figure 5

Factors associated with increased molecular degree of perturbation. (A) The molecular degree of perturbation (MDP) relative to healthy controls was calculated as described in Methods. Left panels show histograms of individuals in each subgroup. Right panel shows individual values with median and IQR per group. Data were compared using one-way ANOVA with Tukey’s multiple comparisons test, with single pooled variance. (B) MDP values were tested for correlations with the indicated parameters in each study subgroup using Spearman correlation ranks. Statistically significant correlations, after Holm-Bonferroni’s adjustment for multiple comparisons, are highlighted (red, positive correlation; blue, negative correlation; black, non-significant correlation). (C) Left panel shows correlation plots between differential WBC counts (neutrophils, monocytes and lymphocytes) and the MDP in all study groups. Dotted lines represent median values for MDP or leukocyte counts/frequency for the entire study population. Shaded areas indicate the individuals with the highest MDP values and the lowest lymphocyte counts/frequency or the highest neutrophil counts/frequency relative to the study population. Right panel shows comparisons of the frequencies in each study subgroup of individuals from the shaded areas on the left. Frequency comparisons were performed using the chi-squared test. (D) Left panel, MDP score was tested for correlation with the radiographic severity score (Spearman rank test). Right panel shows comparisons of MDP values in patients diverging in terms of TB disease distribution (unilateral vs. bilateral lung lesions) or presence of cavitary lesions. Within each study subgroup, values were compared using the Mann-Whitney U test. In all comparisons: *P < 0.05, **P < 0.01, ***P < 0.001.

Figure 6

A blood transcriptional signature of active TB classifies most patients with or without diabetes. We applied a previously defined 393-gene signature of active TB to the 387 genes represented in the microarrays used for the present study. Expression profile of each patient and health subject (column) was ordered by hierarchical clustering (Spearman correlation with average linkage). The condition tree at the top shows 2 main clusters. The number of individuals from all each clinical subgroup (class) who segregated into each cluster is indicated at the top of the figure. Colored blocks at each profile base represent the different classes, molecular perturbation status, MDP score and gender.

Figure 7

Increased activity of pathways associated with epigenetic regulation in diabetic TB patients. (A) Single sample Gene Set Enrichment Analysis (ssGSEA) was performed for each patient using the genes ranked by their MDP score and Reactome pathways as gene sets (P < 0.01, 1,000 permutations). Colors represent increase (red) or decrease (blue) in pathway expression activity vs. healthy controls as given by the normalized enrichment score (NES). Bars at the bottom indicate patients which are perturbed (darker color) or not-perturbed (light color). (B) Mean NES (values inside the circles) of the patients and pathways for each functional group. The size of the circles is proportional to the mean NES and color indicates positive (red) or negative (blue). (C) Pathways associated with “DNA methylation”. Spearman correlation between the expression activity of “DNA methylation” pathway and all other Reactome pathways across all patients. The network shows the pathways with correlation >0.6 (red circles) or <−0.6 (blue circles). The thickness of the edges is proportional to the number of genes in common between two given pathways (nodes).

Figure 8

Integrative analysis of cytokines, gene expression and pathway activity in TBDM. (A) Integration of cytokine levels and microarray gene expression. Sparse group partial least square methods were used to identify the genes whose expression best correlates with most of the cytokines. Colors represent the correlation values between the expression of genes in blood (columns) and the levels of cytokines in plasma (rows). (B) Network containing some of the 114 genes from (A) GeneMania program was used to define the interactions (edges) between genes (nodes), and Gephi program to visualize the network. Colors represent is positive (red) or negative (blue) correlation of gene expression with the cytokines. (C) Integration of gene expression and pathway activity. Circos plot shows 8 selected pathways and 114 genes from (A). Each pathway was assigned to one color (inner circle). The lines connecting pathways and genes represent a high correlation (spearman correlation above 0.6) between the activity of the pathway (defined by the single-sample enrichment score) and the gene expression. Heat maps shown in the outer circles represent the activity of pathways (mean Normalized Enrichment Score) and genes (mean log2 fold-change in patients compared to healthy subjects). Colors indicate increase (red) or decrease (blue) in pathway activity or gene expression fold-change.