Between April 2018 and May 2020, 38 patients with BA and 54 patients with choledochal cysts (CCs) were enrolled in the present study. All patients were diagnosed via laparoscopic bile duct exploration by the same surgical team at Shenzhen Children's Hospital (Shenzhen, China), and liver biopsy tissues were obtained at the time of surgery. The mean age of the patients in the BA group was 72.58±27.31 days, and the group included 15 male and 23 female patients (Table SII). The mean age of the patients in the CC group was 40.32±38.62 months, and the group included 14 male and 40 female patients (Table SIII). The liver tissues were immersed in RNA sample preservation solution (cat. no. R916331; Macklin, Inc.) and cryopreserved at -80˚C. The patients did not receive any treatment before surgery. The present study was approved by the Ethics Committee of Shenzhen Children's Hospital (approval no. SUMC2017-026). The parents of all the subjects provided written informed consent.

TRIzol^® reagent (Invitrogen; Thermo fisher Scientific, Inc.) was used to extract total RNA from BA liver tissues as well as CC tissues. The concentration and purity of RNA was detected using a Nanodrop-1000 (Thermo Fisher Scientific, Inc.) and Qubit RNA HS Assay Kit (cat. no. Q32852; Thermo Fisher Scientific, Inc.), then the yield and quality was evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc.) to test the integrity of the RNA. All RNA integrity numbers were >7 to ensure RNA quality.

Before the construction of the library, ribosomal RNA (rRNA) was removed using a Ribo-Zero Plus rRNA Depletion Kit (cat. no. 20037135; Illumina, Inc.). NEBNext^® Multiplex Small RNA Library Prep Set (cat. no. E7330S; Illumina, Inc.) was used to generate a sequencing library. To map the sequence to each sample, a barcode had to be added. Subsequently, the quality of the library was examined using a RNA high sensitivity chip on an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc.). VAHTS Library Quantification Kit for Illumina (cat. no. NQ101; Vazyme, Inc.) was used to accurately quantify the effective concentration of the library with an ABI StepOnePlus Real-Time PCR system (cat. no. 4376600; Applied Biosystems; Thermo Fisher Scientific, Inc.). The library was diluted to 20 pM as the final concentration. Then a sample cluster was performed on the cBot cluster generation system (cat. no. SY-312-2001; Illumina, Inc.) with TruSeq PE Cluster kit v3-cBot-HS (cat. no. 20015963; Illumina, Inc.), followed by paired-end sequencing of 125-bp reads with the Illumina HiSeq 2500 platform (Illumina, Inc.). All steps followed the manufacturer's protocols.

