Human liver cancer HepG2 cell line was obtained from the Korean Cell Line Bank (Seoul, Korea). HepG2 cells were cultured in Dulbecco's modified Eagle's medium (DMEM) (Gibco; Thermo Fisher Scientific, Inc.) containing 10% fetal bovine serum (FBS) (Gibco; Thermo Fisher Scientific, Inc.), and 100 U/ml penicillin/100 µg/ml streptomycin (P/S) (Gibco; Thermo Fisher Scientific, Inc.) at 37°C of 5% CO₂. Scutellarein (SCU) was purchased from Chengdu Biopurify Phytochemicals Ltd. (Chengdu, Sichuan, China).

HepG2 cells were seeded into 6-well plates at 1.5×10⁵ cells per well and treated with 100 µM of SCU for 48 h at 37°C of 5% CO₂. After incubation, total RNAs were extracted using TRIzol (Ambion, Thermo Scientific, Rockford, IL, USA). The concentration of RNA was determined using a spectrophotometer (BioTek, Winooski, VT, USA). Isolated total RNA was then subjected to sequencing to obtain expression data.

The mRNA sequencing was prepared by Theragenbio (Seongnam-si, Gyeonggi-do) using the following protocol (Table I). The libraries were prepared for 151bp paired-end sequencing using TruSeq stranded mRNA Sample Preparation Kit (Illumina, CA, USA). Utilizing oligo (dT) magnetic beads, mRNA molecules were specifically isolated and fragmented from 1 µg of total RNA. The fragmented mRNAs were synthesized as single-stranded cDNAs through random hexamer priming. Double-stranded cDNA was created by using this as a template for second-strand synthesis. After the sequential process of end repair, A-tailing, and adapter ligation, cDNA libraries were amplified with PCR (Polymerase Chain Reaction). The quality of these cDNA libraries was evaluated with the Agilent 2100 BioAnalyzer (Agilent, CA, USA). According to the manufacturer's library quantification methodology, they were measured using the KAPA library quantification kit (Kapa Biosystems, MA, USA). Following cluster amplification of denatured templates, sequencing was progressed as paired-end (2×151 bp) using Illumina NovaSeq6000 (Illumina, CA, USA).

Using the aligner STAR v.2.7.1a (40) and the ‘quantMode TranscriptomeSAM’ option for estimating transcriptome expression level, filtered reads were mapped to the reference genome associated with the species using ENCODE standard parameters (see to ‘Alignment’ of ‘Help’ section of the HTML report).

By using the ‘strandedness’ option in RSEM v.1.3.1 (41) and taking into account the reads' direction in relation to the library technique, gene expression was estimated. The ‘estimate-rspd’ option was used to increase measurement accuracy. The default settings were used for all other options. FPKM (Fragments Per Kilobase Million) and TPM (Transcripts Per Kilobase Million) values were generated to standardize sequencing depth among samples.

The R package TCC v.1.26.0 was used to identify DEG based on the projected read counts from the preceding stage (42). TCC program compares tag count data using effective normalizing techniques. With the use of the iterative DESeq2 (43)/edgeR (44) approach, normalization factors were computed. Using the p.adjust function of the R package with the default parameter settings, the q-value was determined based on the P-value. Based on the q-value criterion of less than 0.05 for correcting mistakes brought on by multiple testing, the DEG were found.

To predict drug and disease associations for candidate genes, our analysis was carried out in a functional database that was suggested (i.e., Drug Bank) by WEB-based GEne SeT AnaLysis (WebGestalt). Drug Bank and drug were chosen as the functional database and enrichment category, respectively, to predict the most important drugs for our target genes. The gene count was set at ≥ 5, and a P-value with a false discovery rate (FDR) adjusted of less than 0.05 was considered statistically significant. The WebGestalt derives significance results by clustering sets of genes using the Jaccard index as a criterion of similarity, prioritizing sets with significant P-value, and automatically identifying a ‘model’ or representative for each cluster (45).

Mienturnet (MicroRNA ENrichment TURned NETwork) is used to predict the miRNA target related to differentially regulated genes. 463 DEG was taken for prediction (both up and down-regulated). Out of 665 entries, we chose the top 10 entries based on the number of interactions.

