Evolving Up‐regulation of Biliary Fibrosis–Related Extracellular Matrix Molecules After Successful Portoenterostomy

Abstract

Successful portoenterostomy (SPE) improves the short‐term outcome of patients with biliary atresia (BA) by relieving cholestasis and extending survival with native liver. Despite SPE, hepatic fibrosis progresses in most patients, leading to cirrhosis and a deterioration of liver function. The goal of this study was to characterize the effects of SPE on the BA liver transcriptome. We used messenger RNA sequencing to analyze global gene‐expression patterns in liver biopsies obtained at the time of portoenterostomy (n = 13) and 1 year after SPE (n = 8). Biopsies from pediatric (n = 2) and adult (n = 2) organ donors and other neonatal cholestatic conditions (n = 5) served as controls. SPE was accompanied by attenuation of inflammation and concomitant up‐regulation of key extracellular matrix (ECM) genes. Highly overexpressed genes promoting biliary fibrosis and bile duct integrity, such as integrin subunit beta 6 and previously unreported laminin subunit alpha 3, emerged as candidates to control liver fibrosis after SPE. At a cellular level, the relative abundance of activated hepatic stellate cells and liver macrophages decreased following SPE, whereas portal fibroblasts (PFs) and cholangiocytes persisted. Conclusion: The attenuation of inflammation following SPE coincides with emergence of an ECM molecular fingerprint, a set of profibrotic molecules mechanistically connected to biliary fibrosis. The persistence of activated PFs and cholangiocytes after SPE suggests a central role for these cell types in the progression of biliary fibrosis.

Abbreviations
• aHSC
• activated hepatic stellate cell
• aPF
• activated portal fibroblast
• BA
• biliary atresia
• CKAP2L
• cytoskeleton associated protein 2 like
• CLDN19
• claudin 19
• CTHRC1
• collagen triple helix repeat containing 1
• DC
• disease control
• DEG
• differentially expressed gene
• DR
• ductular reaction
• ECM
• extracellular matrix
• ENTPD2
• ectonucleoside triphosphate diphosphohydrolase 2
• FAP
• fibroblast activation protein alpha
• GO
• Gene Ontology
• GSVA
• Gene Set Variation Analysis
• HMMR
• hyaluronan mediated motility receptor
• IHC
• immunohistochemistry
• IQR
• interquartile range
• ISH
• in situ hybridization
• ITGB6
• integrin subunit beta 6
• KC
• Kupffer cell
• KRT22
• keratin 22
• LAMA3
• laminin subunit alpha 3
• LAMC2
• laminin subunit gamma 2
• LT
• liver transplantation
• MMP7
• matrix metallopeptidase 7
• mRNA
• messenger RNA
• MROH7
• maestro heat like repeat family member 7
• NC
• normal control
• NCAPG
• non‐SMC condensin I complex subunit G
• NK
• natural killer
• PE
• portoenterostomy
• RNA‐Seq
• RNA sequencing
• RPKM
• reads per kilobase million
• SAMacs
• scar‐associated macrophages
• SCTR
• secretin receptor
• SPE
• successful portoenterostomy
• ssGSEA
• single‐sample gene‐set enrichment analysis
• STMN2
• stathmin 2
• TGF‐β
• transforming growth factor‐β
• TMEM178B
• transmembrane protein 178B
• VWA2
• von Willebrand factor A domain containing 2
• WNT10A
• Wnt family member 10A

Biliary atresia (BA) is a rare neonatal liver disease characterized by progressive hepatic fibroinflammatory injury leading to cholestasis and fibrotic obstruction of extrahepatic bile ducts. The natural history of BA is grim, as accumulating hepatic damage leads to end‐stage liver disease and death during the first few years of life.⁽¹⁾

The initial surgical treatment of BA is portoenterostomy (PE), in which the obstructed extrahepatic biliary system is removed and replaced with a loop of intestine. Successful PE (SPE), achieved in most patients with BA, normalizes serum bilirubin levels and improves the short‐term outcome by relieving cholestasis and extending survival with native liver.⁽¹⁾ Despite SPE, the hepatic fibro‐obliterative process continues and extends to intrahepatic bile ducts, eventually leading to cirrhosis and deterioration of liver function.⁽²⁾ There are no effective treatments to counteract progressive liver fibrosis after SPE. Consequently, BA remains the most frequent indication for liver transplantation in children worldwide, with 75% of patients requiring this procedure by the age of 20 years.⁽³⁾

