Paraffin-embedded tumor specimens from patients with HCC (n=308) were retrospectively collected from The Affiliated Hospital of Jining Medical University (Jining, China) between June 2013 and August 2019. The clinical and pathological details of the patients were retrospectively reviewed by accessing the hospital's existing electronic medical record system. Ultimately, 266/308 (86%) of patients with HCC with complete clinical information were selected for the subsequent clinical significance analysis. In addition, of the 308 patients, 66 were excluded from the survival analysis due to a lack of follow-up. And additional 40 HCC tissue samples were prospectively collected from The Affiliated Hospital of Jining Medical University between February 2020 and June 2021 for transcriptome sequencing. All samples were histopathologically and clinically confirmed as HCC tissues. The histology, diagnostic methods, and α-fetoprotein (AFP) levels were reviewed for each case of HCC identified to eliminate cases of metastatic cancer or other primary liver cancer types. Microscopically, the HCC cells were arranged in solid nests, trabeculae, acinar or pseudoglandular structures. Papillary structures were occasionally seen. Abundant sinusoidal capillaries could be seen between the tumor alveolus, with steatosis and eosinophilic bodies. Ethics approval was provided by The Ethics Committee of The Affiliated Hospital of Jining Medical University (Jining, China; ethical approval no. 2021C145). Each participant (>18 years of age) in the study provided written informed consent and signed informed consent forms.

The human hepatoblastoma cell line, HepG2, was purchased from ATCC via the Shanghai Huiying Biological Technology Co., Ltd. The cell line was authenticated by human STR profiling and confirmed to be mycoplasma free. HepG2 cells were cultured in DMEM (Nanjing KeyGen Biotech Co., Ltd.) containing 10% fetal bovine serum (FBS; Biological Industries) and 1% penicillin-streptomycin (Nanjing KeyGen Biotech Co., Ltd.), and incubated in a humidified 5% CO₂ incubator at 37°C.

The lentiviral vectors, GV513-ADAMTS16 (Ubi-MCS-CBh-gcGFP-IRES-puromycin) and GV358-BMP2 (Ubi-MCS-3FLAG-SV40-EGFP-IRES-puromycin), and the corresponding control lentiviruses containing the empty vector, CON335 and CON238, respectively, were purchased from Shanghai GeneChem Co., Ltd. Accession numbers for the two genes used in the present study are as follows: ADAMTS16, NM139056; BMP2, NM001200. For lentiviral infection, HepG2 cells were incubated with lentivirus at a MOI of 10 for 16 h, and the stable cell line was selected using 2 µg/ml puromycin for a week, followed by 2 µg/ml puromycin maintenance. The overexpression efficiency of ADAMTS16 and BMP2 were assessed through reverse transcription-quantitative PCR (RT-qPCR).

In total, 40 HCC tissue samples were divided into TB-pos and TB-neg groups by two independent pathologists on the basis of hematoxylin-eosin (HE) staining for transcriptome sequencing, which was conducted by Berry Genomics Co., Ltd. Briefly, using the NEBNext^® UltraTM RNA Library Prep Kit for Illumina^® (E7770S, New England Biolabs, Inc.), total RNA (≥1 µg) from each sample was used for generating sequencing libraries in-line with the manufacturer's instructions. Oligo (dT) magnetic beads were used to enrich and purify poly A-containing mRNA, followed by random interruption of mRNA into short fragments, which served as templates for random-primed cDNA synthesis performed using reverse transcriptase. Subsequently, purified double-stranded cDNA was subjected to end-repair, A-tailing and adapter ligation. Moreover, AMPure XP beads (Beckman Coulter, Beverly, USA) were used to purify cDNA library fragments for selecting the 250-300 bp cDNA fragments. Lastly, PCR amplification was performed by ABI Q3TMReal-Time PCR System, followed by purification of PCR products using the AMPure XP beads for obtaining the cDNA library.

