Identification of Marker Genes and Pathways in Patients with Primary Biliary Cholangitis

Abstract

Primary biliary cholangitis (PBC) is an autoimmune liver disease characterized by cholestasis and cirrhosis, and in which hepatic failure may occur. This study explores the changes in the gene expression profiles of liver tissues during the pathogenesis of PBC. Array dataset GSE79850 was downloaded from the Gene Expression Omnibus database. GeneSpring software was used to analyze differentially expressed genes (DEGs) in liver tissues from PBC patients compared with those from controls. Gene ontology (GO) annotation, the Kyoto Encyclopedia of Genes and Genomes (KEGG), and Reactome pathway enrichment analyses were performed by using Database for Annotation, Visualization and Integrated Discovery (DAVID) software. Cytoscape software was used to construct a protein–protein interaction (PPI) network. Plug-ins Molecular Complex Detection and iRegulon were used for clustering analysis and transcription factors related to key genes with PBC. A total of 77 DEGs, including 47 up- and 30 downregulated genes, were identified. The PPI network was established with 74 nodes and 356 protein pairs. The C-C motif chemokine ligand 5 (CCL5), interleukin 7 receptor (IL7R), and TNF receptor superfamily member 1A (TNFRSF1A) were identified as hub genes in the PPI network and may, therefore, be marker genes for PBC. Further, the upregulated genes CCL5 and IL7R, and downregulated TNFRSF1A were included in immune system processes as a GO term in the category Biological Processes. In conclusion, CCL5, IL7R, TNFRSF1A, and the immune response pathway may have crucial roles in PBC. These genes and pathways may be potential targets for treating PBC.

1. Introduction

Primary Biliary Cholangitis (PBC) is an intrahepatic cholestasis disease involving the interlobular bile duct in the liver. The etiology and pathogenesis of PBC have not been fully clarified. Genetic susceptibility, environmental factors, and immune tolerance deficits are related to the pathogenesis of PBC (Chen et al., 2016). PBC has a hidden onset, and symptoms of patients with PBC are often atypical in the early stage and tend to be ignored and easily misdiagnosed. Once diagnosed, the illness tends to be severe, progressing to a decompensation period of cirrhosis, with refractory ascites, upper gastrointestinal bleeding, and hypersplenism. Anti-mitochondrial antibody (AMA) is a specific diagnostic marker of PBC and is used for the early diagnosis of PBC. The positive rate of the AMA-M2 subtype is 90%–95%. However, many studies have suggested that AMA is not relevant to clinical disease progression and disease classification, and AMA cannot be cleared from the body after treatment or liver transplantation (Luettig et al., 1998; Miyakawa et al., 2001). Because of intrahepatic cholestasis, the most significant changes in biochemical indicators in PBC patient serum are those in alkaline phosphatase (ALP) and gamma-glutamine peptide. A total of 96% of PBC patients demonstrated ALP levels that are typically 2–10 times higher than healthy levels, and this is seen in early and asymptomatic patients but these indicators are not specific.

The diagnosis of PBC needs to meet two of the following three criteria (European Association for the Study of the Liver, 2009; Lindor, 2009; Bowlus and Gershwin, 2014; Zweers et al., 2017): (1) AMA titer >1:40; (2) ALP level is 1.5 times higher than the normal upper limit and lasts for more than 24 weeks; and (3) liver biopsy shows nonsuppurative cholangitis and interlobular bile duct damage. However, patients with AMA negativity have difficulty in early diagnosis and differentiation from other autoimmune liver diseases and to make a definite diagnosis requires an invasive liver biopsy that not only increases the economic burden but also results in psychological and physical harm to PBC patients, and often misses the best timing for treatment. Moreover, ursodeoxycholic acid (UDCA), the only drug approved by the FDA for the treatment of PBC, failed to treat some patients with poor prognosis and such patients later progressed to liver failure. Therefore, early diagnosis and evaluation of PBC are urgently required.

PBC is a complex classic autoimmune disease with obvious genetic components and strong environmental effects. Although the etiology of the disease is unclear, the role of genetic factors is widely known and when combined with exposure to environmental factors, it may lead to the collapse of immune tolerance. It is now acknowledged that patients with PBC lose their tolerance to mitochondrial antigens, and in particular, the E2 component of the pyruvate dehydrogenase complex results in biliary epithelial cell invasion, ultimately leading to bile duct injury, cirrhosis, and liver failure. Genome-wide association studies (GWAS) on PBC patients have been extensive, and many studies have found that there is a significant association between different populations and the onset of PBC. A GWAS study of human leukocyte antigen (HLA) showed that the HLA-DQA1 and HLA-DQB1 loci are strongly correlated with PBC disease susceptibility and, in particular, the HLA-DRB1T, HLA-DQA1, and HLA-DQB1 loci.