Recent studies have profiled histopathological and clinical aspects of post‐SPE fibrosis.⁽⁴⁾ The early neutrophil‐dominated acute inflammation resolves within the first year, but a smaller lymphocyte infiltrate persists, and fibrosis progresses in most patients.⁽⁵⁾ Increased numbers of mesenchymal extracellular matrix (ECM)–producing myofibroblasts and bile duct–forming cholangiocytes are evident in both PE and SPE samples.⁽⁵⁾ A candidate gene analysis revealed that classic molecular markers of active fibrogenesis prevail in BA livers after SPE, despite relief of cholestasis.⁽⁶⁾ However, the effect of SPE on global gene expression in BA liver has not been assessed. Similarly, the molecular factors and cell types driving post‐SPE liver fibrosis remain unknown.

To characterize global changes in liver gene expression and to identify candidate molecules and cell types responsible for the persistent liver injury after SPE, we performed RNA sequencing (RNA‐Seq) analysis on liver samples obtained both at the time of PE and at 1‐year follow‐up in patients who underwent SPE (i.e., normalized serum bilirubin as a consequence of PE). Liver biopsies from patients with other neonatal cholestatic conditions as well as pediatric and adult liver donors served as controls. We hypothesized that molecules highly overexpressed both in the samples obtained at the time of PE and at follow‐up are critical for the progressive liver injury after SPE.

Materials and Methods

Patients and Controls

Human tissue biopsy and blood samples were collected from 13 patients with BA at PE and 8 patients at 1 year following SPE at the Helsinki University Central Hospital. We included liver biopsies from 5 patients with other neonatal cholestatic disorders as disease controls (DCs) and liver biopsies from 2 adult and 2 children organ donors as normal controls (NCs). Patient and control characteristics are provided in Table 1. BA diagnosis was confirmed by histopathological assessment of the portal plate.

1 TABLEBaseline Characteristics of Patients in PE, SPE, DC, and NC Groups

	PE	SPE	DC	NC (Children)	NC (Adults)
Patients, n	13	8	5	2	2
Female, n (%)	6 (46%)	4 (50%)	1 (20%)	0 (0%)	1 (50%)
Type of BA, n (%)
1	1 (8%)	1 (13%)	NA	NA	NA
3	12 (92%)	7 (87%)	NA	NA	NA
Age at PE, days	40 (24‐66)	35 (24‐60)	NA	NA	NA
COJ,* n (%)	10 (77%)	8 (100%)	NA	NA	NA
Time from PE to COJ, days (range)	38 (22‐108)	56 (26‐106)	NA	NA	NA
LT in 5 years, n (%)	4 (31%)	2 (25%)	NA	NA	NA
Age at LT, days (range)	510 (239‐898)	910 (898‐921)	NA	NA	NA
Age at liver biopsy, days (range)	40 (24‐66)	470 (360‐500)†	40 (37‐59)	543 (356‐730)	37 (32‐42), years
Liver biochemistry
Bilirubin, µmol/L	155 (112‐175)	6 (5‐14)†	94 (92‐127)	NA	NA
Conjugated bilirubin, µmol/L	107 (92‐143)	3 (2‐10)†	86 (66‐98)	NA	NA
ALT, U/L	136 (65‐161)	88 (38‐142)	55 (50‐119)	NA	NA
Bile acids, µmol/L	201 (143‐246)	63(30‐135)†	190 (97‐207)	NA	NA

Note::

Variables are expressed as medians and IQRs. Mann‐Whitney U test or Fisher’s exact test was used for comparisons.

^*Bilirubin < 20 µmol/L.

^†P < 0.05 PE versus SPE.

Abbreviations: ALT, alanine aminotransferase; COJ, clearance of jaundice; LT, liver transplantation; NA, not applicable.

Clinical Data and Liver Biopsies