Following library construction, the Qubit 2.0 Fluorometer (Thermo Fisher Scientific, Inc.) was used to quantify the library before diluting to 1.5 ng/µl. Then, the Agilent 2100 Bioanalyzer was used to assess the library insert size. When the expected insert size was obtained, qPCR was conducted to precisely quantify library effective concentration (>2 nM) for ensuing library quality. After the library quality was confirmed, the Illumina platform was used for sequencing to generate 150 bp paired end reads.

The edgeR software package (version 3.28.1; http://bioconductor.org/packages/release/bioc/html/edgeR.html) (22) in the R/Bioconductor environment (Release 3.6.1) was used for differential expression analysis. To control the false discovery rate, the resulting P-values were adjusted according to the Benjamini and Hochberg's approach (23). Genes with log2 (Fold Change) >1 and Q<0.05 were assigned as differentially expressed. To annotate the function of these DEGs, Gene ontology (GO) enrichment analysis of DEG sets were implemented in topGO (version 2.38.1; http://www.bioconductor.org/packages/release/bioc/html/topGO.html) in the R/Bioconductor environment (Release 3.6.1). GO terms with adjusted P<0.05 were considered as significantly enriched by DEGs.

Total RNA from the HepG2 cells was extracted using TRIzol reagent (Ambion; Thermo Fisher Scientific, Inc.). Reverse transcription was conducted using the HiScript III RT SuperMix for qPCR (+gDNA wiper) (Vazyme Biotech Co., Ltd.) according to the manufacturer's protocol. The reverse transcription reactions were performed at 42°C for 2 min, followed by 37°C for 15 min and 85°C for 5 sec. qPCR was performed using ChamQ Universal SYBR qPCR Master Mix (Vazyme Biotech Co., Ltd.) in a CFX Connect Real-time System (Bio-Rad Laboratories Inc.). GAPDH was employed as the internal reference and the 2^−∆∆Cq method was used for quantification (24). Primer sequences for GAPDH, ADAMTS16 and BMP2 were as follows: GAPDH, 5′-CTGACTTCAACAGCGACACC-3′ (forward) and 5′-TGCTGTAGCCAAATTCGTTGT-3′ (reverse); ADAMTS16, 5′-CCGGCCGGTACAAATTTTCG-3′ (forward) and 5′-AACAGCAGCTCCACAATCAGT-3′ (reverse); BMP2, 5′-AGAATGCAAGCAGGTGGGAA-3′ (forward) and 5′-TCTTGGTGCAAAGACCTGCT-3′ (reverse).

The inverse Matrigel invasion assays were performed as reported previously (25). To evaluate the TB capacity of HepG2 cells, 8-µm pore sized Transwell chambers (Corning, Inc.) were used. Briefly, the Transwell chambers were hydrated using serum-free DMEM medium for 30 min. Then, the Matrigel (BD Biosciences) was added to each chamber and solidified in an incubator at 37°C for 1 h. Afterwards, the chambers were inverted and the cells (2×10⁵/well) were seeded onto the upward facing underside of the chamber. Following cell attachment, serum-free DMEM medium was placed in the lower chamber and DMEM medium containing 20% FBS (serving as a chemoattractant) was added to the upper chamber. The cells were allowed to invade through the Matrigel for 3-5 days in a 37°C incubator. Finally, the 3D reconstruction of cell invasion was obtained through Zeiss 800 laser confocal microscopy. The ratio of the fluorescence values of GFP-positive budding cells/total fluorescence value of all GFP-positive cells × 100% was calculated using Image J software (version 1.8.0) to quantify the invaded cells (TB rate).

Spheroid-based sprouting assays were performed as previously described (26), with some optimization. HepG2 spheroids were obtained using the hanging drop method. HepG2 cells were suspended at a density of 1×10⁴ cells/45 µl and then seeded on the lid of 48-well culture plates for 2 days in an incubator at 37°C. Then, all spheroids were harvested and embedded in Matrigel for 2 days in a 37°C incubator. Finally, images were captured using a light microscope (OPTIKA).

In total, 308 formalin-fixed, paraffin-embedded tissue samples from patients with HCC were used to establish HCC tissue microarrays (TMAs) with a 2.0-mm diameter per core for IHC. In brief, TMAs were dewaxed, hydrated and antigen repaired by EDTA antigen repair buffer (PH 8.0), and the endogenous peroxidase activity was then blocked. TMAs were blocked in goat serum (OriGene Technologies, Inc.) at room temperature for 30 min to block non-specific staining. Primary antibodies were incubated with the TMAs overnight at 4°C, including ADAMTS16 (1:200; cat. no. DF9173; Affinity Biosciences) and BMP2 (1:200; cat. no. A0231; ABclonal Biotech Co., Ltd.). Subsequently, a horseradish peroxidase-labeled secondary antibody (KIT-5020; Maxim Co., Ltd., Fuzhou, China) was incubated with the TMAs for 30 min at room temperature. Finally, a chromogenic reaction was developed with DAB kit (DAB-0031 (20×); Fuzhou Maixin Biotech Co., Ltd.), followed by counterstaining with hematoxylin (G1120; Beijing Solarbio Science & Technology Co., Ltd.) for 10-30 sec at room temperature. All IHC staining analyses were assessed using a semi-quantitative scoring approach by two senior pathologists. IHC scores were calculated by multiplying staining intensity with the stained area, in which the staining intensity was scored as follows: 0 (no staining), 1 (light yellow staining), 2 (yellow-brown staining) and 3 (brown staining). The staining area was scored as follows: 1 (1–25%), 2 (26–50%), 3 (51–75%) and 4 (76–100%), according to the percentage of stained area in the field of vision. The median IHC score of ADAMTS16 or BMP2 was used as a cut-off to divide the samples into the low expression and high expression groups. ‘High’ was defined an IHC score higher than the cut-off value, and ‘low’ was defined as an IHC score lower than or equal the cut-off value.

SPSS 25.0 (IBM Corp.) was used for statistical analyses. Differences between two groups were analyzed using the Wilcoxon rank-sum test or unpaired student's t-test. The association of ADAMTS16 and BMP2 expressions with clinicopathological features was analyzed using the chi-square test or Fisher's exact test as appropriate. Kaplan-Meier (KM) curves were generated using the ‘survfit’ function in the survival package of R software (version 3.5.3), while significant differences in survival were compared using the log-rank test or Cramer-von Mises test as appropriate. Cramer-von Mises test was performed to generate the P-values when KM curves crossed over. P<0.05 was considered to indicate a statistically significant difference.

A total of 40 surgical HCC specimens were divided into two groups: TB-pos group (n=21) and TB-neg group (n=19), followed by transcriptome sequencing analysis. As shown in the volcano plot in Fig. 1A, 245 DEGs including 95 upregulated DEGs and 150 downregulated DEGs were obtained for the TB-pos group compared with the TB-neg group. GO enrichment analysis was subsequently performed for the 95 upregulated DEGs. The top 20 upregulated biological processes (BPs) are shown in Fig. 1B. The predominant BPs of the upregulated DEGs are associated with embryonic kidney development.

Results of GO (BP) enrichment analysis for downregulated genes are shown in Fig. 1C. It was noted that downregulated genes were mainly involved in immune-related processes, such as immune response, innate immune response and response to type I interferon. Certain studies have illustrated that the downregulated immune-related genes are associated with tumor immune escape (27–29).

The details of BP enrichment of upregulated genes in the GO analysis are provided in Table SI. In the present study, BPs of embryonic kidney development-related processes related to the upregulated genes of HCC with TB was made the focus as this is a novel area. Genes involved in the BPs associated with embryonic kidney development include ADAMTS16, BMP2, CALB1, FOXD1 and WT1. Furthermore, we found that there were 4 common upregulated genes in these embryonic kidney development-related processes, including ADAMTS16, BMP2, WT1 and FOXD1. ADAMTS16 and BMP2 were selected for further investigation. In addition, we noted that S100A9 which involved in neutrophil aggregation, was a poor prognosis-related gene. The survival analysis result of S100A9 based on TCGA data showed that the higher expression level of S100A9 was also linked with a worse prognosis of HCC patients (P=0.013) (Fig. S1). It was also reported in the literature that S100A9 expression could potentially serve as an independent prognostic marker for HCC (30).

To identify ADAMTS16 and BMP2 expression profiles within HCC tissues with different TB statuses, IHC staining was performed using TMAs of 308 paraffin-embedded HCC tissues. In addition, the DEGs of the budding cells and cancer center were compared. The statistical analysis results for the IHC staining scores of ADAMTS16 and BMP2 expression are summarized in Tables I and II. As shown in Table I, the staining scores of ADAMTS16 expression in the TB-pos HCC tissues were significantly higher than those in the TB-neg HCC tissues (P=0.005). However, the staining scores of ADAMTS16 expression demonstrated no significant difference between the tumor center and budding cells (P=0.174). As shown in Table II, the staining scores of BMP2 expression were significantly higher in TB-pos group compared with that in TB-neg group (P=0.015), and the significantly higher staining scores in budding cells than in tumor center (P=0.042) were observed. The representative results for IHC staining of ADAMTS16 and BMP2 expressions are illustrated in Fig. 2. These results demonstrated that upregulation of ADAMTS16 and BMP2 expression may be related to TB in HCC. Moreover, BMP2 expression was significantly increased in the budding cells compared with the tumor center.

To explore whether ADAMTS16 and BMP2 play a role in the TB of liver cancer, HepG2 cells with stable ADAMTS16-OE or BMP2-OE were constructed using lentiviral vectors. The overexpression efficiency was assessed through RT-qPCR assays (Fig. S2). Subsequently, the TB ability of HepG2 cells was evaluated using in vitro inverse Matrigel invasion and spheroid-based sprouting assays. As shown in Fig. 3A, the ratio of TB in ADAMTS16-OE or BMP2-OE HepG2 cells increased significantly relative to their respective controls. Similarly, the results of the spheroid-based sprouting assays demonstrated that ADAMTS16-OE or BMP2-OE in HepG2 cells resulted in more budding cells compared with the control cells (Fig. 3B). These results demonstrated that ADAMTS16 and BMP2 may be involved in the regulation of TB of liver cancer.

To analyze the association of ADAMTS16 and BMP2 expression levels with HCC clinicopathological characteristics, 266 HCC cases with complete available clinicopathological information for analysis were enrolled. These 266 cases were classified into the low- or high-expression groups by considering the respective median IHC scores for ADAMTS16 and BMP2 expression as the cut-off. The relationship between the clinicopathological characteristics of patients and gene expression was assessed using the chi-square test and Fisher's exact test. Consequently, ADAMTS16 expression was found to be significantly associated with necrosis (P=0.023) and cholestasis (P=0.011) (Table III). As shown in Table IV, BMP2 expression level had a significant association with the Barcelona Clinic Liver Cancer (BCLC) stage (B-C vs. 0-A; P=0.003) and the vessels encapsulating tumor cluster (VETC; P=0.014), which is one of the vessel types in HCC.

A total of 242 HCC cases with available follow-up data were enrolled for survival analysis. These 242 HCC cases were divided into two groups (low or high expression) based on the median IHC scores for ADAMTS16 or BMP2. As shown in Fig. 4A and B, the KM survival curves demonstrated that ADAMTS16 expression was statistically significantly associated with overall survival (OS; ADAMTS16, P=0.046) in patients with HCC. In addition, BMP2 expression was not associated with the OS of patients with HCC (BMP2, P=0.59). Subsequently, the survival curves of these two genes were downloaded from the GEPIA database (http://gepia.cancer-pku.cn/detail.php). As shown in Fig. 4C and D, the survival time of patients with HCC with high BMP2 expression significantly decreased compared with that of patients with low BMP2 expression (P=0.0081). ADAMTS16 expression did not differ significantly between the two groups (P=0.24). The GEO database (accession no. GSE76427; http://www.ncbi.nlm.nih.gov/geo/) was also utilized to perform survival analyses comparing expression levels of ADAMTS16 and BMP2 in HCC. The results demonstrated that neither ADAMTS16 (P=0.28) nor BMP2 (P=0.61) expression were associated with the OS time of patients with HCC (Fig. 4E and F). The potential reasons for this discrepancy were elucidated in the discussion.