However, 80%–90% of PBC patients do not carry these common HLA susceptibility alleles, so HLA cannot explain the entire genetic susceptibility underlying PBC, which indicates that other genes are related to the onset of PBC (Joshita et al., 2017). To date, GWAS have demonstrated that there are dozens of non-HLA sites in the entire genome that are significantly associated with PBC susceptibility (Hirschfield et al., 2009, 2010; Liu et al., 2010, 2012; Juran et al., 2012; Cordell et al., 2015). Moreover, many immune pathways are involved in PBC. For example, IL12-JAK-STAT4 signaling promotes Th1 T cell polarization and is associated with a risk to PBC patients.

In this study, for the purpose of screening the gene expression profiles of the liver tissue of PBC patients, we downloaded the GSE79850 dataset (Hardie et al., 2016) from the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo) (Barrett et al., 2013), and GeneSpring software was used to analyze differentially expressed genes (DEGs). Moreover, gene ontology (GO) and pathway enrichment analyses were carried out. The related key genes and pathways may be potential diagnostic and therapeutic targets for PBC.

2. Methods

2.1. Gene Expression Omnibus data processing

The GEO database is an international public resource library for receiving, registering, archiving, searching, and downloading microarray chip, next-generation sequencing and other forms of high-throughput functional genomics data. We downloaded the GSE79850 array dataset based on the 770 gene Pan-Cancer Immunity panel platform (Nanostring^®; nCounter, Seattle, WA) from the GEO database. The dataset contains 24 liver tissue samples, including 9 high-risk PBC samples, 7 low-risk PBC samples, and 8 control samples, and we combined the high-risk and the low-risk into one group called the PBC group.

2.2. Identification of differentially expressed genes in primary biliary cholangitis

GeneSpring software (version 11.5; Agilent) was used to analyze the array database. The chip model could not be identified, so we used “created Technology” in GeneSpring to create a new platform (nanoString GPL19965 platform) and then imported matrix and comment files. We applied “Experiment grouping” to divide the data into “PBC” and “control” groups, respectively, in accordance with the sample information. The principal component analysis was used to control probe quality. “Filter ProbeSets by Expression” was set with intensity values less than 20% in GeneSping for filtering out the quality of the probes. An unpaired t-test was used to identify DEGs with >2.0-fold change in expression, and p < 0.05 was defined as statistically significant. Mev software was used to construct a heat map.

2.3. Gene ontology and pathway enrichment of differentially expressed gene analyses

