The RNA-seq dataset GSE122340 was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) (15). The downloaded RNA-seq data was comprised of liver tissue from 171 BA patients and 7 normal controls (NC). The BA data was grouped as publisher previously reported, including the discovery cohort (n = 121) and the validation cohort (n = 50). Detailed dataset were available in Supplementary Table S1.

The HSPGs and genes related to liver fibrosis and ductular reaction were uploaded to the STRING database (https://www.string-db.org/) to analyze protein functional interaction relationships.

The heterogeneity of BA and NC tissues based on HSPGs expression levels was distinguished by PCA analysis. Results were visualized using R package “scatterplot3d” (version 0.3-42).

Immunocytes abundance of BA and NC liver transcriptomic data was calculated by Microenvironment Cell Population-counter (MCPcounter) (16). Using R package “MCPcounter” (version 1.2.0), the abundance of 10 kinds of immunocytes, including T cells, CD8 T cells, Cytotoxic lymphocytes, NK cells, B lineage, Monocytic lineage, Myeloid dendritic cells, Neutrophils, Endothelial cells and Fibroblasts, was estimated. Results were visualized using the R package “ggplot2” (version 3.3.5).

The difference of the expression levels of each HSPGs, fibrosis and ductular reaction markers and the abundance of immunocytes between BA and NC were analysed using the R function “wilcox.test”. The results were visualized using the R package “ggplot2” (version 3.3.5).

A non-negative matrix factorization (NMF) clustering was used by R package “NMF” (version 0.23.0) to identify BA subgroup based on HSPGs (17). The optimal value of k was obtained when the magnitude of the cophenetic correlation coefficient began to fall drastically in the discovery cohort. R function “consensusmap” was used to plot the heatmap of clustering with the value k from 2 to 10.

In this study, Gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of the differentially expressed genes (DEGs) were carried out using R package “clusterProfiler” (version 4.0.5) (18, 19). The visualization was handled via the R package “ggpubr” (version 0.4.0) and “GOplot” (version 1.0.2).

GSEA was conduucted using GSEA software (version 3.0), which was downloaded from http://software.broadinstitute.org/gsea/index.jsp. For pathway enrichment analysis, the c2.cp.kegg.v7.4.symbols.gmt was downloaded from the Molecular Signatures Database (http://www.gseamsigdb.org/gsea/downloads.jsp). FDR < 0.25 was considered statistically significant.

Machine learning algorithm least absolute shrinkage and selection operator (LASSO) regression was applied to screen hub HSPGs genes (20). R package “glmnet” (version 4.1-3) was used to conduct LASSO regression.

Correlation analyses were conducted using Spearman's method via the R function “cor.test”. The results were plotted using the R package “corrplot” (version 0.90), “circlize” (version 0.4.15) circlize and “ggplot2” (version 3.3.5).

The study enrolled patients with BA (n = 39) first treated at the Children's Hospital of Nanjing Medical University, Nanjing, China. The diagnosis of BA was based on intraoperative biliary tract exploration combined with histological features of liver biopsies. This study was approved by the Ethics Committee of the Children's Hospital of Nanjing Medical University and all participants’ parents signed informed consents.

Clinical characteristics were assessed for all patients, including sex, the age when underwent the Kasai portoenterostomy, level of direct bilirubin (DBil) before surgery and DBil level 3 months after surgery. Based on jaundice clearance (DBil level) at 3 months post-Kasai surgery, all patients were divided into 2 groups, jaundice group (DBil > 34.2 μmol/L) and jaundice-free group (DBil ≤ 34.2 μmol/L).

Total RNA was obtained from biopsy liver tissue using TRIzol reagent as described by the manufacturer (Invitrogen Life Technologies Co, USA). Quantitative RT-PCR was performed to determine the expression levels of SCD4 and GPC3. For mRNA detection, total RNAs (500 ng) were reverse transcribed using the reverse transcription kit (Takara, Tokyo, Japan) under 37°C for 15 min and 85°C for 30 s. RT-PCR was performed on the ABI Prism 7900HT (Applied Biosystems, USA). β-actin was used as an endogenous control. Forward and reverse primer's sequences of SCD4 and GPC3 were as follows: SDC4 former primer: GGACCTCCTAGAAGGCCGATA, reverse primer: AGGGCCGATCATGGAGTCTT; GPC3 former primer: CCTTTGAAATTGTTGTTCGCCA, reverse primer: CCTGGGTTCATTAGCTGGGTA. PCR was performed for 5 s at 95°C and for 30 s at 60°C for 40 cycles.

The specimens were fixed with 4% paraformaldehyde (PFA) and embedded in paraffin after tissue processing. After deparaffinization, rehydration, and antigen retrieval, 2-µm-thick sections were blocked with 10% normal donkey serum for 1 h at room temperature and stained with primary antibodies. The following primary antibodies were used in immunohistochemical staining: anti-SDC4 (1:500; Proteintech, Wuhan, China) and anti-GPC3 (1:500; Proteintech, Wuhan, China). Biotinylated anti-rabbit (Beyond, Shanghai, China) antibody was used as secondary antibodies according to the manufacturer's protocol. 3,3N-Diaminobenzidine Tertrahydrochloride (Beyond, Shanghai, China) was used as a chromogen according to the manufacturer's protocol. For double-immunofluorescence analysis, sections were pretreated in the same way as for IHC. The following primary antibodies were used: anti-SDC4 (1:50; mouse; Santa Cruz Biotechnology) and anti-KRT19 (1:500; rabbit; Abcam); anti-GPC3 (1:500; rabbit; Abcam) and anti-ACTA2 (1:500; mouse; Abcam). All slides were incubated with primary antibody overnight, which was diluted with 10% normal donkey serum after 1-h serum blocking. After DAPI (4’6-diamidino-2-phenylindole) staining, images were acquired using an Olympus BX51 microscope.

Continuous variables were presented as mean ± standard deviations (SD), and qualitative variable data were presented as a percentage with a 95% confidence interval. The Student's t-test was used for the comparison of measurement data between the 2 groups. For comparisons of enumeration data, Fisher's exact test was used. p-value < 0.05 was considered statistically significant. All analyses were executed by GraphPad Prism 8 software (GraphPad Software Inc., California, USA).

PCA analysis showed the heterogeneity of BA and NC tissues based on HSPGs expression levels (Figure 1A). The enrichment score of HSPGs was higher in BA using ssGSEA analysis (Figure 1B). SDC1, SDC3, AGRN, GPC2, GPC3, GPC4, HSPG2, COL18A1 and TGFBR3 were significantly upregulated in BA (Figure 1C and Table 1).

Transcriptional levels of liver fibrosis markers (ACTA2, COL1A1, COL1A2, COL3A1) and ductular reaction markers (EPCAM, KRT7, KRT19) were significantly upregulated in BA (Figure 2A).

HSPGs were shown in a PPI network with liver fibrosis and ductular reaction markers (Figure 2B). The expressions of AGRN, GPC1, GPC2, GPC4, HSPG2 were positively correlated with fibrosis and ductular reaction markers in BA, while SDC3 was only with fibrosis markers (Figure 2C). Although GPC3 expression was higher in BA than NC, it was negatively correlated with markers of these 2 features. The similar results were confirmed in the validation cohort (Supplementary Figure S1).

Based on the median of HSPG ssGSEA enrichment score, 121 BA patients in the dataset were divided into 2 group, HSPG_high (n = 61) and HSPG_low (n = 60). Via GSEA, several KEGG pathways associated with immune reaction were ranked in HSPG_high group, such as B_CELL_RECEPTOR_SIGNALING_PATHWAY, T_CELL_RECEPTOR_SIGNALING_PATHWAY, ECM_ RECEPTOR_INTERACTION and some singaling pathways (Figure 3A). The GSEA results indicated that HSPGs, reported as regulators of immune activation, might participate in the immunocyte infiltration in BA.

MCPcounter was applied to estimate the relative abundance of immunocytes in BA and NC microenvironment. The abundance of 10 types of immunocytes were shown in Figure 3B, including T cells, CD8 T cells, cytotoxic lymphocytes, B lineage, NK cells, monocytic lineage, myeloid dendritic cells, neutrophils, endothelial cells and fibroblasts. The abundance of T cells, CD8 T cells, B lineage, monocytic lineage, myeloid dendritic cells, endothelial cells and fibroblasts were upregulated in BA, while neutrophils was downregulated. The correlation analysis between immunocytes in BA patients showed the abundance of neutrophils was negatively correlated with some of other immunocytes, such as T cells, CD8 T cells, cytotoxic lymphocytes and fibroblasts (Figure 3C). Figure 3D showed the abundance of most immunocytes were positively correlated with the expression levels of fibrosis and ductular reaction markers in BA, while neutrophils was negatively correlated with these markers. Figure 3E showed that HSPGs were widely correlated with immunocytes in BA patients. Specifically, GPC3 was the only HSPG negatively correlated with the abundance of fibroblasts.

HSPGs were then applied to identify molecular subtypes using NMF consensus clustering. NMF heatmap and cophenetic correlation coefficients were shown in Supplementary Figure S2, S3. Optimal number of clusters were chosen as 3 for further analyses (Figure 4A). Then 121 BA patients in discovery cohort were divided into 3 clusters including C1 (n = 14), C2 (n = 63) and C3 (n = 44). Fibrosis and ductular reaction markers in C3 were lower than other clusters, indicating that patients in C3 might have less serious liver injury (Figure 4B–H). SDC2, GPC3 and GPC5 were significantly higher in C3 (Supplementary Figure S4A–N). The abundance of most immunocytes were lower in C3, such as T cells, cytotoxic lymphocytes, fibroblasts (Supplementary Figure S5A–J). GO enrichment analysis revealed that differential expressed genes of C1 were mainly enriched in immune response, cytokine response, cell adhesion, however, metabolic processes was mainly enriched in C3 (Figure 4I). Consistent with GO analyses, KEGG identified enrichment of immune related pathways in C1, such as PI3K-Akt, MAPK, Hippo, TNF signaling pathways, while metabolic pathways were enriched in C3 (Figure 4J). Interestingly, C3 subtype in BA seemed to exhibit relatively modest liver injury, suggesting that C3 patients might have better prognosis.

Heterogeneity of BA patients was observed based on the expression of HSPGs. HSPGs were screened via LASSO regression for the hub HSPGs of C1 and C3 subtype. For C1, SDC2, SDC3, SDC4, GPC1, GPC3, GPC4 and COL18A1 were calculated by LASSO (Figure 5A). For C3, SDC1, SDC2, GPC1, GPC2, GPC3, HSPG2 and COL18A1 were the selected genes (Figure 5B). The ROC curve of these genes showed that SDC4 was the suitable HSPGs for identifying C1 subgroup (AUC = 0.905) (Figure 5C). GPC3 was the suitable HSPGs for C3 subgroup (AUC = 0.962) (Figure 5D). The expression of SDC4 was significantly higher in C1 than other clusters, and GPC3 was higher in C3 (Figure 5E,F). The validation cohort showed GPC3 was negatively correlated with all liver fibrosis markers ACTA2, COL1A1, COL1A2, COL3A1 and ductular reaction markers KRT7 and KRT19, while SDC4 was positively correlated with ductular reaction index KRT19 (Supplementary Figure S6).

Based on jaundice clearance standard 3 months after Kasai surgery, 39 patients were divided into 2 groups, jaundice (n = 11) and jaundice-free (n = 28) (Figure 6A). The clinical characteristics of the patients in this study were similar in the 2 groups (Table 2). PCR results showed that SDC4 expression of liver tissues was significantly higher in jaundice group compared with jaundice-free group (Figure 6B). While GPC3 expression was significantly higher in jaundice-free group (Figure 6C).

To further validate the above observation, the expression of SDC4 and GPC3 were evaluated in BA liver tissue. SDC4 was strongly positive in portal areas and weakly positive in hepatocytes, while GPC3 had strong positivity in hepatocytes via immunohistochemical (IHC) staining (Figure 7A,B). Immunofluorescence assay showed the expression of SDC4 were mostly restricted with KRT19+ cells (Figure 7C). GPC3, however, was high-expressed in hepatocytes rather than in portal area, which indicates it may have function in hepatocytes (Figure 7D). Finally, correlation analysis showed the expression levels of GPC3 in 39 BA patients have significantly negative correlation with alanine aminotransferase (ALT) and aspartate aminotransferase (AST) level (Figure 7E,F).