GO database provides a set of hierarchical controlled vocabulary classified into 3 categories: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). For functional characterization of the DEG, a GO-based trend test was carried out using the R package called GOseq (46) through the Wallenius non-central hypergeometric distribution. Selected genes of P-value <0.05 following the test were regarded as statistically significant.

Molecular docking of Erlotinib and PTEN were performed using PyMol (The PyMOL Molecular Graphics System, Version 1.2r3pre, Schrödinger, LLC.) and USCF (47) chimera with the default parameters. Molecular docking of histone H1 and PTEN were performed using HDOCK server (http://hdock.phys.hust.edu.cn/) with the default parameters. Analysis for molecular interaction between histone H1 and PTEN was conducted through Discovery Studio 2018 (BIOVIA, San Diego, CA, USA).

The protein-protein interaction for up-regulated and down-regulated genes was performed using an online tool, Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) (Accessed on 24, Jun 2021). Up-regulated gene interaction showed about 61 nodes, whereas the expected number of nodes was 57 (P-value; 1×10¹⁶). Down-regulated gene interaction showed about 26 nodes, whereas the expected number of nodes was 2 (P-value; 7.47×10¹¹). Thus, the interaction was more significant than expected. In up-regulated genes, the significant pathways enriched were chosen based on KEGG, Wiki, and Reactome pathway analyses.