GO was used to analyze the biological function of the DEGs, including three aspects of biology: biological processes (BP), cellular components (CC), and molecular function (MF) (Gene Ontology Consortium, 2006). The Kyoto Encyclopedia of Genes and Genomes (KEGG) is one of the most commonly used biological information databases. Its main function is to provide information at a molecular level, especially for high-throughput sequencing technology used in improving the function of biological systems (Kanehisa and Goto, 2000). The Database for Annotation, Visualization and Integrated Discovery (DAVID) 6.8 (https://david.ncifcrf.gov) was used to analyze GO terms, KEGG, and Reactome pathways. p < 0.05 was defined as statistically significant.

2.4. Analysis of protein–protein interaction network

We used the Search Tool for the Retrieval of Interacting Genes (STRING, version 10.5) database, an online tool, to evaluate the interactive relationships among DEGs (Szklarczyk et al., 2011, 2015). All of the DEGs in this study were uploaded to STRING, and the threshold value was set with a combined score of >0.4. To construct protein–protein interaction (PPI) networks, Cytoscape software was used (Shannon et al., 2003). Hub genes were selected by using the cytoHubba plug-in with degrees of interaction ≥10. Molecular Complex Detection (MCODE) in Cytoscape was used to filter the clusters of the PPI network with scores >3.

GenCLiP 2.0 (http://ci.smu.edu.cn/GenCLiP2/analysis.php), the online tool for literature mining (Wang et al., 2014), was used to analyze the hub genes obtained from the PPI network. We used the “Literature Mining Gene Networks” option to construct the co-citation network of hub genes. Subsequently, the related biological function of these hub genes was analyzed by using the “Gene Cluster with Literature Profiles” option with p ≤ 1 × 10⁻¹⁸ and hit ≥18.

2.5. Analysis of marker genes and transcription factors related to primary biliary cholangitis

The comparative toxicogenomics database (CTD) is an online tool for formalizing, harmonizing, and centralizing gene and protein data for different species (Davis et al., 2015). We evaluated whether any hub genes among the DEGs had been confirmed as marker genes of PBC in the CTD database in previous studies. “PBC” was used as the keyword. Then, the plug-in iRegulon in Cytoscape (Janky et al., 2014) was used to calculate transcription factors related to our marker genes. Marker genes and their corresponding transcription factor pairs with Normalized Enrichment Scores (NES) >5 were selected, and the regulatory transcription factor network was the output.

3. Results

3.1. Identification of differentially expressed genes

The GeneSpring analysis of GSE79850 in liver tissue from PBC patients compared with control patients yielded 77 DEGs in total (p < 0.05 and >2.0-fold change in expression). This included 47 upregulated genes and 30 downregulated genes (as shown in Fig. 1).

FIG. 1.

Heat map of differentially expressed genes. Light grey color represents a lower expression level, and dark grey color represents a higher expression level. Black color indicates no differential expression.

3.2. Gene ontology term enrichment analyses

We uploaded upregulated or downregulated DEGs to DAVID to analyze the GO terms, KEGG, and Reactome pathways. The GO analysis results demonstrated that all the DEGs were significantly enriched in BPs. The upregulated DEGs were enriched in genes involved in the immune response, response to cytokines, and cellular responses to cytokine stimuli (Table 1); the downregulated DEGs were enriched in genes involved in the defense response, inflammatory response, and response to cytokines (Table 1). The GO-CC analysis also indicated that the upregulated DEGs were significantly enriched in major histocompatibility complex (MHC) class II protein complex, MHC protein complex, and cell membrane components; downregulated DEGs were enriched in extracellular space, cell surface, and intrinsic components of the plasma membrane (Table 1). In addition, for MF, the upregulated DEGs were enriched in antigen binding, MHC class II receptor activity, and signal transducer activity (Table 1); the downregulated DEGs were enriched in signal transducer activity, macromolecular complex binding, and identical protein binding (Table 1).

Table 1.

Gene Ontology Analysis of Differentially Expressed Genes (p < 0.05)

GO ID	Term	No. of genes	p
Upregulated genes
GO-BP terms
GO:0006955	Immune response	35	7.94E-26
GO:0034097	Response to cytokine	26	2.08E-21
GO:0071345	Cellular response to cytokine stimulus	25	2.45E-21
GO:0019221	Cytokine-mediated signaling pathway	23	7.36E-21
GO:0070887	Cellular response to chemical stimulus	37	1.08E-20
GO-CC terms
GO:0042613	MHC class II protein complex	7	3.90E-11
GO:0042611	MHC protein complex	7	3.05E-10
GO:0098552	Side of membrane	11	6.28E-07
GO:0098553	Lumenal side of endoplasmic reticulum membrane	5	1.65E-06
GO:0071556	Integral component of lumenal side of endoplasmic reticulum membrane	5	1.65E-06
GO-MF terms
GO:0003823	Antigen binding	10	3.90E-10
GO:0032395	MHC class II receptor activity	5	8.31E-08
GO:0004871	Signal transducer activity	20	1.41E-07
GO:0038023	Signaling receptor activity	16	6.05E-06
GO:0004872	Receptor activity	17	1.03E-05
Downregulated genes
GO-BP terms
GO:0006952	Defense response	27	2.05E-24
GO:0006955	Immune response	26	2.30E-22
GO:0034097	Response to cytokine	21	2.34E-20
GO:0071345	Cellular response to cytokine stimulus	19	3.46E-18
GO:0009607	Response to biotic stimulus	20	4.61E-18
GO-CC terms
GO:0005615	Extracellular space	4	2.61E-05
GO:0009986	Cell surface	12	3.84E-04
GO:0031226	Intrinsic component of plasma membrane	8	2.65E-03
GO:0005887	Integral component of plasma membrane	10	7.88E-03
GO:0044440	Endosomal part	9	8.45E-03
GO-MF terms
GO:0004871	Signal transducer activity	11	3.40E-04
GO:0044877	Macromolecular complex binding	9	9.88E-04
GO:0042802	Identical protein binding	9	1.28E-03
GO:0005102	Receptor binding	9	2.06E-03
GO:0004888	Transmembrane signaling receptor activity	8	4.77E-03

BP, biological process; CC, cellular component; GO, gene ontology; MF, molecular function; MHC, major histocompatibility complex.

3.3. Pathway analysis

The KEGG pathways of upregulated DEGs were significantly enriched in the type A influenza and herpes simplex infection (A in Table 2); correspondingly, downregulated DEGs were significantly enriched in the tumor necrosis factor (TNF) signaling pathway and the nuclear factor-kappa B (NF-κB) signaling pathway (A in Table 2). In the Reactome pathway of the DEGs, the upregulated genes were significantly enriched in interferon gamma signaling and MHC class II antigen presentation (B in Table 2); whereas the downregulated genes were enriched in TAK1 activates NF-κB by phosphorylation and activation of I-kappa-B-kinase (IKK) complex and the nucleotide-binding oligomerization domain (NOD)1/2 signaling pathway (B in Table 2).

Table 2.

Kyoto Encyclopedia of Genes and Genomes and Reactome Pathway Enrichment Analyses for Differentially Expressed Genes (p < 0.05)

ID	Description	No. of genes	p
(A) KEGG pathway
Upregulated
05164	Influenza A	14	2.20E-12
05168	Herpes simplex infection	14	4.21E-12
05145	Toxoplasmosis	12	1.32E-11
05416	Viral myocarditis	9	4.13E-10
05321	Inflammatory bowel disease	9	1.08E-09
05140	Leishmaniasis	9	2.53E-09
04672	Intestinal immune network for IgA production	8	3.74E-09
05310	Asthma	7	8.34E-09
05152	Tuberculosis	11	1.81E-08
05330	Allograft rejection	7	3.18E-08
05320	Autoimmune thyroid disease	7	2.63E-07
05150	Staphylococcus aureus infection	7	3.31E-07
05166	HTLV-I infection	11	5.88E-07
04514	Cell adhesion molecules	9	6.09E-07
05332	Graft-vs.-host disease	6	7.08E-07
04940	Type I diabetes mellitus	6	2.45E-06
04612	Antigen processing and presentation	7	2.57E-06
05323	Rheumatoid arthritis	7	6.08E-06
04620	Toll-like receptor signaling pathway	7	1.79E-05
04060	Cytokine–cytokine receptor interaction	9	2.23E-05
04151	PI3K-Akt signaling pathway	10	6.23E-05
05162	Measles	7	6.47E-05
05322	Systemic lupus erythematosus	7	6.74E-05
05215	Prostate cancer	6	9.45E-05
04621	NOD-like receptor signaling pathway	5	1.95E-04
04062	Chemokine signaling pathway	7	4.07E-04
05160	Hepatitis C	6	6.51E-04
Downregulated
04668	TNF signaling pathway	8	6.27E-08
04064	NF-κB signaling pathway	6	1.23E-05
05164	Influenza A	7	2.78E-05
05142	Chagas disease (American trypanosomiasis)	6	2.93E-05
05160	Hepatitis C	6	9.55E-05
05152	Tuberculosis	6	3.65E-04
05168	Herpes simplex infection	6	4.26E-04
04621	NOD-like receptor signaling pathway	4	9.70E-04
(B) Reactome pathway
Up-regulated
877300	Interferon gamma signaling	8	7.57E-08
2132295	MHC class II antigen presentation	7	1.18E-05
202430	Translocation of ZAP-70 to immunological synapse	4	1.07E-05
202427	Phosphorylation of CD3 and TCR zeta chains	4	1.58E-04
389948	PD-1 signaling	4	1.78E-04
202433	Generation of second messenger molecules	4	4.75E-04
202424	Downstream TCR signaling	5	8.82E-04
933543	NF-κB activation through FADD/RIP-1 pathway mediated by caspase-8 and -10	3	1.16E-03
5673001	RAF/MAP kinase cascade	5	1.48E-03
Downregulated
445989	TAK1 activates NF-κB by phosphorylation and activation of IKKs complex	5	1.17E-06
168638	NOD1/2 signaling pathway	4	1.24E-04
446652	Interleukin-1 family signaling	4	1.24E-04
2559582	Senescence-associated secretory phenotype	5	2.12E-04
140534	via Death receptors in the presence of ligand	3	7.48E-04
937072	TRAF6-mediated induction of TAK1 complex within TLR4 complex	3	9.66E-04

FADD, Fas-associated via death domain; HTLV-I, human T cell leukemia virus-I; IKK, I-kappa-B-kinase; KEGG, Kyoto Encyclopedia of Genes and Genomes; NF-κB, nuclear factor kappa B; MAP, mitogen-activated protein; NOD, nucleotide-binding oligomerization domain; PD-1, programmed death-1; PI3K, phosphoinositide-3-kinase, AKT serine/threonine kinase; RAF, Raf-1 proto-oncogene, serine/threonine kinase; RIP, receptor interacting serine/threonine; TAK1, TGF-beta-activated kinase; TCR, T cell receptor; TLR4, toll-like receptor 4; TRAF6, TNF receptor-associated factor 6; ZAP, zeta chain of T cell receptor-associated protein kinase.

3.4. Protein–protein interaction network analysis

Based on the analysis of the STRING database, data with a PPI score of >0.4 were selected. We obtained 74 nodes and 356 protein pairs (Fig. 2). The top 10 hub genes with higher degrees included intercellular adhesion molecule 1 (degree = 29), signal transducer and activator of transcription 3 (degree = 27), mitogen-activated protein kinase 3 (degree = 26), C-C motif chemokine ligand 2 (degree = 24), Caspase 8 (degree = 24), CD40 Molecule (degree = 23), C-C motif chemokine ligand 5 (CCL5) (degree = 23), Fos proto-oncogene, AP-1 transcription factor subunit (degree = 23), Janus Kinase 2 (degree = 23), and interferon regulatory factor 1 (degree = 22).

FIG. 2.

Protein–protein interaction network of differentially expressed genes. Dark grey nodes represent upregulated genes, and light grey nodes represent downregulated genes.

All the selected nodes and edges were analyzed by using the MCODE module. We obtained two clusters with MCODE scores >3. Two significant modules called “cluster 1” and “cluster 2” were selected (Fig. 3). Genes involved in the two clusters were functionally annotated (Fig. 3a). Enrichment analysis demonstrated that the genes in cluster 1 were mostly enriched in the TNF signaling pathway, Influenza A, Chagas disease (American trypanosomiasis), and Toxoplasmosis; whereas the genes in cluster 2 were enriched in herpes simplex virus, antigen processing and presentation, asthma, and graft-versus-host disease.

FIG. 3.

Analysis of 74 nodes and 356 edges by plug-ins MCODE. (a) Cluster1 and (b) cluster 2 were the top two significant clusters. MCODE, Molecular Complex Detection.

Literature mining of 29 genes with a degree >10 in the PPI network was performed. A total of 27 previously reported genes were included to construct the co-citation network (as shown in Fig. 4a) and previous studies regarding the top 10 of these genes are shown in Table 3. The enrichment analysis showed that these 27 genes were enriched in terms of cell activation, the innate immune response, and the immune system (Fig. 4b).

FIG. 4.

(a) Co-citation network of 27 genes with a higher degree of expression according to previous studies. The numbers on the linking lines indicate the number of studies co-cited. (b) Significant enrichment results of these 27 genes from previous studies.

Table 3.

Hub Genes Identified by the Present Study in the Literature

Gene	Co-genes	Co-cite	Total
MAPK3	16	1729	23983
STAT3	16	2289	13379
FOS	13	775	19906
CASP8	12	1165	9235
CCL5	11	2622	7812
TNFSF10	10	1876	7696
CD40	10	894	9389
ICAM1	10	2047	21539
JAK2	10	1417	6440
NFKBIA	19	504	5534

3.5. Marker genes and transcription factors related to primary biliary cholangitis analyses

Three previously reported marker genes associated with PBC that are listed in the CTD, including CCL5, interleukin 7 receptor (IL7R), and TNF receptor superfamily member 1A (TNFRSF1A), were identified among the DEGs in this study. The relative expression levels of these three genes in the liver tissues of PBC and control patients are shown in Figure 5a. The transcription regulatory network of these marker genes is shown in Figure 5b. There are five transcription factors with an NES >5, including interferon regulatory factor 8 (IRF8, NES = 7.701), B Cell CLL/Lymphoma 3 (BCL3, NES = 7.061), TATA-box binding protein (TBP, NES = 6.878), ETS transcription factor (ELK1, NES = 6.625), and Basonuclin 1 (BNC1, NES = 5.744).

FIG. 5.

(a) Expression levels of three marker genes in liver biopsies from PBC and control that had been listed as liver cirrhosis, biliary-associated marker genes in the comparative toxicogenomics database. (b) Regulatory transcription factor network of the marker genes. PBC, primary biliary cholangitis.

4. Discussion

In this study, a total of 77 DEGs (47 upregulated and 33 downregulated genes) were identified in the liver tissues of PBC patients compared with controls within the GSE79850 array dataset. The analysis showed that CCL5 had a high degree of interaction in the PPI network and was a marker gene in PBC. Moreover, upregulated DEGs were significantly enriched in the Reactome pathways of interferon gamma signaling, and in various associated terms of the GO-BP category with more significant p-values. In addition, CCL5 was a specific type of upregulated DEG that was enriched in the immune system processes of the GO-BP category. CCL5 has chemotactic effects on T cells, eosinophils, and basophils, and it plays an active role in the recruitment of leukocytes to inflammatory sites. It has been shown that CCL5 is associated with many autoimmune diseases, such as Sjögren syndrome (Hjelmervik et al., 2005), and rheumatoid arthritis (Shadidi et al., 2003). Nischalke et al. (2004) detected upregulation of CCL5 in liver tissue, whereas Zhang et al. (2014) performed a quantitative real-time (qRT)-PCR analysis of CCL5 expression and confirmed a downregulation of expression in peripheral blood mononuclear cells. In this study, the relative expression of CCL5 was upregulated in the liver tissues of patients with PBC compared with control patients. This inconsistency between peripheral blood and liver tissue possibly indicates that CCL5 chemokines accumulate during the process of inflammation in the liver. It also provides a basis for further study of the immune mechanism of PBC.

Moreover, TNFRSF1A also had a high degree of interaction in the PPI network and was a marker gene associated with PBC in this study. It plays a part in cell survival, apoptosis, and inflammatory processes. It also plays a role in autoimmune diseases, and mutations in the TNFRSF1A gene may be associated with multiple sclerosis (Kulakova et al., 2016). Single-nucleotide polymorphism of TNFRSF1A (rs1800693) is associated with PBC (Mells et al., 2011; Dong et al., 2015) in both the United Kingdom and Chinese populations. In this study, TNFRSF1A is involved in the TNF signaling pathway, in which downregulated genes were significantly enriched in the KEGG pathways. This gene encodes one of two receptors for TNFα, and compared with wild-type mice, the degree of liver fibrosis after administering a potent hepatotoxin was attenuated in c^−/− mice (Kitamura et al., 2002). Therefore, TNFRSF1A may be a key gene associated with PBC.

IL7R was identified as another PBC marker gene in this study. IL7R (rs860413) has been confirmed as a candidate gene by the U.K. PBC Consortium (Mells et al., 2011). The rs6890853 single-nucleotide polymorphism located in the IL7R gene showed significant associations with PBC in the Japanese population (Nakamura et al., 2012). However, no association was observed for IL7R (rs6897932) in a Han Chinese PBC cohort (Dong et al., 2015). IL7R mutation is a type of heterogeneous genetic disorder that causes Omenn syndrome, characterized by reduced T cell mutations, immune disorders, and autoimmune and graft-versus-host disease (Liston et al., 2008). Therefore, IL7R is another key gene associated with PBC and a molecular mechanism analysis is a starting point to study the pathogenesis of PBC.

In addition, the first four pathways of upregulated DEGs identified in this study were significantly enriched in viral and parasitic infection pathways in the KEGG pathways. A molecular simulation mechanism was demonstrated, which is associated with a pathogen-triggered autoimmune reaction mechanism. The mechanism shows that the pathogen epitope (i.e., peptides, carbohydrates or DNA sequences) shares a homology with the pathogen itself to varying degrees, and in response to pathogens and their hosts results in a mixed antibody- and cell-mediated immune response (Van de Water et al., 2001). To date, specific microbes that have been reported as potentially related to PBC have been mainly bacteria, but they also include viruses, parasites, and fungi (Selmi et al., 2010). Therefore, in this study, the upregulated DEGs may indicate that those genes associated with PBC were cross-reactive with microbial infection.

CCL5, IL7R, and TNFRSF1A were all enriched in the term “immune response” in the GO-BP category in this study. Therefore, we speculate that these three genes may be involved in PBC by way of the immune response pathway. The main immunological topics focus on the loss of immune tolerance, resulting in T cell-mediated injury, an imbalance in the regulation of immune cells, and defects in the immune response to foreign antigens. Subsequently, minimal fibrosis, bile duct injury, and portal vein inflammation are caused. Therefore, to better understand the pathogenesis of PBC, it is necessary to study the key genes and pathways just mentioned.

In conclusion, the CCL5, IL7R, and TNFRSF1A genes; the immune response pathway; as well as molecular mimicry may have crucial roles in PBC. These genes and pathways may be potential new targets for treating PBC. However, further mechanistic studies of these genes and pathways are required.

Footnotes

Acknowledgment

The authors would like to acknowledge Dr. Zhenyu Yin for his critical reading and editing.

Author Disclosure Statement

The authors declare they have no competing financial interests.

References

Barrett

, Wilhite

S.E.

, Ledoux

, et al. 2013. NCBI GEO: Archive for functional genomics data sets—Update. Nucleic Acids Res. 41, D991–D995.

Bowlus

C.L.

, and Gershwin

M.E.

2014. The diagnosis of primary biliary cirrhosis. Autoimmun. Rev. 13, 441–444.

Chen

L.P.

, Zhao

, Lyu

, et al. 2016. [Environmental factors and primary biliary cirrhosis]. Zhonghua Gan Zang Bing Za Zhi, 24, 541–544.

Cordell

H.J.

, Han

, Mells

G.F.

, et al. 2015. International genome-wide meta-analysis identifies new primary biliary cirrhosis risk loci and targetable pathogenic pathways. Nat. Commun. 6, 8019.

Davis

A.P.

, Grondin

C.J.

, Lennon-Hopkins

, et al. 2015. The Comparative Toxicogenomics Database's 10th year anniversary: Update 2015. Nucleic Acids Res. 43, D914–D920.

Dong

, Li

, Tang

, et al. 2015. Multiple genetic variants associated with primary biliary cirrhosis in a Han Chinese population. Clin. Rev. Allergy Immunol. 48, 316–321.

European Association for the Study of the Liver. 2009. EASL Clinical Practice Guidelines: Management of cholestatic liver diseases. J. Hepatol. 51, 237–267.

Gene Ontology Consortium. 2006. The Gene Ontology (GO) project in 2006. Nucleic Acids Res. 34, D322–D326.

Hardie

, Green

, Jopson

, et al. 2016. Early molecular stratification of high-risk primary biliary cholangitis. EBioMedicine, 14, 65–73.

10.

Hirschfield

G.M.

, Liu

, Han

, et al. 2010. Variants at IRF5-TNPO3, 17q12-21 and MMEL1 are associated with primary biliary cirrhosis. Nat. Genet. 42, 655–657.

11.

Hirschfield

G.M.

, Liu

, Xu

, et al. 2009. Primary biliary cirrhosis associated with HLA, IL12A, and IL12RB2 variants. N. Engl. J. Med. 360, 2544–2555.

12.

Hjelmervik

T.O.

, Petersen

, Jonassen

, et al. 2005. Gene expression profiling of minor salivary glands clearly distinguishes primary Sjogren's syndrome patients from healthy control subjects. Arthritis Rheum. 52, 1534–1544.

13.

Janky

, Verfaillie

, Imrichova

, et al. 2014. iRegulon: From a gene list to a gene regulatory network using large motif and track collections. PLoS Comput. Biol. 10, e1003731.

14.

Joshita

, Umemura

, Tanaka

, et al. 2017. Genetic contribution to the pathogenesis of primary biliary cholangitis. J. Immunol. Res. 2017, 3073504.

15.

Juran

B.D.

, Hirschfield

G.M.

, Invernizzi

, et al. 2012. Immunochip analyses identify a novel risk locus for primary biliary cirrhosis at 13q14, multiple independent associations at four established risk loci and epistasis between 1p31 and 7q32 risk variants. Hum. Mol. Genet. 21, 5209–5221.

16.

Kanehisa

, and Goto

2000. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 28, 27–30.

17.

Kitamura

, Nakamoto

, Akiyama

, et al. 2002. Pathogenic roles of tumor necrosis factor receptor p55-mediated signals in dimethylnitrosamine-induced murine liver fibrosis. Lab. Invest. 82, 571–583.

18.

Kulakova

O.G.

, Bashinskaya

V.V.

, Tsareva

E.Y.

, et al. 2016. [Association analysis of cytokine receptors' genes polymorphisms with clinical features of multiple sclerosis]. Zh. Nevrol. Psikhiatr. Im. S. S. Korsakova, 116, 10–15.

19.

Lindor

K.D.

2009. A new section: Image of the month. Hepatology, 49, 344.

20.

Liston

, Enders

, and Siggs

O.M.

2008. Unravelling the association of partial T-cell immunodeficiency and immune dysregulation. Nat. Rev. Immunol. 8, 545–558.

21.

Liu

J.Z.

, Almarri

M.A.

, Gaffney

D.J.

, et al. 2012. Dense fine-mapping study identifies new susceptibility loci for primary biliary cirrhosis. Nat. Genet. 44, 1137–1141.

22.

Liu

, Invernizzi

, Lu

, et al. 2010. Genome-wide meta-analyses identify three loci associated with primary biliary cirrhosis. Nat. Genet. 42, 658–660.

23.

Luettig

, Boeker

K.H.

, Schoessler

, et al. 1998. The antinuclear autoantibodies Sp100 and gp210 persist after orthotopic liver transplantation in patients with primary biliary cirrhosis. J. Hepatol. 28, 824–828.

24.

Mells

G.F.

, Floyd

J.A.

, Morley

K.I.

, et al. 2011. Genome-wide association study identifies 12 new susceptibility loci for primary biliary cirrhosis. Nat. Genet. 43, 329–332.

25.

Miyakawa

, Tanaka

, Kikuchi

, et al. 2001. Detection of antimitochondrial autoantibodies in immunofluorescent AMA-negative patients with primary biliary cirrhosis using recombinant autoantigens. Hepatology, 34, 243–248.

26.

Nakamura

, Nishida

, Kawashima

, et al. 2012. Genome-wide association study identifies TNFSF15 and POU2AF1 as susceptibility loci for primary biliary cirrhosis in the Japanese population. Am. J. Hum. Genet. 91, 721–728.

27.

Nischalke

H.D.

, Nattermann

, Fischer

H.P.

, et al. 2004. Semiquantitative analysis of intrahepatic CC-chemokine mRNas in chronic hepatitis C. Mediators Inflamm. 13, 357–359.

28.

Selmi

, De Santis

, Cavaciocchi

, et al. 2010. Infectious agents and xenobiotics in the etiology of primary biliary cirrhosis. Dis. Markers, 29, 287–299.

29.

Shadidi

K.R.

, Aarvak

, Henriksen

J.E.

, et al. 2003. The chemokines CCL5, CCL2 and CXCL12 play significant roles in the migration of Th1 cells into rheumatoid synovial tissue. Scand. J. Immunol. 57, 192–198.

30.

Shannon

, Markiel

, Ozier

, et al. 2003. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504.

31.

Szklarczyk

, Franceschini

, Kuhn

, et al. 2011. The STRING database in 2011: Functional interaction networks of proteins, globally integrated and scored. Nucleic Acids Res. 39, D561–D568.

32.

Szklarczyk

, Franceschini

, Wyder

, et al. 2015. STRING v10: Protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 43, D447–D452.

33.

Van de Water

, Ishibashi

, Coppel

R.L.

, et al. 2001. Molecular mimicry and primary biliary cirrhosis: Premises not promises. Hepatology, 33, 771–775.

34.

Wang

J.H.

, Zhao

L.F.

, Lin

, et al. 2014. GenCLiP 2.0: A web server for functional clustering of genes and construction of molecular networks based on free terms. Bioinformatics, 30, 2534–2536.

35.

Zhang

, Ma

, Li

, et al. 2014. Gene expression profiles of peripheral blood mononuclear cells in primary biliary cirrhosis. Clin. Exp. Med. 14, 409–416.

36.

Zweers

S.J.

, de Vries

E.M.

, Lenicek

, et al. 2017. Prolonged fibroblast growth factor 19 response in patients with primary sclerosing cholangitis after an oral chenodeoxycholic acid challenge. Hepatol. Int. 11, 132–140.