Routine blood samples were used to measure plasma levels of bilirubin, alanine aminotransferase, and total bile acids. Wedge liver biopsies were taken during PE, and ultrasound‐guided core‐needle biopsies were obtained under general anesthesia during endoscopic variceal surveillance at follow‐up. Part of each biopsy was fixed in formalin and embedded in paraffin to be used for immunohistochemistry (IHC) staining, and the remainder was placed in RNAlater solution (Ambion; Thermo Fisher Scientific, Waltham, MA) and stored at −80°C.

IHC Staining

The formalin‐fixed, paraffin‐embedded (FFPE) sections were deparaffinized, hydrated, and treated with Target retrieval solution pH 9 (Dako, Glostrup, Denmark), as reported.⁽⁷⁾ Commercially available antibody for laminin subunit alpha 3 (LAMA3) (rabbit polyclonal, PA5‐52665; Thermo Fisher Scientific) was used at a dilution of 1:2,000 along with the Novolink Polymer Detection system (Leica Biosystems, Newcastle, United Kingdom). Primary antibody was incubated overnight at 4°C. Images were generated using a 3DHISTECH Panoramic 250 FLASH II digital slide scanner (3DHISTETCH Kft., Budapest, Hungary) at the Genome Biology Unit, supported by HiLIFE and the Faculty of Medicine, University of Helsinki, and Biocenter Finland.

RNA In Situ Hybridization

