Differential expression profile of genes involved in the immune response associated to progression of chronic Chagas disease

Abstract

Background: Patients with chronic Chagas disease present marked clinical and immunological heterogeneity. During the disease, multiple immune mechanisms are activated to fight the parasite. The purpose of this study was to investigate the expression patterns of genes involved in relevant immunological processes throughout the disease in patients with chronic Chagas disease.

Methodology/principal findings: High-throughput RT-qPCR with QuantStudio 12K Flex real-time PCR system was used to evaluate the expression of 106 immune-related genes in PBMC from a cohort of cardiac Chagas disease patients (CCC I), asymptomatic patients (IND) and healthy donors (HD) after being stimulated with T. cruzi soluble antigens. Principal component analysis (PCA), cluster analysis and volcano plots were used to identify differentially expressed genes. In addition, gene set enrichment analysis (GSEA) was employed to identify the enriched immunological pathways in which these genes are involved. PCA revealed the existence of a statistically divergent expression profile of the 36 genes correlated with PC1 between CCC I patients and HD (p < 0.0001). Differential gene expression analysis revealed upregulation of 41 genes (expression fold-change > 1.5) and downregulation of 14 genes (expression fold-change < 0.66) (p = 8.4x10-13 to p = 0.007) in CCC I patients versus HD. Furthermore, significant differences in the expression level of specific genes have been identified between CCC I and IND patients (8 up and 1 downregulated). GSEA showed that several upregulated genes in CCC I patients participate in immunological pathways such as antigen-dependent B cell activation, stress induction of HSP regulation, NO2-dependent IL12 pathway in NK cells, cytokines-inflammatory response and IL-10 anti-inflammatory signaling.

Conclusions: Cardiac Chagas disease patients show an antigen-specific differential gene expression profile in which several relevant immunological pathways seem to be activated. Assessment of gene expression profiles reveal unique insights into the immune response that occurs along chronic Chagas disease.

Fig 1. Principal component analysis (PCA) of gene expression profiles in cardiac Chagas disease and healthy subjects.

PCA was applied on NRQ values of expression level of 106 analyzed genes from cardiac Chagas disease patients (CCC I, orange triangles) and healthy donors (HD, blue circles). (A) PCA score plot of Principal Components 1 (PC1) and 2 (PC2) on the x and y-axis, respectively. (B) PCA score plot of Principal Components 1 (PC1) and 3 (PC3) on the x and y-axis, respectively. (C) 3D graphic representing the three Principal Components 1, 2 and 3 (PC1, PC2 and PC3). The proportion of total variance related to each principal component is given as a percentage and indicated on the axis next to the corresponding principal component.

Fig 2. Differential gene expression level in cardiac Chagas disease patients versus healthy donors.

(A) Volcano plot showing the differential expression level of the 106 analyzed genes in chronic Chagas disease patients with cardiac alterations (CCC I) (n = 18) and healthy donors (HD) (n = 30). The x-axis represents log₂ of expression fold-change between CCC I and HD (Log₂FC), where FC is calculated as the ratio between the average gene expression in the two groups (CCC I/HD). The y-axis corresponds to the statistical significance, expressed as the negative logarithm of the p-value (-Log₁₀ p-value). The purple horizontal lines indicate the cut-off for the statistical significance (corresponding to p = 0.05, p = 0.01, p = 0.001 and p = 0.0001). The grey vertical lines represent the log₂FC of -0.6 and 0.6 (corresponding to FC of 0.66 and 1.5, respectively) and log₂FC of −1 and 1 (corresponding to FC of 0.5 and 2, respectively) used as biological thresholds to identify differentially expressed genes. Negative values correspond to downregulated genes (blue dots) and positive values to upregulated genes (red dots) in CCC I compared to HD subjects. Black dots represent non-differentially expressed genes between CCC I and HD subjects. *Since the CCR5 gene did not show any detectable expression in the HD group, to apply the formulas, the value 0 was replaced by 10⁻⁶, resulting in a Log₂FC value of 19.35. (B) Expression level of 39 genes differentially expressed in cardiac Chagas disease patients versus healthy donors. Each plot bar corresponds to the gene referred at the bottom of the graph (x-axis). The y-axis represents log₂ of expression fold-change (Log₂FC) for each gene. Positive values (red bars) indicate upregulated genes and negative values (blue bars) indicate downregulated genes in CCC I compared to HD subjects. The grey horizontal lines highlight the genes with the largest differences in expression of at least more than four or less than a quarter-fold (Log₂FC > 2; Log₂FC < -2). *Since the CCR5 gene did not show any detectable expression in the HD group, to apply the formulas, the value 0 was replaced by 10⁻⁶, resulting in a Log₂FC value of 19.35.