HepG2 cells were treated with indicated concentrations of SCU (0, 25, and 100) µM and lysed using radioimmunoprecipitation assay (RIPA) buffer (iNtRON Biotechnology, Seoul, Korea) containing phosphatase and protease inhibitor cocktail (Thermo Scientific, Rockford, IL, USA). Protein quantification were determined using a Pierce™ BCA assay (Thermo Fisher Scientific, Rockford, IL, USA). An equal quantity of protein (10 µg) from each sample was electrophoresed on (8–15)% SDS-polyacrylamide gels and transferred to a polyvinylidene difluoride (PVDF) membrane (ATTO Co., Ltd., Tokyo, Japan). Membranes were blocked with 5% (bovine serum albumin (BSA) in Tris-buffered saline containing 1% Tween 20 (TBS-T, pH 7.4) at room temperature for 1 h, and incubated overnight at 4°C with primary antibodies. The membranes were washed with TBS-T buffer for every 15 min in five times at room temperature, further they were incubated with 1:5,000 dilution of HRP-conjugated secondary antibody for 2–3 h at room temperature. The obtained proteins were detected by an electrochemiluminescence (ECL) detection system (Bio-Rad Laboratory, Hercules, CA, USA), and analyzed using the Image Lab 4.1 (Bio-Rad) program. The densitometry readings of the protein bands were normalized by comparison with the expression of β-actin as control, using the ImageJ software program (U.S. National Institutes of Health, Bethesda, MD, USA). Antibodies of histone H1-4 (Cat. no. 41328S), PTPRC (Cat. no. 72787S), and β-actin (Cat. no. 4970S) were purchased from Cell Signaling Technology (Danvers, MA, USA). Horseradish peroxidase (HRP)-conjugated secondary antibodies to anti-rabbit (cat. no. A120-101P) and anti-mouse (cat. no. A90-116P) were obtained from Bethyl Laboratories, Inc (Montgomery, USA).

Western blot experimental data were analyzed using GraphPad Prism version 8.0.2 (GraphPad Software). The findings were presented as the means ± standard deviation (SD) of triplicate samples. The unpaired Student's t-test was used to analyze the results, and a P-value of <0.05 was considered statistically significant.

In our previous study on scutellarein (SCU) treatment in human liver cancer HepG2 cells, we revealed that SCU could inhibit cell proliferation and metastasis through PTEN activation in HepG2 cells (30). Herein, SCU-treated group were compared against the SCU-untreated group which unraveled about total of 463 DEG (|log₂ FC|>1.0 and P-value <0.05), including 288 up-regulated and 175 down-regulated DEG (Table II). The DEG was represented as a volcano plot using an R-Bioconductor, which consists of a total of 60,676 variables (Fig. 1).

To identify the relationship between DEG and already available drugs we performed a drug-disease association analysis. After drug and disease association analyses (Fig. 2, Table III), we took the top 10 drugs which were significantly associated with DEG (FDR <0.05, gene count was set at ≥5). Among these drugs, the most effective cancer-specific drug is Erlotinib and which is a PTEN regulatory target drug (48). Furthermore, Erlotinib is associated with liver cancer, and it is used for the clinical treatment of liver cancer patients as a sole treatment or combination treatment with Sorafenib (49). To confirm the interaction of Erlotinib with PTEN, we performed drug and protein molecular docking (Fig. 3). UCSF Chimera software was used to validate the ligand-protein structures of these molecules. According to molecular docking studies, as shown in the image with the co-crystallized ligand, the two ligands occupied the active site. It was also discovered that many residues assisted in fitting the ligands into the binding pocket. The interacting amino acid residues involved in the bound complex of Erlotinib with PTEN were found to be PRO169, AGR172, ARG173, TYR176, TYR177, TYR180, PHE279, ILE280, PRO281, LEU318, THR319, LEU320, ASP324, and ASN323 (Table IV). The molecular dock scores suggest that the compound demonstrated PyMOL final intramolecular energy of −7.45 kcal/mol. In summary, SCU is associated with the drug Erlotinib, which is used as a clinical drug to treat PTEN-related diseases such as liver cancer. The development and use of peptides or inhibitors that target PTEN-related kinases, transcription factors, and cellular proteins could have a significant therapeutic benefit in the treatment of PTEN-related diseases (48). Therefore, we strongly suggest a strong relationship between DEG and PTEN.

Target identification of miRNA concerning DEG was performed with the use of an online tool Mienturnet (MicroRNA ENrichment TURned NETwork). Mienturnet predicts the target miRNA based on statistical analysis, network-based visualization, and analysis. Around 665 putative target miRNA were predicted concerning 463 DEG. The top 10 miRNA targets were selected based on the number of interactions and plotted in a graph. miR-335-5p and miR-26b-5p showed a high number of target interactions which is about 34 and 14 (Fig. 4).

Functional enrichment analysis is a method for locating gene or protein classes that are disproportionately represented in a big collection of genes or proteins that may be associated with anti-tumor properties. When these hub genes were subjected to enrichment analysis to identify significant gene ontology (GO) terms, a list of up-regulated and down-regulated genes involved in nucleosomes, cellular content, and molecular function were shown in Fig. 5. The highly expressed function in each category is listed in Table V. Additionally, Linker histone H1 (Histone H1) was involved in all obtained enriched functional analyses.

Fig. 6A and B show the protein-protein interaction (PPI) network of up-regulated and down-regulated genes. The histone has a higher degree of the up-regulated PPI network, including 61 nodes and 293 edges. In contrast, the down-regulated PPI network, including 26 nodes and 17 edges, shows no significant relations. The figure summarizes the network of anticipated connections for a specific protein group. The edges represent the anticipated functional connections based on seven types of evidence: fusion evidence, neighborhood evidence, concurrence evidence, experimental evidence, text mining evidence, database evidence, and co-expression evidence. In addition, the pathway enrichment analysis of the potential targets were identified using different database (Table VI). Table VI lists the results of the screening. The KEGG database were found to be necroptosis, drug metabolism, pentose and glucuronate interconversions, porphyrin and chlorophyll metabolism, and phagosome pathway. Wiki database was found to be histone modifications, codeine and morphine metabolism, tamoxifen metabolism, and translation factors pathway. Reactome database shown to be HDACs deacetylate histones, DNA methylation, PRC2 methylates histones and DNA, Transcriptional regulation by small RNAs, deubiquitination, apoptosis induced DNA fragmentation, and caspase activation via death receptors in the presence of ligand pathway.

The interactions among these three nodes and their first adjacent nodes were used to construct the sub-network (Fig. 7A-D). In our previous study, we have shown that SCU anti-tumor activity is regulated by PTEN (30), and transcriptomics analysis showed the role of histone is evident. To unravel the interaction between PTEN and histone H1, we performed molecular docking.

After examining the STRING analysis, the two crucial targets were validated by western analysis on SCU-treated HepG2 cells. These selected target proteins histone H1-4 (H1-4) and protein tyrosine phosphatase receptor type C (PTPRC) upon treatment with 100 µΜ of SCU in HepG2 cells. Fig. 8 shows that SCU up-regulated the expression of histone H1-4 (H1-4) and PTPRC respectively. These results further certify the role of SCU in regulating the expression of crucial targets in HepG2 cells.