RNA in situ hybridization (ISH) was performed as described.⁽⁷⁾ In short, FFPE tissue sections were processed for target detection using RNAscope Multiplex Fluorescent Reagent Kit version 2 (#323100; Advanced Cell Diagnostics, Newark, CA) according to the manufacturer’s instructions. Hs‐LAMA3 (530681), Hs‐keratin 7–C3 (550151‐C3), negative control probe diaminopimelate biosynthetic (#320871), and positive control probe (#320861) were used. TSA Plus Cyanine 3 at a dilution of 1:1,500 (NEL744001KT; PerkinElmer, Waltham, MA) was used for signal detection. Images were taken using a 3DHISTECH Panoramic 250 FLASH II digital slide scanner at the Genome Biology Unit, University of Helsinki.

Messenger RNA Extraction

Fresh frozen tissue samples immersed in RNAlater solution were homogenized using the Precellys lysing kit and tissue homogenizer (Bertin Technologies, France). Total RNA was extracted using the RNeasy Mini Kit (Qiagen, Hilden, Germany) according to the supplier’s instructions. RNA integrity was measured using Agilent 2100 Bioanalyzer (or TapeStation) (Agilent Technologies, Santa Clara, CA). Messenger RNA (mRNA) was purified using magnetic beads, and next‐generation sequencing was performed using the Illumina NextSeq Sequencing System and the NEBNext Ultra Directional RNA Library Prep Kit for Illumina for complementary DNA (cDNA) library preparation.

cDNA Sequencing and Differential Expression Analysis

Sequence read data that passed quality control checks were mapped to the human transcriptome using HISAT2 (hierarchical indexing for spliced alignment of transcripts) spliced aligner,⁽⁸⁾ and gene expression level was quantified using StringTie.⁽⁹⁾ Differential expression analysis was performed using DESeq2 Bioconductor package.⁽¹⁰⁾ To gain further insights into the biology and functions of the differentially expressed genes (DEGs), the Gene Ontology (GO) pathway analysis tool (GOstats) was used to identify biological pathways that are significantly enriched.⁽¹¹⁾ The clusterProfiler package was used to visualize the GO analysis results.⁽¹²⁾ Hypergeometric tests with the Benjamini‐Hochberg false discovery rate were used to adjust the P value.

Single‐Sample Gene‐Set Enrichment Analysis

Single‐sample gene‐set enrichment analysis (ssGSEA) implemented in the Gene Set Variation Analysis (GSVA) R Bioconductor package was used to estimate the relative abundance of individual cell types from the RNA‐Seq data.⁽¹²⁾ Published gene‐signature markers for activated hepatic stellate cells (aHSCs),⁽¹³⁾ activated portal fibroblasts (aPFs),⁽¹⁴⁾ cholangiocytes,⁽¹⁵⁾ natural killer (NK) cells and CD8 + T cells,^(16,17) and Kupffer cells (KCs) and scar‐associated macrophages (SAMacs)⁽¹⁸⁾ were used for ssGSEA.

RNA‐Seq data were filtered to exclude genes in which fewer than 10 samples had reads per kilobase million (RPKM) values > 0.5. Log‐transformed RPKM value was used for ssGSEA. The ssGSEA score for each cell type signature was obtained through the GSVA R Bioconductor package.⁽¹²⁾

Statistical Analysis

Mann‐Whitney U test or Fisher’s exact test was used for comparisons. Comparison of changes among multiple groups was performed using either a Kruskal‐Wallis test, Dunn test, Tukey test, or repeated measures one‐way analysis of variance. P value < 0.05 was considered statistically significant. All IHC and ISH stainings were repeated in a minimum of 3 patients, and representative images are displayed in the figures.

Ethical Permission

The study protocol was approved by the hospital ethics committee (#345/03/1372008) and complies with the ethical guidelines of the 1975 Declaration of Helsinki. Institutional approval was renewed on July 21, 2017 (§68 HUS/149/2017). Informed consent was obtained from all participants or their legal guardians. No donor organs were obtained from executed prisoners or institutionalized persons.

Results

Patient Characteristics

A total of 13 newborns with BA underwent PE at a median age of 40 days (Table 1); 8 of these patients who cleared their jaundice and reached normalized serum bilirubin (<20 µmol/L) after SPE underwent a follow‐up native liver biopsy at the median age of 470 days. The PE and SPE groups had similar baseline characteristics (Table 1). The PE and DC groups underwent liver biopsy at the same age and had similar increases in serum bilirubin and bile acid levels.

ECM Genes Are Activated in BA at the Time of PE

To survey global gene expression in BA liver at the time of PE, we performed RNA‐Seq on a total of 22 liver samples (PE [n = 13], DC [n = 5], and NC [n = 4]). We identified 6,189 statistical DEGs in BA liver biopsies obtained at PE as compared with NC. To exclude the nonspecific effects of cholestasis, we omitted 1,384 DEGs with overlapping differential expression between DC and NC samples, resulting in 4,805 DEGs with BA‐specific liver expression. To further constrain the number of DEGs and identify changes that were more likely to be clinically relevant, we focused on DEGs with more than a five‐fold change in all of our subsequent analyses (Fig. 1A).

1 FIG.. Activation of ECM genes at PE. (A) Flow chart of DEG analysis. DEGs with differential expression between DC and NC samples as well as DEGs with < five‐fold difference between PE and NC were omitted. (B) Volcano plot of 4,805 BA‐specific DEGs. DEGs with adjusted P value < 0.05 and >five‐fold up‐regulation (log2 > 2.32) were selected for subsequent GO term analysis. (C) GO term analysis of the 273 BA‐specific up‐regulated DEGs. The 20 most significantly affected GO pathways are shown, arranged by descending GeneRatio (number of affected genes divided by the number of genes in the pathway) on the x‐axis. Dot size signifies the number of affected genes in the GO pathway and color the adjusted P value of GO pathway enrichment.

Volcano plot analysis of BA‐specific genes identified 321 DEGs with more than a five‐fold change as compared with NC samples (Fig. 1B). These 321 DEGs were subjected to GO enrichment analysis, assigning genes to predefined pathways depending on their functional and spatial characteristics. The up‐regulated DEGs were central components of ECM production and structure, confirming that ECM genes are abnormally activated in BA at the time of PE (Fig. 1C), a finding that agrees with a recent report by Luo et al.⁽¹⁹⁾

Molecular Analysis of SPE Samples Demonstrated Changes in Liver Metabolism and Attenuation of Inflammation

Next, we leveraged our unique follow‐up biopsies to assess changes in gene expression elicited by SPE by comparing RNA expression levels of BA‐specific genes in biopsies obtained at PE (n = 13) and after SPE (n = 8). Initial principal component analysis suggested a differential gene‐expression profile between SPE and PE samples (Fig. 2A). Further analysis identified 276 BA‐specific genes exhibiting more than a five‐fold difference in expression between PE and SPE samples (Fig. 2B).

2 FIG.. Effect of SPE on global liver gene expression. (A) Principal component analysis showing the global gene expression of SPE follow‐up (blue dots) and at‐PE samples (red dots). (B) Volcano plot showing DEGs between SPE and PE livers. DEGs with adjusted P value < 0.05 and >five‐fold up‐regulation (log2 > 2.32) were selected for subsequent GO term analysis. (C) GO term analysis of up‐regulated DEGs between SPE and PE livers. (D) GO term analysis of down‐regulated DEGs between SPE and PE livers. In each dot plot, the 20 most significantly affected GO pathways are shown, arranged by descending GeneRatio. Dot size signifies the number of affected genes in the GO pathway and color the adjusted P value of GO pathway enrichment. Abbreviations: PC, principal component; OHC, organic hydroxy compound.

GO term analysis was performed separately on DEGs, depending on whether their expression was up‐regulated (n = 88) or down‐regulated (n = 188) following SPE. Analysis of the up‐regulated genes showed activation of pathways involved with key liver metabolic functions, such as hydroxylation (organic hydroxy compound metabolism—a pivotal initial reaction in a multitude of metabolic processes, including bile acid metabolism) and energy metabolism (Fig. 2C). This may reflect the fact that follow‐up samples were obtained from older individuals than those at PE. On the other hand, it is tempting to speculate that BA livers possess inherent abnormalities in metabolic process and energy metabolism, as suggested by studies showing increased energy expenditure after PE.⁽²⁰⁾

The DEGs with down‐regulated expression after SPE primarily included pathways involved in immune function and defense response, consistent with attenuation of acute inflammatory response (Fig. 2D). These findings recapitulate the outcomes obtained from the histological and candidate gene analysis of follow‐up samples after SPE.^(5,6)

Molecular Fingerprint of ECM Genes Highly Expressed in SPE Livers

To further characterize the effects of SPE on hepatic gene expression, we focused on BA‐specific genes with aberrant expression after SPE. We identified a total of 182 DEGs with more than a five‐fold difference as compared with NC (Fig. 3A). As with the PE samples, GO term analysis of up‐regulated DEGs showed pronounced activation of ECM genes in the SPE samples (Fig. 3B).

3 FIG.. Genes with aberrant expression after SPE. (A) Volcano plot showing DEGs between SPE and NC livers. DEGs with adjusted P value < 0.05 and >five‐fold up‐regulation (log2 > 2.32) were selected for subsequent GO term analysis. (B) GO term analysis of down‐regulated DEGs between SPE and NC livers. The 20 most significantly affected GO pathways are shown, arranged by descending GeneRatio. Dot size signifies the number of affected genes in the GO pathway and color the adjusted P value of GO pathway enrichment.

Notably, expression of ECM‐related genes evolved as the result of SPE. As compared with the 26 ECM‐related DEGs at PE (Supporting Fig. S1), 12 genes (collagen type X alpha 1 chain [COL10A1], collagen type XV alpha 1 chain [COL15A1], collagen type XXII alpha 1 chain [COL22A1], collagen type IX alpha 2 chain [COL9A2], matrix metallopeptidase 7 [MMP7], integrin subunit beta 6 [ITGB6], collagen triple helix repeat containing 1 [CTHRC1], ectonucleoside triphosphate diphosphohydrolase 2 [ENTPD2], Wnt family member 10A [WNT10A], von Willebrand factor A domain containing 2 [VWA2], LAMA3, and laminin subunit gamma 2 [LAMC2]) maintained their high expression after SPE (Fig. 4). In addition to collagens (Fig. 4A), established markers, and constituents of fibrosis, this list included genes previously associated with biliary fibrosis and BA (Fig. 4B), such as MMP7^(21,22) and ITGB6, a potential promoter of portal fibrosis in BA by cholangiocyte transforming growth factor‐β (TGF‐β) signaling.⁽²³⁾ Other genes previously connected to liver fibrosis (Fig. 4C) include CTHRC1, a putative activator of HSCs and promotor of liver fibrosis by activation of TGF‐β signaling⁽²⁴⁾; ENTPD2, expressed in periductal myofibroblasts and HSCs in models of liver fibrosis^(25,26); and WNT10A, a gene regulating liver fibrosis by activation of HSCs.⁽²⁷⁾ The potential role of VWA2 in liver fibrosis remains unknown (Fig. 4D). We hypothesize that this “ECM molecular fingerprint” associated with TGF‐β signaling and biliary fibrosis plays a central role in the persistent and progressive fibrosis after SPE.

4 FIG.. Box plots showing mRNA expression of ECM genes highly expressed in SPE, PE, and NC groups: Collagens (A), MMP7 and ITGB6 (B), CTHRC1, ENTPD2, and WNT10A (C), VWA2 (D), and LAMA3 and LAMC2 (E). Expression levels are depicted as RPKM. P values were calculated by Wilcoxon rank‐sum test. Black dots represent expression of independent samples; boxes represent the interquartile range (IQR); and whiskers represent the first and fourth quartile. Solid lines show the mean expression. *P < 0.05.

Persistently Highly Overexpressed Genes Are Candidate Drivers of Liver Injury After SPE

Next, we focused on the genes with highest expression at follow‐up. When the fold change limit was set to 10 between SPE and NC samples, the number of DEGs decreased to 13 (LAMA3, ITGB6, WNT10A, fibroblast activation protein alpha [FAP], secretin receptor [SCTR], stathmin 2 [STMN2], hyaluronan mediated motility receptor [HMMR], cytoskeleton associated protein 2 like [CKAP2L], claudin 19 [CLDN19], keratin 22 [KRT222], maestro heat like repeat family member 7 [MROH7], non‐SMC condensin I complex subunit G [NCAPG], and transmembrane protein 178B [TMEM178B]). Interestingly, all of these genes were already more than 10‐fold overexpressed at PE.

Three of these genes (LAMA3, ITGB6, and WNT10A) were part of the ECM molecular fingerprint described previously. Of the remaining 10 DEGs, at least four were previously connected to liver fibrosis (Fig. 5A): (1) FAP, a fibroblast activator physically interacting with ITGB6⁽²⁸⁾; (2) SCTR, which promotes liver fibrosis through changes in TGF‐β1 by stimulating proliferating cholangiocytes⁽²⁹⁾; (3) STMN2, a gene playing a role in the activation of HSCs⁽³⁰⁾; and (4) HMMR, a gene overexpressed in portal myofibroblasts of cirrhotic liver.⁽³¹⁾ The possible contribution of the remaining six genes (CKAP2L, CLDN19, KRT222, MROH7, NCAPG, and TMEM178B) to liver injury has not been assessed (Fig. 5B).

5 FIG.. Box plots of mRNA expression of DEGs highly overexpressed after SPE. (A) DEGs previously connected to liver fibrosis. (B) Other DEGs. Expression levels are depicted as RPKM. P values were calculated by Wilcoxon rank‐sum test. Black dots represent expression of independent samples; boxes represent the IQR; and whiskers represent the first and fourth quartile. Solid lines show the mean expression. *P < 0.05.

Expression and Localization of LAMA3 Suggests Its Significance for Post‐SPE Fibrosis

One of the more intriguing candidates was LAMA3, a gene coding the α subunit of laminin‐322, an ECM component critical for maintaining bile duct structure under conditions of cholangiocyte stress.⁽³²⁾ In our analysis, LAMA3 was more than 10‐fold up‐regulated both at PE and after SPE (Fig. 4E). In addition to LAMA3, another component of laminin‐322, LAMC2, was more than five‐fold up‐regulated in both PE and SPE samples (Fig. 4E).

To delineate the expression pattern of LAMA3, we performed ISH and IHC analysis. LAMA3 was specifically expressed in ductular cholangiocytes and in areas of ductular reaction (DR) in both PE and SPE samples as compared with normal liver controls (Fig. 6A‐C). The specificity of LAMA3 expression in ductular cholangiocytes was verified with LAMA3‐CK7 double‐ISH showing a consistent co‐expression of LAMA3 with the established cholangiocyte marker CK7 (Fig. 6A‐C).

6 FIG.. Abnormal expression of LAMA3 in neo‐ductular areas of BA livers. Liver samples from NC (A), PE (B), and SPE (C) were subjected to ISH or IHC. Expression of LAMA3 protein and LAMA3 as well as CK7 mRNA are shown in the left and right panels, respectively. Increased LAMA3 expression with a consistent co‐expression of CK7 in ductular cholangiocytes can be seen within areas of DR in the PE and SPE samples. Brown indicates the positive staining in IHC panels. The black arrowhead indicates increased LAMA3 protein expression, and the white arrowhead indicates the LAMA3 mRNA expression in the ductular cholangiocytes. Bars: IHC, 100 μm; ISH, 10 μm.

Relative Abundance of aHSCs and Liver Macrophages Decreases Following SPE

Next, we used ssGSEA to estimate the changes in the relative abundance of cell types associated with liver fibrosis among NC, PE, and SPE samples. PE samples exhibited increased relative abundance of effector cell types in liver fibrosis, such as aPFs and aHSCs (Fig. 7A), with a concomitant decrease in the number of CD8 + T cells (Fig. 7B). These findings are in keeping with recent results by Luo et al.⁽¹⁹⁾ However, we did not observe significant changes in the relative abundance of NK cells (Fig. 7C). In addition, we extended our analysis to cells of the mononuclear phagocyte system, known to play a role in liver fibrosis.⁽¹⁸⁾ These analyses revealed an increase in the relative abundance of SAMacs at PE (Fig. 7D).

7 FIG.. SPE induces changes in relative abundance of individual cell types. ssGSEA was used to estimate the relative abundance of individual cell types from the RNA‐Seq data. Box plots show changes in relative abundance of aHSCs and aPFs (A), CD8+ T cells (B), NK cells (C), KCs, SAMacs 1, and SAMacs 2 (D), and cholangiocytes (E) in NC, PE, and SPE groups. Previously published gene‐expression profiles for aHSCs,(13) aPFs,(14) cholangiocytes,(15) NK cells and CD8+ T cells,(16,17) and KCs and two subpopulations of SAMacs(18) were used for ssGSEA using the GSVA R Bioconductor package.(12) Enrichments are represented as ssGSEA scores. Black dots represent expression of independent samples; boxes represent the IQR; and whiskers represent the first and fourth quartile. Solid line shows the mean expression. P values were calculated by Wilcoxon rank‐sum test. *P < 0.05; **P < 0.01.

Finally, we assessed the effects of SPE on the relative abundance of cell types. Interestingly, we observed a statistically significant reduction in relative abundancies of aHSCs, SAMacs, and KCs as compared with PE samples (Fig. 7A,D), whereas the relative abundance of aPFs (Fig. 7A) did not significantly decrease after SPE. The relative abundance of cholangiocytes after SPE trended progressively upward as compared with PE and NC (Fig. 7E).

Discussion

In this study, we performed RNA‐Seq on a unique set of follow‐up biopsies and compared their mRNA expression to liver biopsies obtained at the time of PE and normal liver controls after excluding DEGs overlapping with other neonatal cholestatic disorders from the analysis. Results were analyzed with a bioinformatic approach using GO term and ssGSEA to assess the effects of SPE on global BA‐specific liver gene‐expression patterns, affected pathways, and relative abundance of cell types crucially involved with liver fibrosis.^(33‐35)

Earlier histopathological studies reported a decrease in neutrophil‐mediated inflammation following SPE.⁽⁵⁾ Our results confirm this finding at a molecular level insofar as four of the five most negatively affected GO pathways were directly related to neutrophil function. In addition, ssGSEA showed that the relative abundance of SAMacs and KCs was also decreased by SPE. However, the attenuation of inflammation was not restricted solely to neutrophils. Leukocyte migration, humoral immune response, and chemotaxis were among the most negatively affected GO pathways. Interestingly, the relative abundance of CD8 + T cells showed marked variation among the SPE samples. This may reflect variations between T helper 1/2–mediated responses playing a role in the rhesus rotavirus type A mouse model for BA.⁽³⁶⁾ Together, these results suggest that SPE reduces the intense immunological response seen at the time of PE.

Previous IHC assessments have confirmed that proliferation of ECM‐producing myofibroblasts and neo‐bile duct–forming cholangiocytes in DR persist after SPE.⁽⁵⁾ An earlier analysis of candidate genes demonstrated that molecular markers of active fibrogenesis are persistently expressed in BA livers after SPE, despite relief of cholestasis.⁽⁶⁾ Our results validate and extend these observations. As compared with NCs, 26 ECM‐related genes were highly differentially expressed at the time of PE. After SPE, 12 of these genes remained highly expressed, indicating their role in post‐SPE fibrosis. Together, these results suggest that progression of post‐SPE liver fibrosis is a dynamic, evolving process that differs from PE and involves significant changes in the molecular activity of ECM.⁽³⁷⁾

One theoretical approach to improve outcomes of BA is through tailored antifibrotic therapy administered to patients who clear their jaundice after PE. In this study, we identified an ECM molecular fingerprint, a set of 12 ECM genes constantly and highly differentially expressed from PE until 1 year after SPE. We hypothesize that these genes may be part of an “ECM fibrogenic niche,” and these genes may not only be a consequence of but also promote the fibrotic process.⁽³⁸⁾ LAMA3 and ITGB6 are particularly interesting in this respect.

As a part of the bile duct basement membrane, laminin‐322 has been proposed to play a crucial role in the DR niche sustaining the biliary phenotype of cholangiocytes.^(32,39) To our knowledge, expression of LAMA3 has not been evaluated previously in BA. We found overexpression of LAMA3 in cholangiocytes within areas of DR both at PE and after SPE. The role of LAMA3 in the formation of DR and fibrogenesis in BA will be the subject of future studies.

In addition to laminins, integrins play a role in liver fibrosis.⁽⁴⁰⁾ The unique β6 chain (encoded by ITGB6 gene) of integrin αvβ6 is exclusively expressed on activated cholangiocytes featuring fibrogenic progenitor cells during biliary fibrosis also in BA.⁽⁴¹⁾ Orally active small molecule inhibitors of integrin αvβ6 have been shown to potently attenuate fibrogenesis by inhibition of cholangiocytes.⁽⁴²⁾ We found overexpression of ITGB6 both at PE and after SPE, supporting the concept of integrin αvβ6 inhibitors as antifibrotic therapy after SPE. Another ECM‐related gene overexpressed both at PE and SPE was MMP7, a tissue‐remodeling proteinase expressed primarily by cholangiocytes. MMP7 has emerged as one of the most promising diagnostic markers for BA.^(21,43) IHC and candidate gene studies have previously shown overexpression of MMP7 after SPE,⁽²²⁾ offering validation of our RNA‐Seq results.

Our ssGSEA showed that the relative abundance of both neutrophils and tissue monocytes, such as SAMacs and KCs, decreases following SPE. These results, together with GO pathway analysis, suggest diminishing importance of immune cell–mediated fibrosis after SPE. Recently, ssGSEA by the Bezerra Lab associated both aHSCs and aPFs to BA fibrosis by showing their increased relative abundance in BA livers at the time of PE.⁽¹⁹⁾ In addition to confirming these results, we demonstrate differential changes in the abundance of these cells after SPE. Although the relative abundance of aHSCs was significantly decreased as a result of SPE, abundance of aPFs showed no similar significant decrease. This increase, together with the several persistently overexpressed genes mechanistically involved with biliary fibrosis and the observed progressive (albeit not statistically significant) in the relative abundance of cholangiocytes highlights their potential role as drivers of fibrosis after SPE and reinforces the previously reported role of cholangiocytes as promoters of the BA fibrosis,^(36,44,45) postulating that aPFs function as “myofibroblasts for cholangiocytes.”⁽⁴⁶⁾ Finally, the differences in ECM‐related gene expression and relative cell type abundancies between PE and SPE underscore the notion that “at PE” and “after SPE” are molecularly distinct environments, and therefore the factors driving DR and fibrosis differ.^(47,48)

Weaknesses of this study include the limited number of patients analyzed, the use of homogenized samples instead of more precise single‐cell sequencing methods, and the limited number of non‐age‐matched controls. To compensate for the low number of samples, we focused on the DEGs with more than a five‐fold difference as compared with normal liver controls. Notably, setting the threshold to 10‐fold did not change the results of affected pathways, as assessed by GO term analysis.

In summary, we report the effects of SPE on the BA liver transcriptome. We demonstrate attenuation of inflammation and evolving up‐regulation of ECM molecules mechanistically connected to biliary fibrosis as a result of SPE. Finally, persistently and highly overexpressed genes emerged as crucial candidates controlling liver fibrosis after SPE.

Acknowledgment

The authors thank the Biomedicum Functional Genomics Unit for helping with the RNA sequencing.