The accuracy of the sequencing results of the extracted RNA from liver BA and CC tissues was verified by filtering the sequencing data. Using Trimmomatic (23), the reads containing adapters were removed, then the low-quality sequences at the 5' and 3' ends were trimmed, and the reads containing >5% N bases were removed. This produced high-quality clean reads for all downstream analyses. The reference genome as well as gene annotation files were downloaded from the Ensembl genome browser (Ensembl GRCh37 release 110-July 2023; https://grch37.ensembl.org/index.html). The sequencing reads were mapped to the human genome using the HISAT2 software, and circRNAs were identified using ‘circRNA_finder’ analysis. Additionally, the expression of known miRNAs was compared with the precursor and mature miRNA sequences in miRbase (version 22) (24) using default parameters (25). The differentially expressed mRNAs, miRNAs and circRNAs were identified using the edgeR software package (26), with P<0.05 and log2 fold-change (FC)>1 as selection parameters. The pheatmap R package (https://cran.r-project.org/web/packages/pheatmap/pheatmap.pdf) was used to cluster the samples. CircRNA sequencing was performed by Vazyme Biotech Co., Ltd. qPCR was performed to verify the expression of circRNA in tissue samples and all primer sequences are presented in Table SI.

R software (version 4.0.2; https://cran.r-project.org/doc/contrib/Liu-R-refcard.pdf) was used to estimate DECs between BA and CC samples. GO (https://geneontology.org/) term enrichment analysis was performed, including molecular function, cellular component and biological process, along with KEGG (http://www.genome.ad.jp/kegg/) analysis. DECs were identified using the clusterProfiler (http://www.bioconductor.org/packages/release/bioc/html/clusterProfiler.html), org.Hs.eg.db (http://www.bioconductor.org/packages/release/data/annotation/html/org.Hs.eg.db.html), enrichplot (http://www.bioconductor.org/packages/release/bioc/html/enrichplot.html) and ggplot2 packages (https://cran.rstudio.com/bin/windows/contrib/4.2/ggplot2_3.4.2.zip) in Bioconductor, which is an R package used to perform GO functional and KEGG pathway enrichment analysis.

The expression of the circRNAs was validated using BA and CC tissues. Liver tissues were frozen in liquid nitrogen and then crushed into a homogenate. Total RNA was extracted using TRIzol^® reagent and transcribed into cDNA using an rtSTAR^™ First-Strand cDNA synthesis kit (cat. no. AS-FS-003-02; Arraystar, Inc.). Specific primers (presented in Table SI) were designed with Primer Premier 5.0 (Premier Biosoft), and synthesized by Vazyme Biotech Co., Ltd. Following the manufacturer's protocols, Arraystar SYBR^® Green Real-time qPCR Master Mix (Arraystar Inc.) was used for qPCR. The cycling conditions were 5 min at 95˚C for the initial denaturation period, then 15 sec at 95˚C for denaturation and 1 min at 60˚C for annealing and extension, repeated for 40 cycles. Expression levels were normalized to endogenous control (taqman endogenous controls FG, Human GAPDH; cat. no. 4352934E; Applied Biosystem; Thermo Fisher Scientific, Inc.), and the FC relative to circRNA expression levels in the CC group was calculated using the 2^-ΔΔCq method according to a previous study (27). All steps followed the manufacturers' protocols.

StarBase (v2.0) (28), TargetScan (https://www.targetscan.org/vert_80/) and miRanda software (https://cbio.mskcc.org/miRNA2003/miranda.html) were used to predict the downstream miRNAs of the circRNAs. The target mRNAs of candidate miRNAs were further analyzed using the miRDB (https://mirdb.org/custom.html), miRTarBase (https://mirtarbase.cuhk.edu.cn) and TargetScan databases. Subsequently, functional enrichment analysis of these mRNAs was carried out as a Venn diagram using R software. All databases were used according to default parameters.

In the present study, the online tool Search Tool for the Retrieval of Interacting Genes/Proteins (https://string-db.org/) was used to analyze the protein-protein interaction (PPI) of the predicted target genes, and Cytoscape software (version 3.4.0) (29,30) was used to construct the PPI network. Through node observation, the key nodes of the PPI network were examined. The Wilcoxon rank sum test was used to test for significant differences in topological properties between the BA and CC groups.

To evaluate the impact of differential gene expression on the disease status of BA, ROC curve analysis was used. This was conducted by plotting the ROC curve using gene expression data juxtaposed with the sample state (with or without BA), thereby enabling an assessment of gene expression accuracy (31). The pROC package in R software was used to generate these ROC curves (31). The entropy weight method was used to determine the entropy weight of each gene, following which the ROC curves for five genes with significant differences between the BA and CC groups were plotted. The accuracy of each biomarker was determined by the area under the curve (AUC) derived from the ROC curve analysis.

R software was used to integrate and analyze the data. Continuous variables are expressed as the mean ± standard deviation (at least 3 experimental repeats). An independent samples t-test was used to compare the continuous variables between BC and CC groups, as the samples were independent from each other. The figures were prepared using GraphPad Prism 8.0 (GraphPad Software; Dotmatics). P<0.05 was considered to indicate a statistically significant difference.

In the present study, the sequencing reads were first mapped to the human genome, and the circRNAs were then systematically identified and annotated using ‘circRNA_finder’ analysis. In total, 7,349 circRNAs were identified and three types of circRNAs were revealed, including exon, intergenic and intron circRNAs. Among them, most circRNAs (83.4%) were of the exon type, 6.4% were intergenic and 7.9% were intron type (a small percentage were not annotated). circRNA transcripts were distributed in the majority of chromosomes (Chr) (Fig. 1A). circRNAs from Chr1, Chr2, Chr3 and Chr13 accounted for 9.49, 7.48, 6.18 and 5.58%, respectively. These chromosomes corresponded to more than half of the RNAs of interest. In addition, the length of the circRNAs from these four choromosomes ranged from 152-9,637 bp, and the distribution frequency was 69.3% for circRNAs ranging from 152-1,000 bp and 16.4% for circRNAs >2,000 bp (Fig. 1B).

To identify the DECs associated with BA and CC, high-throughput analysis was performed on liver tissues from BA and CC cases (3 cases each). The R software package was used to analyze the differential expression, and a list of upregulated and downregulated DECs was obtained. According to the criteria of log2FC >2 and P<0.05, there were 78 DECs, including 16 upregulated circRNAs and 62 downregulated circRNAs (Fig. 2A). The top five upregulated genes and the top three downregulated genes (using the volcano map in Fig. 2B) were selected for subsequent RT-qPCR verification. Table I presents the names of all upregulated and downregulated circRNAs.

GO functional and KEGG pathway enrichment analysis were performed on the host genes of the 78 DECs using R software. The results of GO analysis indicated that the main biological processes of these host genes were ‘positive regulation of catabolic process’, ‘negative regulation of catabolic processes’, ‘regulation of microtubule motor activity’ and ‘cellular response to alcohol’. The primary cellular component category consisted of the categories ‘intracellular part’, ‘organelle part’, ‘plasma membrane region’ and ‘cytoplasm’. Finally, the main molecular functions included ‘enzyme binding’, ‘GTPase activating protein binding’, ‘glucocorticoid receptor binding’ and various enzyme activities, such as ‘transferase activity’ and ‘phosphotransferase activity, alcohol group as acceptor’ (Fig. 3A). The KEGG pathway analysis of these genes was primarily enriched in ‘pyruvate metabolism’, ‘ABC transporters’, ‘intestinal immune network for IgA production’, ‘viral myocarditis’, ‘leishmaniasis’, ‘Staphylococcus aureus infection’, ‘hematopoietic cell lineage’, ‘toxoplasmosis’, ‘cell adhesion molecules’, ‘systemic lupus erythematosus’ and ‘phagosome’ (Fig. 3B).

Liver tissue samples from 38 patients in the BA group and 54 patients in the CC group were analyzed using RT-qPCR. A total of five significantly upregulated circRNAs (hsa_circ_0006137, hsa_circ_0079422, hsa_circ_0007375, hsa_circ_0005597 and hsa_circ_0006961) and three significantly downregulated circRNAs (hsa_circ_0081171, hsa_circ_0084665 and hsa_circ_0075828) from the R software analysis were selected for RT-qPCR to verify the expression of these DECs. The RT-qPCR results demonstrated that the expression levels of hsa_circ_0006137, hsa_circ_0079422 and hsa_circ_0007375 were significantly increased (Fig. 4A-C), while the expression levels of hsa_circ_0081171 and hsa_circ_0084665 were significantly reduced (Fig. 4F and G) in patients with BA compared with the CC group. However, there was no significant difference in the expression levels of hsa_circ_0005597 and hsa_circ_0006961 between the two groups (Fig. 4D and E). Hsa_circ_0075828 also exhibited a significant increase in patients with BA compared with the CC group (Fig. 4H), contrary to the previous screening results. Therefore, five validated DECs were used for bioinformatics analysis.

An increasing number of studies have demonstrated that circRNAs can increase the expression levels of downstream genes by binding to miRNAs as molecular sponges (32,33). Therefore, 244 potential target miRNAs of hsa_circ_0006137, hsa_circ_0079422, hsa_circ_0007375, hsa_circ_0081171 and hsa_circ_0084665 were predicted through starBase (v2.0). According to competitive endogenous RNA (ceRNA) theory, there is a negative correlation between a circRNA and its target miRNAs (34). Therefore, through a literature search, seven miRNAs were selected as the target miRNAs of the circRNAs (hsa_circ_0006137/miR-26a-5p, hsa_circ_0006137/miR-145-5p, hsa_circ_0079422/miR-593-3p, hsa_circ_0007375/miR-1206, hsa_circ_0007375/miR-1208, hsa_circ_0081171/miR-18a-5p and hsa_circ_0084665/miR-22-5p) for further analysis. Subsequently, 430 target mRNAs were predicted to correspond to these seven miRNAs through the miRDB, miRTarBase and TargetScan databases (Fig. 5A).

To examine the potential functional role of the five circRNAs, GO and KEGG pathway enrichment analysis on the target mRNAs was carried out. As presented in Fig. 5B, these genes were significantly enriched in the forward transcriptional regulation of ‘protein serine/threonine kinase activity’, ‘RNA polymerase II promoter’, ‘DNA-binding transcription activator activity, RNA polymerase II-specific’, ‘DNA-binding transcription factor binding’ and ‘SMAD binding’. The associated pathways obtained using KEGG analysis were fewer, but ‘TGF-β signaling pathway’ and ‘EGFR tyrosine kinase inhibitor resistance’ were included in the enriched pathways (Fig. 5C). In summary, these functional analysis results suggested that the circRNA network may regulate the development of BA through the TGF-β and EGFR signaling pathways, which supports previous study results.

To further investigate the diagnostic potential of the aforementioned circRNAs, ROC analysis was used to evaluate the detection sensitivity and specificity. As presented in Fig. 6A-E, the AUC of hsa_circ_0006137, hsa_circ_0079422, hsa_circ_0007375, hsa_circ_0081171 and hsa_circ_0084665 in the differential diagnosis of BA compared with CC was >0.8, indicating that these circRNAs have a relatively high sensitivity and specificity for BA. These findings suggested that these circRNAs may serve as potential indicators for distinguishing BA from CC and could offer important insights for clinical research.