Fig 3. Differential gene expression level in cardiac Chagas disease patients versus healthy donors.

(A) Volcano plot representing the differential gene expression level of the 106 analyzed genes in chronic Chagas disease patients with indeterminate form (IND) (n = 39) and patients with cardiac alterations (CCC I) (n = 18). The x-axis represents log₂ of expression fold-change between CCC I and IND (Log₂FC), where FC is calculated as the ratio between the average gene expression in two groups (CCC I/IND). The y-axis corresponds to the statistical significance, expressed as the negative logarithm of the p-value (-Log₁₀ p-value). The purple horizontal line indicates the cut-off for the statistical significance (corresponding to p = 0.05). The grey vertical lines represent the log₂FC of −0.6 and 0.6 (corresponding to FC of 0.66 and 1.5, respectively) used as biological thresholds to identify differentially expressed genes. Negative values correspond to downregulated genes (blue dots) and positive values to upregulated genes (red dots) in CCC I compared to IND patients. Black dots represent non-differentially expressed genes between CCC I and IND patients. (B) Relative expression level of 9 genes differentially expressed in cardiac versus indeterminate Chagas disease patients. Each plot bar corresponds to the gene referred at the bottom of the graph (x-axis). The y-axis represents the log₂ of expression fold-change (Log₂FC) for each gene. Positive values (red bars) indicate upregulated genes and negative values (blue bars) indicate downregulated genes in CCC I compared to IND patients. The grey horizontal line indicates the cut-off to highlight the genes with the greatest differences in expression (Log₂FC > 1).

Fig 4. Comparison of gene expression level in Chagas disease patients (cardiac and indeterminate) and healthy donors.

Expression level of CCR5, IFNG, IL1B, IL5, IL5RA, IL6, IL10, IL10RB and GZMM genes measured as mean Normalized Relative Quantities (NRQ) in patients with chronic Chagas disease with indeterminate form (IND), patients with cardiac alterations (CCC I) and healthy donors (HD). Statistically significant differences determined by two-tail Mann-Whitney test or two-tailed unpaired t-test, as appropriate, are indicated (* p < 0.05, ** p < 0.01, *** p < 0.001 and **** p < 0.0001).

Fig 5. Gene set enrichment analysis (GSEA) in cardiac Chagas disease patients versus healthy donors.

GSEA heatmaps of representative gene sets from (A) Antigen dependent B cell activation pathway (BIOCARTA_ASBCELL_PATHWAY), (B) Stress induction of HSP regulation (BIOCARTA_HSP27_PATHWAY), (C) Cytokines and inflammatory response (BIOCARTA_INFLAM_PATHWAY), (D) NO2-dependent IL12 pathway in NK cells (BIOCARTA_NO2IL12_PATHWAY), and (E) IL-10 anti-inflammatory signaling pathway (BIOCARTA_IL10_PATHWAY) enriched in CCC I patients versus HD subjects. The values of the expression level of each gene are represented as colors, ranging from dark red to dark blue, based on the highest and lowest normalized relative quantity (NRQ) value of each gene, respectively. Parameters set for GSEA were: Molecular Signatures Database (MSigDB) BioCarta gene set collection, permutations = 100,000, permutation type: phenotype, enrichment statistic: weighted, metric for ranking genes: t-test, max size: 500, min size: 3.