This article is mentioned in:

Introduction

Biliary atresia (BA) is a rare hepatobiliary disorder in infants with poor prognosis caused by an obstruction of the bile ducts, leading to a disruption in bile flow, which can rapidly progress to biliary cirrhosis (1,2). The incidence of BA in Asia is higher than in western countries, with reports indicating a rate of 1/9,600 to 1/16,700 of live births (3). BA is reported to occur in half of all cases of neonatal cholestasis with a greater incidence in females, and associated congenital malformations occur in 3–10% of cases of BA (4). Without treatment, BA progresses into liver cirrhosis, hepatic failure and mortality in <2 years following birth (5,6). Natural mortality of BA is 100%, therefore liver transplantation is a necessary requirement for late stage BA (7); however, current medical management of patients with BA is unsatisfactory, placing an economic burden on society and family members (5,7).

The etiology of BA remains to be fully elucidated; however, previous studies have identified three factors suggested to contribute to the formation of BA. One is congenital genetic abnormalities, for example, the mutation of INVS gene, which leads to abnormal functioning of inversin, and may terminate the development of biliary ducts during the embryonic period (8). Another factor identified by various studies is viral infection (9,10), with cytomegalovirus, Epstein-Barr virus and human papilloma virus having all been isolated from the livers of infants with BA (11–13). The third factor identified was immune disorders, which appeared following the viral infection (14). Previous studies have supported the hypothesis that the viral infection triggering host immune responses may be the basis for the pathogenesis of BA (15). Leonhardt et al (16) and Schreiber and Kleinman (17) additionally agreed that BA is a virus-induced, immune-mediated disease that occurs in genetically susceptible individuals. Previous studies have demonstrated that BA was associated with T helper cell response type I-mediated portal tract inflammation and that interleukin (IL)-2, interferon γ, IL-12, tumor necrosis factor α and macrophage activation all served important roles (18,19). The current study aimed to elucidate the pathogenesis of BA using biological information analysis in relation to immunology.

Microarray analysis has been reported as an effective method to monitor alterations in gene expression and identify genes that are associated with biological processes induced by BA. Leonhardt et al (16) identified several targeted genes that were implicated in the development of BA, including CCL2, CCL5, CCR5, CXCL9, CXCL10 and IL1F5 (20). The majority of these genes were associated with the immune system or with modulating inflammation.

Previous studies have demonstrated that the expression levels of cell signaling and transcription regulatory genes, such as CIT or ActRII were upregulated during BA progression (21), and BA was observed to be a necrotizing process caused by viral cholangitis at the embryonic stage and the expression level of the type II major histocompatibility complex was increased in the livers of patients (22). However, there is little evidence indicating the relevance between human leukocyte expression and BA. Thus, the cause and pathogenetic mechanism of BA remain largely unknown.

In the present study, microarray analysis was used to monitor differentially expressed genes (DEGs) in BA samples at different time points (3, 7 and 14 days) following induction with rhesus rotavirus (RRV) compared with sham-RRV-control. Comprehensive bioinformatics analysis was used, which combined gene ontology (GO) and Kyoto Enrichment of Genes and Genome (KEGG) pathway information, in order to provide a more thorough understanding of the biological mechanisms of BA. The current study aimed to predict the hub genes that are most likely to be associated with BA, and to identify their molecular mechanisms in BA. This may aid in the development of novel therapeutic strategies for BA.

Materials and methods

Affymetrix microarray data

The microarray data were obtained from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) database from the National Center for Biotechnology Information (Bethesda, MD, USA) (23), which is the largest public database of its type. The gene expression profiles were extracted from the GEO database (ID: GSE46967) based on the GPL6246 Affymetrix Mouse Gene 1.0 ST array (Affymetrix, Inc., Santa Clara, CA, USA). The experimental mice (18 tissue samples) were injected intraperitoneally with 1.5×106 units of rhesus rotavirus (RRV) or 0.9% normal saline (NS; control) within 24 h. The microarray data were collected from the experimental group and the control group at 3, 7 and 14 days after RRV or saline injection. Each experimental condition contained three sets of samples (three treated with RRV three with normal saline).

