Regulation of immune responses in primary biliary cholangitis: a transcriptomic analysis of peripheral immune cells

Abstract

Background aims: In patients with primary biliary cholangitis (PBC), the serum liver biochemistry measured during treatment with ursodeoxycholic acid-the UDCA response-accurately predicts long-term outcome. Molecular characterization of patients stratified by UDCA response can improve biological understanding of the high-risk disease, thereby helping to identify alternative approaches to disease-modifying therapy. In this study, we sought to characterize the immunobiology of the UDCA response using transcriptional profiling of peripheral blood mononuclear cell subsets.

Methods: We performed bulk RNA-sequencing of monocytes and TH1, TH17, TREG, and B cells isolated from the peripheral blood of 15 PBC patients with adequate UDCA response ("responders"), 16 PBC patients with inadequate UDCA response ("nonresponders"), and 15 matched controls. We used the Weighted Gene Co-expression Network Analysis to identify networks of co-expressed genes ("modules") associated with response status and the most highly connected genes ("hub genes") within them. Finally, we performed a Multi-Omics Factor Analysis of the Weighted Gene Co-expression Network Analysis modules to identify the principal axes of biological variation ("latent factors") across all peripheral blood mononuclear cell subsets.

Results: Using the Weighted Gene Co-expression Network Analysis, we identified modules associated with response and/or disease status (q<0.05) in each peripheral blood mononuclear cell subset. Hub genes and functional annotations suggested that monocytes are proinflammatory in nonresponders, but antiinflammatory in responders; TH1 and TH17 cells are activated in all PBC cases but better regulated in responders; and TREG cells are activated-but also kept in check-in responders. Using the Multi-Omics Factor Analysis, we found that antiinflammatory activity in monocytes, regulation of TH1 cells, and activation of TREG cells are interrelated and more prominent in responders.

Conclusions: We provide evidence that adaptive immune responses are better regulated in patients with PBC with adequate UDCA response.

FIGURE 1

Weighted gene co-expression network analysis (WGCNA) of RNA-sequencing data from CD14 cells. (A) Dendrogram showing gene co-expression networks (‘modules’). (B) Heatmap showing the strength of correlation between each module (rows) and each trait (columns). Asterisks indicate that the module is associated with the trait at q<0.05, where q is the false discovery rate-adjusted p-value. Rows outlined in black identify modules with at least moderate correlation (r≥0.5) with 1 or more traits and enriched with Hallmark gene sets at q<0.05. (C) Boxplots of the module eigengenes (ie, the first principal component of gene expression in the module; y-axis) in each group of participants (x-axis). Asterisks indicate the statistical difference between the module eigengene significance (based on the correlation between the module eigengene and trait of interest) and the trait comparison groups. Asterisks indicate the corresponding p values. (D) GeneMANIA protein-protein interaction plot showing protein-coding–hub genes in the salmon module for nonresponders versus controls. Protein-coding–hub genes are shown with cross-hatched circles of a uniform size, while those that were added as relevant genes by GeneMANIA are shown with solid circles whose size is proportional to the number of interactions they have. Lines correspond to the type of interactions. Purple: Genes known to be co-expressed in existing gene databases. Pink: Proteins known to be linked. Turquoise: Genes present in a shared annotated pathway. Blue: Genes expressed in the same tissue. Orange: Predicted functional relationships between genes. Green: Predicated genetic interactions. Olive: Shared protein domains. Abbreviations: WGCNA, Weighted gene co-expression network analysis.

FIGURE 2

Multi-omics latent factor analysis (MOFA) of modules identified by weighted gene co-expression network analysis (WGCNA) for all cell types. (A) Heatmap to show the fitted MOFA model, displaying the percentage of variance explained for each factor (rows) in each cell type. (B) Violin plots representing the distribution of the factor values for each patient group and each latent factor (1–7). ANOVA was used to identify latent factors associated with trait status; the p-value for each latent factor is shown. Abbreviations: MOFA, Multi-omics latent factor analysis; WGCNA, Weighted gene co-expression network analysis.

FIGURE 3

UDCA nonresponders show immunological differences compared to responders. Schematic illustrating the significant pathways correlated with nonresponders versus controls (A) and responders versus controls (B) overall highlighting the key findings of this study: (1) that monocytes are activated and proinflammatory in nonresponders; (2) T_H1 and T_H17 cells are activated in nonresponders and responders, but there is stronger regulation of both in the latter; (3) T_REG cells exhibit greater activation counterbalanced by greater regulation in responders; and (4) B cells are activated in responders.