The BA transcriptome dataset GSE46960 was obtained from the public Gene Expression Omnibus (GEO) database; it included liver tissue from 64 cases of BA, 14 non-BA patients with intrahepatic cholestasis, and 7 normal samples (control). The differentially expressed genes of BA liver tissues (64 cases) and control group liver tissues (21 cases) were compared and analyzed with online GEO2R software. Histograms for the two groups were drawn according to the expression levels of PPARα, OCTN2, CPT1A, CYP7A1, MDR3, and NR1H4 (FXR). Differences were considered significant at P<0.05.

According to previous studies, free carnitine in BA patients can be used to diagnose BA, while PPARα and NR1H4 (FXR) are key genes in the PPAR signaling pathway and bile acid metabolism (7) with significantly reduced expression in BA. Gene set (h.all.v7.4.symbols.gmt) was used for gene set enrichment analysis (GSEA), which is a computational method that determines whether an a priori defined set of genes shows statistically significant, concordant differences between two biological states (15), which was performed for PPARα and NR1H4 (FXR) expression. Potentially involved signaling pathways were evaluated, and their molecular mechanisms in BA development were investigated.

Gene expression was divided into high and low groups according to median PPARα and NR1H4 (FXR) expression levels, and 1000 times set sequences were performed for each analysis. The expression levels of PPARα and NR1H4 (FXR) are given as the normalized enrichment score (NES), normalized significance level (nominal P value), and adjusted multiple hypothesis tests (false discovery rate [FDR] q value). Pair enrichment pathways for each phenotype are classified. Gene sets with NES≥1.0, nominal P value ≤ 0.05, and FDRq value ≤ 0.25 were identified as meaningful.