A total of 28,816 probes were mapped with the gene names labelled in the GPL6246 annotation platform, and a log2 transformation was performed on 25,299 genes (24). The limma package (http://bioinf.wehi.edu.au/limma) in R language was used to select the DEGs of the case-group compared with the control-group (25). The P-values were adjusted by the Benjamin and Hochberg method, and the FDR-value based on the multtest package of 0.05 was used as the cut-off criterion (26).

Hiarchical clustering of DEGs

The expression values of DEGs in each group were selected according to the probe information from the downloaded file. The R language pheatmap package (http://cran.r-project.org/web/packages/pheatmap/index.html) was used for hierarchical clustering analysis (27) of DEGs based on Euclidean distance (28), and the results were presented in heat maps. Genes with similar expression levels could be collected together by hierarchical clustering for a convenient search. Samples clustered based on the gene expression value can aid in determining whether the screening DEGs have sample specificity or not, which means the selected sample types (case or control) can be identified based on the differential gene expression information.

Common DEGs of three time points

The expression levels of selected genes from samples at 3, 7 and 14 days were compared, in order to identify the common DEGs (29).

GO analysis and pathway enrichment analysis

GO analysis has become a commonly used approach for functional studies of large-scale genomic or transcription data (30). The KEGG pathway database contains information of how molecules or genes are networked, which is complementary to the majority of the existing molecular biology databases containing the information of individual genes (31). The Database for Annotation, Visualization and Integrated Discovery (DAVID) bioinformatics resources consist of an integrated biological knowledge base and analytical tools aimed at systematically extracting biological meaning from large gene or protein lists (32). DAVID was used to analyze the enriched functions and pathways of the common DEGs at the three time points with P<0.05.

In addition, the expression values of DEGs that participated in the most significant GO functions and KEGG pathways were compared.

Results

Screening and hierarchical clustering analysis

In order to investigate DEGs in mice with BA at 3, 7 and 14 days following viral induction, the publicly available microarray dataset GSE46967 from GEO was used. A total of 306 DEGs were selected in the 3 day samples, 721 DEGs in the 7 day samples and 370 DEGs in the 14 day samples (Fig. 1). From the heat maps (Fig. 2), the results of hierarchical clustering analysis revealed that the case- and control-groups can be separated by the selected DEGs in each group, which means the expression patterns of DEGs screened were significant, and they can be used to distinguish case and control samples. A total of 74 common DEGs were identified between the samples from the different time points (Fig. 3A), and the expression values of all of these common DEGs were upregulated, as compared with in the control samples(Fig. 3B).

Figure 1 Cassette figure of expression values of DEGs. Case-group is marked by a red horizontal line in each sample.

Figure 2 Hierarchical clustering heat maps of the DEGs. Heat maps of DEGs at (A) 3, (B) 7 and (C) 14 days. Color gradient from blue-orange represents expression level (case-group and control-group) from up- to downregulation. Case-group is marked by a red horizontal line in each sample. DEGs, differentially expressed genes.

Figure 3 Expression values of the selected 74 common DEGs at three time points. (A) Venn diagram of DEGs from the three time points. (B) Hierarchical clustering heat maps of DEGs. Orange represents upregulation, and blue represents downregulation. DEGs, differentially expressed genes.

Functional annotation analysis of common DEGs

To investigate the alterations in the function of selected common DEGs, DAVID was used to identify the significant GO categories in the biological process. A total of 12 GO terms were observed to be enriched among the 74 DEGs (Table I), of which several significantly enriched functions were associated with immune function, including the defense response, leukocyte migration, cell chemotaxis and leukocyte chemotaxis. C3AR1, CCL3, C5AR1, CXCL5, FPR1, CXCR2, CCL5, ITGAM, CCL6, CXCL13, CSF3R, IL1B and FCER1G were the genes enriched in the significant functions.

Table I Function enrichment analysis of the common DEGs.

Table I

Function enrichment analysis of the common DEGs.

Term	Count	FDR
GO:0042330~taxis	13	1.67×10−13
GO:0006935~chemotaxis	14	1.67×10−13
GO:0006955~immune response	21	1.67×10−13
GO:0007626~locomotory behavior	14	4.38×10−09
GO:0006952~defense response	17	5.54×10−09
GO:0007610~behavior	14	3.07×10−06
GO:0030595~leukocyte chemotaxis	6	6.90×10−05
GO:0060326~cell chemotaxis	6	6.90×10−05
GO:0006954~inflammatory response	10	1.98×10−04
GO:0030593~neutrophil chemotaxis	5	6.03×10−04
GO:0009611~response to wounding	11	7.79×10−04
GO:0050900~leukocyte migration	6	7.86×10−04

[i] DEGs, differentially expressed genes; GO, gene ontology; FDR, false discovery rate.

Pathway analysis of common DEGs

To further understand the pathways of the DEGs screened in the current study, DAVID was used to identify the significant pathways. Two pathways (Table II), were identified to be enriched, the 'cytokine-cytokine receptor interaction' and 'chemokine signaling pathway'. The 'cytokine-cytokine receptor interaction' pathway was the most significant pathway with ten genes involved, as follows: IL1R2, CCL3, CXCL5, CXCL13, CCR1, IL1B, CSF3R, CXCR2, CCL5 and CCL6.

Table II Pathways involved with common DEGs.

Table II

Pathways involved with common DEGs.

Kyoto enrichment of genes and genomes	Gene number	FDR
mmu04060: Cytokine-cytokine receptor interaction	10	0.003256
mmu04062: Chemokine signaling pathway	9	0.003469

[i] DEGs, differentially expressed genes; FDR, false discovery rate.

Expression values of DEGs participating in significant functions and pathways

A total of 13 DEGs were identified to be associated with the significant functions, and 10 DEGs with the significant pathways. Following the comparison of the two types of DEGs, the following 6 common DEGs were identified: CCL3, CXCL5, CXCL13, CXCR2, CCL5 and CCL6. The expression profile of these 6 DEGs in each of the three samples were analyzed (Fig. 4A). It is observed in the figure that all 6 of the DEGs in the three sample types were significantly overexpressed, and CXCL13 was identified to have a consistent overexpression from the 3rd day to the 14th day. CCL3, CCL5 and CCL6; however, demonstrated the highest overexpression on the 7th day, but were decreased on the 14th day. CXCL5 and CXCR2 presented the highest overexpression on the 3rd day (Fig. 4B).

Figure 4 Expression values of the selected 6 DEGs that participate in the significant functions and pathways. (A) Hierarchical clustering heat maps of the 6 DEGs. (B) Expression values of the 6 DEGs at different time points. DEGs, differentially expressed genes.

Discussion

BA is an infantile disorder in which there is obliteration or discontinuity of the extrahepatic biliary tree. If left untreated, it progresses to cirrhosis, liver failure and eventually to mortality, and the current treatment strategies are unsatisfactory (2,5). Previous studies using histological sections of extrahepatic biliary tract tissue from patients, in addition to experimental studies using the BA mice model, have demonstrated that an inflammatory process appears to contribute to the pathogenesis of BA (18,33). Previous studies that further identify the relevant inflammatory cascade revealed an upregulation in pro-inflammatory cytokines (20). A total of 306 DEGs were identified in the 3 day samples, 721 in the 7 day samples and 370 in the 14 day samples in mice with BA following injection with rhesus rotavirus. In total, 74 common DEGs were selected from the three sample types at different time points, which function predominantly in the immune biological process, including the defense response, leukocyte migration, cell chemotaxis and leukocyte chemotaxis. The 'cytokine-cytokine receptor interaction' and 'chemokine signaling pathway' were identified as the most significantly enriched pathways. A total of six significantly overexpressed genes (CCL3, CXCL5, CXCL13, CXCR2, CCL5 and CCL6) associated with the immune response and cytokine pathways at the three time points were screened, which were investigated in the significant functions and pathways.

Chemokines are small proteins that are vital in immune and inflammatory reactions. Induction of leukocyte migration is the predominant function of the chemokines, which additionally affect angiogenesis, collagen production and the proliferation of hematopoietic precursors (34). CCL3 (MIP-1α) is a chemokine produced by activated T-cells, basophilic granulocytes and type 2 T helper cells (35). The chemokine CCL3 has several functions, including serving a role in cell lysis mediated by natural killer cells, participating in the inflammatory response associated with immunological processes, and mediating the release of the other cytokines (36). Seki et al (37) reported that MIP-1α had the ability to gather the hepatic stellate and inflammatory cells migrating to the site of injury, so as to function on the liver fibration. Zerfaoui et al (38) demonstrated that CCL3 was able to activate T cells thus producing CCL5 (RANTES). CCL3 and CCL5 have the ability to assemble and to activate eosinophilic granulocytes to release histamine and leukotrienes so as to induce inflammatory cell adherence to the vessel wall (38). Increased levels of CCL5 have been observed in patients with primary biliary cirrhosis, and the therapeutic benefits of fibrates in this biliary disease may be due to the inhibition of this CCL5, with its resultant effects on inflammatory cell migration (39). From the results of the present study, it is suggested that the expression of CCL3 and CCL5 are elevated in a time-dependent manner until 7 days, indicating that these two molecules are critical for the pathogenesis of BA in mice. Furthermore, previous studies have demonstrated that the two cytokines are associated with the migration and infiltration of monocytes, memory T lymphocytes and natural killer cells (40–42), so it was speculated that they may be involved in the early recruitment of different subsets of leukocytes in the formation of BA.

The CXCR2 gene was first cloned by Murphy and Tiffany (43) and the chemokine CXCR2, which is produced on the surface of neutrophils and natural killer cells, is also a member of the G-protein superfamily. A previous study demonstrated that CXCR2 may combine with the chemokine CXCL5 in order to activate the downstream effector molecules, and may subsequently mediate signal transduction and cell migration via the allosteric G-protein (44). CXCL5 is a small-molecule type secretory protein, which has a direct chemotaxis for neutrophils and is able to mediate the infiltration of neutrophils, thus promoting the metastasis of liver cancer cells (45). Li et al (46) reported that CXCL5 contributed to the formation of tumor vessels by activating certain signaling pathways, such as the protein kinase B, signal transduction and activating transcription factor pathways (47). In the present study, the variation tendency of the CXCL5 expression level was similar with that of CXCR2 and their expression peaked at 3 days. This phenomenon may indicate that the chemokine receptor CXCR2 and CXCL5 may combine, in order to affect inflammation in mice with BA by collecting leukocytes.

CXCL13 is a chemokine produced by B-cells and is critical for targets seeking B-cell and B-cell zone formation in lymphatic tissue (48). A previous study indicated that modulation of the inflammatory signaling cascade, particularly cytokine secretion, may attenuate the severity of BA in the murine model. Shivakumar et al demonstrated that RRV-infected, interferon γ-deficient Bulb/c mice had a lower mortality rate than wild type animals (48). Anesl et al (49) demonstrated that the natural antibody production in mice deficient in CXCL13 was markedly reduced, and no response was triggered to a bacterial antigen by intraperitoneal injection in this case. In addition, the roles of CXCL13 have been observed in other immune system-associated diseases. Hjelmstrom et al (50) observed that CXCL13 was overexpressed in transgenic mice suffering from chronic inflammation induced by lymphotoxin-α (50). Based on the data from the current study, the expression levels of CXCL13 were increased until 14 days in mice with BA. Elevated levels of CXCL13 is not unique to active BA, but it may serve as a marker of disease activity or response to treatment in combination with other diagnostic markers.

Notably, low expression levels of CCL6 were observed in BA-positive mice. It was hypothesized that the chemokine CCL6 may serve an important function in the progression of BA, however further study is required in order to further elucidate the mechanisms.

In conclusion, CCL3, CCL5, CXCL5, CXCR2, CXCL13 and CCL6, are suggested to be involved in the pathogenesis of BA, but serve different roles. The current study provided a basis for future experiments with the potential use of RRV-injected mice, which may enhance understanding of the pathogenesis of BA and lead to the development of novel therapeutic strategies.

References