PPI networks in meaningful signaling pathways were obtained through the online analysis website STRING (https://string-db.org/) and exported in.TSV format. The obtained source file was imported into the open source software platform Cytoscape (https://cytoscape.org/) for visual analysis, and the plug-in cytoHubba was used for hub gene analysis to screen for the key gene SMAD family member 3 (SMAD3).

Extraction and clustering of SMAD3, PPARα, and NR1H4 expression was performed with the R package Consensus Cluster Plus (16). The sample was divided into two clusters.

CIBERSORT is a method of calculating cell composition based on an expression profile and can be used for almost any tissue. This algorithm was used to calculate the proportions of 22 immune cells in each BA patient. The sum of the 22 immune cell group fractions for each sample is 1 (17).

The degrees of infiltration of 28 immune cells were calculated based on the levels of gene expression in the 28 published immune cell genomes by applying the R Package GSVA Single Sample GSEA (ssGSEA) method. Data set GSE15235 is used for external validation.

This study was reviewed and approved by the Institutional Review

Board of the Shanghai Children’s Hospital. A total of 30 of BA、28 of choledochus cyst (CC) patients’ liver were used for verification. Liver samples of 10 BA and 10 CC patients were prepared as paraffin sections for immunohistochemistry, the remaining 38 were frozen and stored in liquid nitrogen for real-time PCR and western blot analyses.

Six paraffin sections approximately 4-μm thick were excised from each tissue block and used for immunohistochemical labeling with six antibodies purchased from Abcam (Cambridge, UK): recombinant anti-PPAR alpha [EPR21244], anti-solute carrier family 22 member 5 antibody, anti-CPT1A antibody, anti-CYP7A1 antibody, anti-ABCB4 antibody targeting the N-terminal, and human NR1H4 (FXR)/NR1H4 antibody.

Immunohistochemical staining was performed for PPARα, NR1H4 (FXR), OCTN2, CPT1A, CYP7A1, and MDR3 (ABCB4) using the above-mentioned antibodies. Since PPARα and NR1H4 (FXR) are expressed in the nucleus, immunohistochemical staining shows positive expression in the nuclei of hepatocytes. OCTN and MDR3 (ABCB4) were expressed on the cell membrane, and brownish-yellow staining of the hepatic cell membrane was considered as positive expression. Since CPT1A is expressed in the outer mitochondrial membrane and CYP7A1 is expressed in the endoplasmic reticulum and mitochondria, immunohistochemical staining revealed a positive brownish yellow color in the cytoplasm of the liver. All samples were processed with positive and negative controls.

Real-time PCR was performed to detect the expression of NR1H4 (FXR), PPARα, and SMAD3 in liver tissues from both groups of patients. Ten frozen BA specimens were used as the experimental group, and fourteen with choledochal cysts were used as the control group. Liver tissue was minced and mixed with 0.5 ml of TRIzol to completely dissolve the samples, which were centrifuged at 12,000 g for 5 min at 4°C. The supernatant was collected and added to 250 μl of chloroform, total RNA was extracted by the conventional TRIzol method, and the absorbance was measured in order to calculate the amounts of 1-μg RNA samples. We prepared a well-mixed 20 μl reverse transcription reaction system that was placed in a PCR machine for reverse transcription. We prepared a 10-μl real-time fluorescent quantitative PCR reaction system. The primer sequences for the target genes and the control glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene are listed in Table 1 . Relative expression of the target gene was calculated by the 2-ΔΔCt method. Th amplification conditions were 95°C for 20 min, 95°C for 30 s, 58°C for 30 s, 40 cycles, 72°C for 40s.

Western blotting was performed to detect the expression of NR1H4 (FXR), PPARα, and SMAD3 in liver tissue. Approximately 100 mg of frozen tissue was thawed in 1 ml of protein lysis buffer, and total protein was extracted by an ultrasonic method. Protein concentrations were determined using the bicinchoninic acid method. Equal samples were loaded for gel electrophoresis, the proteins were then transferred to membranes, which were blocked prior to adding the primary antibody for overnight incubation overnight at 4°C. The next day, the membranes were incubated with an appropriate secondary antibody for 1 hour at room temperature. Finally, the blot was incubated in a coloring solution, and the signal was developed after rinsing.

GraphPad Prism 9 (GraphPad Software, Inc, San Diego, CA, USA) was used to generate graphs, and SPSS 26.0 software (IBM Corp, Armonk, NY, USA) was used for data processing, measurement data (mean ± standard deviation), and immunohistochemical semi-quantitative analysis in the BA group and common bile duct cyst group. The integrated optical density levels of the two groups were non-normally distributed continuous measurement data. Differences were compared using rank-sum tests, and P<0.05 was considered statistically significant.

The GSE46960 dataset was downloaded from the GEO database and processed with GEO2R software. It included samples of between BA liver tissue (64 patients) and control liver tissue (14 non-BA patients with intrahepatic cholestasis, 7 normal liver patients). Gene expression levels were analyzed, and genes with adjusted P-values (Padj)<0.05 were considered as differentially expressed ( Figure 1A ). Histograms were drawn according to the expression levels of PPARα, OCTN2, CPT1A, CYP7A1, MDR3, and NR1H4 (FXR). PPARα, CYP7A1 and NR1H4 (FXR) were significantly differentially expressed between the two groups ( Figure 1B ).

The gene set h.all.v7.4.symbols.gmt was used for GSEA to assess response pathways that may be involved in PPARα and NR1H4 (FXR) signaling and BA development and progression. Gene expression data were divided into high and low expression groups according to PPARα and NR1H4 (FXR) levels in GSE46960, and the results were validated in 47 BA patients in the GSE15235 dataset. Results that were not statistically significant in both datasets were excluded, and the same proportions of results from both groups were retained. For PPARα and NR1H4 expression we used liver biopsy dataset GSE65359 (84 cases of hepatitis B), fetal and adult liver sample dataset GSE61276 (106 cases), and cholestasis gene expression excluding BA. Data were also selected from the GSE46960 (14 cases) and the GSE112790 liver cancer dataset (183 cases). The results showed that more BA-specific pathways were enriched in BA patients with low PPARα and NR1H4 expression. When PPARα expression is low, signaling is enhanced in the transforming growth factor-β (TGFβ), p53, and phosphoinositide 3-kinase (PI3K)-AKT-mammalian target of rapamycin (mTOR) pathways ( Figures 2A, B, E, F ). Low NR1H4 expression is associated with inflammatory responses, interleukin (IL)6/STAT3 signaling, early estrogen response, IL2 STAT5 signaling, TGFβ signaling, apoptosis, and angiogenesis ( Figures 2C, D, G, H ). TGFβ signaling was significant in both groups. We found that the same signaling pathway was not screened in the PPAR and FXR high expression groups. The TGFβ signaling pathway in the low expression group was significant in both groups ( Table 2 ).

The online analysis website STRING was used to generate a network of interactions between TGFβ pathway proteins and PPARα, NR1H4, and other proteins. The resulting PPI network had a total of 57 nodes, 345 edges, and PPI enrichment P<1.0e-16. Visual analysis was performed with Cytoscape ( Figures 3A, B ), and cytoHubba was used for hub gene analysis. The potential major hub gene was identified as SMAD3 ( Figure 3C and Table 3 ). We then screened all genes contributing to TGFβ signaling in these two data sets of biliary atresia and confirmed that SMAD3 played an important role in both data sets ( Figure 3D ). We speculated that SMAD3 activated the TGFβ pathway by interacting with PPARα and NR1H4, thereby promoting BA-related biliary fibrosis.

A correlation heat map was created for PPARα and NR1H4, which are the major genes of the TGFβ pathway ( Figure 4A ). Considering the important role of SMAD3 in the TGFβ pathway and its close relationship with PPARα and NR1H4, we included PPARα, NR1H4 (FXR), and SMAD3 in the follow-up analysis. We performed consistent clustering of 64 BA samples based on PPARα, NR1H4 (FXR), and SMAD3 expression matrices and split the samples into two clusters by median. The heat map showed low expression of SMAD3 and high expressions of PPARα and NR1H4 (FXR) in cluster 1 (n=36), and high expression of SMAD3 and low expression of PPARα and NR1H4 (FXR) in cluster 2 (n=28) ( Figure 4B ). Next we investigated the Spearman correlations between the three genes separately. We observed a positive correlation between PPARα and NR1H4 (R=0.48, P<0.001), and negative correlation between SMAD3 and PPARα (R=-0.41, P<0.001), as well as SMAD3 and NR1H4 (R=-0.67, P<0.001). Consistent clustering was also carried out for the samples in the validation group GSE15235 according to the three genes, and similar heat maps were obtained. The 47 samples were divided into groups 1 cluster 1 (n=27) and cluster 2 (n=20) ( Figure 4C ).

We performed CIBERSORT analysis and ssGSEA to better understand differences in immune function. CIBERSORT analysis showed that cluster 2 had high proportions of activated mast cells (MCs) and neutrophils ( Figure 5A ). These were not confirmed in the validation set ( Figure 5B ).

ssGSEA showed 28 immune cell subtypes. The results revealed that activated dendritic cells, CD56dim natural killer (NK) cells, monocytes, neutrophils, Type 17 T helper cells were more common in cluster 2 ( Figure 5C ). It was confirmed in the validation set ( Figure 5D ) that cluster 2 might have stronger immune infiltration than cluster 1. Suggesting that it may have stronger immune cell infiltration than cluster 1. (ns, not significant, *P<0.05, **P<0.01, ***P<0.001).

Samples from the control group showed complete hepatic lobule structure was with hepatocellular cords lined up in an orderly manner. In the BA group, the extrahepatic bile duct was dilated and occluded to varying degrees, and the hepatic lobule and fibrous tissue around the bile duct were hyperplastic. Hepatocyte degeneration and necrosis were observed. The hepatic lobule structure and hepatic plate arrangement were disturbed, and intracellular cholestasis was noted ( Figure 6A ).The samples did not meet the requirements of normality testing (P<0.05), so differences were compared with Mann–Whitney U tests. Quantitative comparisons of IHC results are shown in Figure 6B and Table 4 .

According to the results of PCR, we confirmed NR1H4 (FXR) expression in the BA group was 1.257±0.878 compared to 2.352±1.276 in the control group, which was statistically significant (P<0.01). PPARα expression levels were similar in BA and control groups (0.152±0.116 and 0.318±0.2, respectively; P<0.05). SMAD3 expression levels were significantly different between the BA and control groups (0.153±0.086 and 0.061±0.03, respectively; P<0.01) ( Figure 6C and Table 5 ).

The WB results showed that the mean levels of NR1H4 (FXR) expression in the BA and control groups were 0.555±0.376 and 1.407±0.482, respectively. The mean level of PPARα expression in the BA group was 0.339±0.185, compared to 0.914±0.366 in the control group. The mean levels of SMAD3 expression were 1.29±0.414 and 0.258±0.189 in the BA and control groups, respectively. The observed variables in each group were close to normal distribution (P>0.05) so t-tests were selected. The results showed that control group NR1H4 (FXR) expression was significantly higher than in the BA group (t=3.556, P<0.01). The same pattern was observed for PPARα expression (t=3.988, P<0.01). However, the expression level of SMAD3 was higher in the experiment(t=-4.708, P<0.01). Such results coincided with mRNA expression ( Figures 6D, E and Table 6 ).