Bioinformatic Analysis and In Vitro and In Vivo Experiments Reveal That Fibrillarin Participates in the Promotion of Lung Metastasis in Hepatocellular Carcinoma

Lung metastasis, the most frequent metastatic pattern in hepatocellular carcinoma, is an important contributor to poor prognosis. However, the mechanisms responsible for lung metastasis in hepatocellular carcinoma remain unknown. Aiming to explore these mechanisms, weighted gene coexpression network analysis (WGCNA) was firstly used to find hub genes related to lung metastasis. Then, we obtained 67 genes related to lung metastasis in hepatocellular carcinoma which were mainly related to ribosomal pathways and functions, and a protein interaction network analysis identified that fibrillarin (FBL) might be an important hub gene. Furthermore, we found that FBL is highly expressed in hepatocellular carcinoma and that its high expression increases the rate of lung metastasis and indicates a poor prognosis. Knockdown of FBL could significantly reduce proliferation and stemness as well as inhibiting the migration and invasion of hepatocellular carcinoma cells. Moreover, we found that FBL might be involved in the regulation of MYC and E2F pathways in hepatocellular carcinoma. Finally, we demonstrated that the knockdown of FBL could suppress hepatocellular carcinoma cell growth in vivo. In conclusion, ribosome-biogenesis-related proteins, especially Fibrillarin, play important roles in lung metastasis from hepatocellular carcinoma.


Introduction
With the characteristics of an insidious onset, rapid progress, and poor prognosis, hepatocellular carcinoma is one of the leading causes of cancer-associated deaths worldwide [1]. Although more and more treatments for hepatocellular carcinoma are being developed, the 5-year overall survival rate of hepatocellular carcinoma patients has not improved significantly in the past 10 years, and one of the most important reasons for this is that hepatocellular carcinoma is very prone to metastasis [2][3][4]. The metastatic sites of hepatocellular carcinoma include the lungs, bones, adrenal glands, brain, lymph nodes, and peritoneum, among which lungs are the most frequent site [5,6]. However, the mechanism underlying lung metastasis from hepatocellular carcinoma is still not very clear.
Weighted gene coexpression network analysis (WGCNA) is a bioinformatics analysis method in which a coexpression network is constructed across the entire microarray or RNA sequencing data. It aims to find gene modules with coexpressed genes and to explore the relationship between gene modules and phenotypes of interest as well as the core genes in the network [7]. The Cancer Genome Atlas (TCGA) is one of the most widely used databases in cancer research today, and from this, many researchers have used the WGCNA method to discover a variety of molecular mechanisms involved in tumor development [8,9]. Although many studies have used hepatocellular carcinoma data Bioengineering 2022, 9,396 2 of 24 from the TCGA to perform WGCNA analysis, no relevant research has investigated lung metastasis from hepatocellular carcinoma [10,11].
In our study, WGCNA was firstly used to find the hub genes associated with hepatocellular carcinoma and lung metastasis by comparing hepatocellular carcinoma cases with and without lung metastasis in the TCGA. Then, we used the Gene Expression Omnibus (GEO) database and our own cohort to verify the impact of the FBL gene on the prognosis. In addition, in vitro or in vivo experiments were used to verify the effects of the FBL gene on the proliferation, stemness, migration, and invasion of hepatocellular carcinoma cells ( Figure 1). Our findings may provide new ideas relevant to the exploration of the molecular mechanism responsible for lung metastasis in hepatocellular carcinoma. core genes in the network [7]. The Cancer Genome Atlas (TCGA) is one of the most widely used databases in cancer research today, and from this, many researchers have used the WGCNA method to discover a variety of molecular mechanisms involved in tumor development [8,9]. Although many studies have used hepatocellular carcinoma data from the TCGA to perform WGCNA analysis, no relevant research has investigated lung metastasis from hepatocellular carcinoma [10,11].
In our study, WGCNA was firstly used to find the hub genes associated with hepatocellular carcinoma and lung metastasis by comparing hepatocellular carcinoma cases with and without lung metastasis in the TCGA. Then, we used the Gene Expression Omnibus (GEO) database and our own cohort to verify the impact of the FBL gene on the prognosis. In addition, in vitro or in vivo experiments were used to verify the effects of the FBL gene on the proliferation, stemness, migration, and invasion of hepatocellular carcinoma cells (Figure 1). Our findings may provide new ideas relevant to the exploration of the molecular mechanism responsible for lung metastasis in hepatocellular carcinoma.

Data Acquisition and Screening
We downloaded gene expression profile data and clinical data collected from patients with hepatocellular carcinoma from the TCGA database (https://portal.gdc.cancer.gov/ projects/TCGA-LIHC, accessed on 27 February 2020). Transcriptome sequencing data of 370 hepatocellular carcinoma tissues and 47 normal liver tissues were obtained. We focused on the differences between samples of lung metastasis after hepatectomy and samples without recurrence or metastasis within 5 years after hepatectomy. According to the followup data, a total of 32 cases screened met the study's requirements (including 18 cases of lung metastasis after hepatectomy and 14 cases of no recurrence or metastasis within 5 years after hepatectomy). We also downloaded the GSE14520 dataset and corresponding clinical information from the GEO database (https://www.ncbi.nlm.nih.gov/geo/, accessed on Bioengineering 2022, 9,396 3 of 24 4 April 2020), which contained 247 hepatocellular carcinoma samples and 241 normal liver tissue samples adjacent to cancer. Besides, hepatocellular carcinoma proteomes containing 165 samples were obtained from Proteomic Data Commons PDC000198 (https: //pdc.cancer.gov/pdc/, accessed on 19 August 2020).

Construction of the Gene Coexpression Network
The WGCNA R package (version 1.69) was used to analyze the data [7]. Firstly, we obtained the top 5000 genes according to the median absolute deviation (MAD) size of each gene in the samples. Then, we clustered the samples, deleted outliers (TCGA.DD.A73G), and performed a gene coexpression network analysis on samples that met the criteria. A topological overlap matrix (TOM) was constructed with a suitable soft threshold β value (β = 9, scale-free R 2 = 0.85). Genes were classified according to the hybrid dynamic shearing tree method, the minimum number of genes in each module was set at 30, and finally the differences of each module were calculated. A heatmap plot was used to visualize the topological overlap matrix (TOM) for all genes.

Identification of Significant Relevant Modules and Hub Genes
The gene modules related to lung metastasis were determined by the correlations between the modules and lung metastasis. In our study, modules whose p value was <0.01 were considered statistically significant. For statistically significant modules, we used the "signedKME" function in the WGCNA package (version 1.69) to calculate the module membership (MM) and gene significance (GS) in the module. The genes with MM > 0.8 and GS > 0.4 were screened as genes related to lung metastasis. Further, we used the pheatmap package (version 1.0.12) to display the differences in the expression of lung-metastasis-related genes in hepatocellular carcinoma tissues and normal liver tissues as well as the differences between hepatocellular carcinoma tissues with lung metastasis and those without metastasis.

Gene Ontology Annotations and Kyoto Encyclopedia of Genes and Genomes Pathway Analysis of the Hub Genes
The Clusterprofiler package (version 3.16.0) was used to annotate and enrich the signal pathway of lung-metastasis-related genes using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database and the Gene Ontology (GO) database. This included biological processes (BPs), cellular components (CCs), and molecular functions (MFs). A p adjusted value of <0.05 was considered statistically significant [12].

Protein-Protein Interaction Network Analysis of the Hub Genes
The Search Tool for the Retrieval of Interacting Genes (STRING) database (version 11.0) (https://www.string-db.org, accessed on 16 April 2021) was used to analyze the interactions between the proteins of target genes to map the Protein-Protein Interaction network with a high level of confidence (minimum required interaction score = 0.700) [13]. The cytoHubba plug-in was used to screen the top 10 genes according to the Degree score, and the MCODE plug-in was used to screen the core subnetworks to determine the genes in the core area of the interaction network with the Cytoscape software.

Expression and Prognosis Analysis of the FBL in TCGA, GSE14520, and Our Cohort
We used the ggpubr package (version 0.4.0) in R to construct violin plots or box plots to show the difference between the expression of FBL in hepatocellular carcinoma tissues versus that in adjacent normal liver tissues. The Wilcoxon rank sum test was used for statistical testing, and p < 0.05 was judged to be statistically significant.
We used the survival package (version 3.1-12) and the survminer package (version 0.4.8) in R to calculate the cumulative survival time of patients through the Kaplan-Meier (KM) method, and the built-in Peto-Peto's modified survival estimate method in the survminer package was used for statistical testing. The overall survival time was defined as the time from the start of surgery to death due to any cause or until the last follow-up. Disease-free survival was defined as the time from the start of the operation until recurrence of the disease or death (for any reason).

Clinical Specimens
All HCC tissues and normal liver tissues were collected from HCC patients who underwent liver resection at the Sun Yat-Sen Memorial Hospital, Sun Yat-Sen University. All procedures were conducted with the approval of the Ethical Committee of the Sun Yat-Sen Memorial Hospital, Sun Yat-Sen University. Patient consent was obtained before the collection of tissue samples.

Cell Lines and Culturing Conditions
The human hepatocellular carcinoma cell lines Huh7 and PLC/PRF/5, and the immortalized liver cell line LO2 were purchased from the America Type Culture Collection. The cell lines were cultured in Dulbecco's modified Eagle's medium (BI, Biological Industries Israel Beit Haemek Ltd., Beit Haemek, Israel) supplemented with 10% fetal bovine serum (GIBCO, Life Technologies Corporation, New York, NY, USA), 100 µg/mL penicillin, and 100 µg/mL streptomycin at 37 • C in a humidified incubator containing 5% CO 2 .

Western Blotting
Cells and tissue samples were harvested in RIPA lysis buffer containing protease and phosphatase inhibitors to obtain protein extracts. A bicinchoninic acid (BCA) assay was performed to measure the protein concentrations. Then, equal amounts of the extracts were separated in an SDS-polyacrylamide gel, electrophoretically transferred to the polyvinylidene fluoride membrane (Millipore, Merck KGaA, Darmstadt, Germany), and incubated sequentially with primary and secondary antibodies after being blocked with 5% BSA. Finally, signals of proteins on the membrane were detected by chemiluminescence reagents (Millipore, Merck KGaA, Darmstadt, Germany) as per the manufacturer's instructions and imaged with an Imaging System. The relative densitometry was calculated with Image J (1.50d version). The antibodies used for the Western blot analysis were FBL antibody (ab166630, Abcam plc, Cambridge Biomedical Campus, United Kingdom; 1:1000) and GAPDH antibody (2118s, Cell Signaling Technology, Inc., Danvers, MA, USA; 1:2000). Except for human specimens, Western blotting was performed in three independent replicates to confirm the reproducibility of the experiments.

CCK8 Proliferation Assay
Cells transfected with siRNA for 48 h were seeded into a 96-well plate at a density of 1000 (Huh7) or 1500 (PLC/PRF/5) cells/well. After the planked cells grew adherent to the wall, CCK8 proliferation assay kits (Yeasen Biotechnology Co, Ltd., Pudong New Area, Shanghai, China) were used to test the cell proliferation at different time points (day 0,1,3,5) based on the manufacturer's instructions. The absorbance at a wavelength of 450 nm was measured on a TS Microplate Reader (Tecan Group Ltd., Seestrasse, Männedorf, Switzerland). These experiments were performed in triplicate.

EdU Proliferation Assay
A total of 8 × 10 3 huh7 cells or 1.5 × 10 4 PLC/PRF/5 cells transfected with siRNA for 72 h were seeded into a 96-well plate. After 24 h, as the normal cell growth cycle was restored, an appropriate EdU solution was added to the 96-well plate with the cells mentioned above. After incubation for 2 h, cells were fixed with 4% paraformaldehyde for 15 min, and the EdU incorporated into the DNA of cells was detected by fluorescent antibodies in accordance with the manufacturer's instructions. Finally, the proportion of proliferating cells was tested through a High Content Analysis on the ImageXpress Micro Confocal (Molecular Devices, San Jose, CA, USA). These experiments were performed in triplicate.

Colony Formation Assay
Huh7 or PLC/PRF/5 cells transfected with siRNA for 48 h were seeded into 6-well plates in complete growth medium. Initially, for Huh7, there were 1500 cells per well, and for PLC/PRF/5, there were 5000 cells per well. The culture medium was changed every 3 days. After 12 days, the cells were stained with 0.1% crystal violet after being fixed with 4% paraformaldehyde. The colony number was counted using ELISPOT Reader (AID GmbH, Penzberg, Germany). These experiments were performed in triplicate.

Sphere Formation Assay
A total of 500 cells transfected with siRNA for 72 h were seeded into low-adhesion 96-well plates in serum-free epithelial basal medium supplemented with B27 (50×) (GIBCO, Life Technologies Corporation, New York, NY, USA), 20 ng/mL EGF (PeproTech Inc., Cranbury, NJ, USA), and 20 ng/mL FGF (PeproTech Inc., Cranbury, NJ, USA) in each well. Fresh microsphere medium was regularly added for supplementary nutrients. Finally, we counted the spheres under a microscope after incubation for 7 days. These experiments were performed in triplicate.

Scratch Wound Healing Assay
When the cells transfected with siRNA for 72 h in a 6-well plate reached 95% confluence as a monolayer, we gently and slowly scratched the monolayer with a new 200 µL pipette tip across the center of the well to create a cross in each well. After scratching, serum-free medium was replenished into each well. Finally, we took photos of the scratches around the cross in each well under a microscope at 0 h (the time at which the scratching of the monolayer was completed) as well as 48 h later. The wound healing area of each scratch was quantitatively evaluated at different times using ImageJ. These experiments were performed in triplicate.

Transwell Migration and Invasion Assay
For the Transwell migration assay, 1 × 10 5 huh7 cells or 2 × 10 5 PLC/PRF/5 cells transfected with siRNA for 72 h were seeded with serum-free medium in the upper compartments of Transwell inserts (8 µm pore size, Corning, Glendale, AZ, USA), while medium containing 10% FBS was added to the lower bottom compartment. After 24 (huh7) or 36 h (PLC/PRF/5) of incubation, the inserts were taken out. Cells that remained in the upper chamber were gently removed with a cotton swab, and then cells on the lower chamber were fixed with 4% paraformaldehyde for 15 min and stained with 0.1% crystal violet for 20 min. Finally, we counted the number of cells on the lower chamber under a microscope. We randomly chose 5 different views and took the average number. These experiments were performed in triplicate.
For the Transwell invasion assay, 1 × 10 5 huh7 cells or 2 × 10 5 PLC/PRF/5 cells transfected with siRNA for 72 h were seeded with serum-free medium in the upper compartment of Transwell inserts (8 µm pore size, Corning, Glendale, AZ, USA) where the diluted Matrigel (diluted ratio = 1:9) was laid in advance, while medium containing 10% FBS was added to the lower bottom compartment. After 36 h of incubation, the inserts were taken out. Cells that remained on the upper chamber were gently removed with a cotton swab, and then cells on the lower chamber were fixed with 4% paraformaldehyde for 15 min and stained with 0.1% crystal violet for 20 min. Finally, the number of cells on the lower chamber was counted under a microscope. We randomly chose 3 different views and took the average number. These experiments were performed in triplicate.

Gene Set Enrichment Analysis
According to the median expression of FBL in hepatocellular carcinoma tissues in the TCGA and GSE14520 datasets, the samples were divided into high-expressing FBL samples and low-expressing FBL samples. Then, we imported the grouped sample data into the GSEA software (version 4.0.3) and selected "c7.all.v7.1.entrez.gmt" in the Molecular Signatures Database (MSig DB) on the GSEA website as the reference gene set to evaluate the impact of FBL on each reference set [14]. The default weighted enrichment statistical method was used to perform the Gene Set Enrichment Analysis (GSEA), and the number of random combinations was designed to be 1000. Gene sets that met the criteria of NES ≥ 1, NOM p-val < 0.05, and FDR q-val < 0.25 were defined as statistically significant enriched gene sets.

Establishment of the FBL Knockdown Stable Transfectant in Huh7
Huh7 cells were transfected with the FBL knockdown lentivirus and the corresponding negative control lentivirus (IGE Biotechnology Ltd., Guangzhou International Bio Island, Guangzhou, China) at a confluence of 50-70%, respectively. Multiplicity of infection (MOI) was set at 10. The fresh medium was replaced after transfection for 8 h. After 3 days, culture medium with 4 µg/mL puromycin was used to select successfully transfected cells. Finally, Western blotting was used to confirmed the efficacy of the infection 7 days later.

Xenograft Mouse Model
The 10 male nude mice (6-8 weeks old) were randomly divided into two groups. Then, 5 × 10 6 Huh7 cells stably transfected with the negative control lentivirus or the FBL knockdown lentivirus were injected subcutaneously into mice with 100 µL PBS-Matrigel mixture (PBS: Matrigel = 1:1) respectively. Tumor sizes and volumes were measured regularly, and the tumor volumes were calculated as follows: tumor volume = width/2 × width × length. After 4-5 weeks, the mice were sacrificed, and the subcutaneous tumors were excised, imaged, and weighted. We strictly followed the guidelines for laboratory animal care of Sun Yat-Sen University, and the following conditions were regarded as the humane end point of this experiment: the maximum diameter of the subcutaneous tumor in nude mice exceeded 1.5 cm, or the weight loss of nude mice was greater than 10%, or there were obvious endangered manifestations in mice (such as significantly reduced food intake, significantly reduced activity, obvious respiratory depression, etc.).

Statistical Analysis
In this study, SPSS 24.0 was used for the statistical analysis, and the GraphPad Prism 7.0 software, Adobe Illustrator software, and R language were used for mapping. All in vitro experiments were performed in triplicate to confirm the reproducibility of the results obtained. All the experimental data are presented as the mean ± standard error (S.E.M.). Comparisons between two or more groups were statistically analyzed by the Student's t test, Wilcoxon rank sum test, paired t-test, or random block design analysis of variance depending on the corresponding data type (Student's t test: Western blotting analysis of FBL protein expression in hepatocellular carcinoma and adjacent liver tissues, tumor growth curve chart comparison, weight of gross tumors comparison; Wilcoxon rank sum test: FBL expression in hepatocellular carcinoma and adjacent liver tissues in TCGA and GEO; paired t-test: RT-qPCR analysis of FBL mRNA expression in hepatocellular carcinoma and adjacent liver tissues; random block design analysis: CCK8 proliferation assay, EdU Proliferation Assay, colony formation assay, sphere formation assay, scratch wound healing assay, Transwell migration and invasion assay between the control and two experimental groups). Four-grid table data were statistically analyzed by the χ 2 test. p < 0.05 was considered statistically significant. Log-rank test was used for overall survival and disease-free survival without curved intersection. Peto-Peto's modified survival estimate test was used for overall survival and disease-free survival with curved intersection. Pearson correlation test was performed to analyze the correlation amongst FBL and other genes.

Coexpression Network Construction and Module Identification
Hierarchical clustering with the Euclidean distance was used to determine whether there were sample outliers. As shown in Figure 2A, one sample (TCGA.DD.A73G) was found to deviate significantly from the group, so it was removed from the analysis. We constructed a weighted gene coexpression network for the remaining samples. As shown in Figure 2B, the optimal soft threshold β value was 9, which could make the gene distribution conform to the scale-free network. According to the hybrid hierarchical clustering and dynamic cut tree method, 17 modules were identified with a different color for each module ( Figure 2C). Afterward, we constructed a heatmap of the network to visualize the topological overlap matrix (TOM) among all the genes and analyze the interactive relations among the 17 modules ( Figure 2D).
Hierarchical clustering with the Euclidean distance was used to determine whether there were sample outliers. As shown in Figure 2A, one sample (TCGA.DD.A73G) was found to deviate significantly from the group, so it was removed from the analysis. We constructed a weighted gene coexpression network for the remaining samples. As shown in Figure 2B, the optimal soft threshold β value was 9, which could make the gene distribution conform to the scale-free network. According to the hybrid hierarchical clustering and dynamic cut tree method, 17 modules were identified with a different color for each module ( Figure 2C). Afterward, we constructed a heatmap of the network to visualize the topological overlap matrix (TOM) among all the genes and analyze the interactive relations among the 17 modules ( Figure 2D).

Identification of Significant Modules and Hub Genes
A correlation analysis was performed for the 17 modules obtained and the traits of interest (lung metastasis) to look for significant associations. As shown in Figure 3A, the light cyan module and the turquoise module with p value < 0.01 were identified as significant modules. The correlation coefficients between the two modules mentioned above and lung metastasis were 0.5 and 0.48, respectively. These values were higher than the absolute values of the correlation coefficients of other modules, indicating that these two modules had the closest relationship with lung metastasis. We further analyzed the internal genes of the two modules mentioned above. As shown in Figure 3B,C, the correlation coefficient between GS and MM in the light cyan module was 0.43, and that for the turquoise module was 0.55, indicating that the two modules were suitable for identifying genes related to lung metastasis. The genes with GS >0.4 and MM >0.8 in significant modules were identified as hub genes. In the light cyan module, there were Bioengineering 2022, 9, 396 9 of 24 8 hub genes; in the turquoise module, there were 61 hub genes. After taking the union of the two, we obtained 67 hub genes related to lung metastasis (Table S1). Next, we used the transcriptome sequencing data in the TCGA to show the difference in expression of these 67 genes related to lung metastasis in hepatocellular carcinoma tissues and normal liver tissues as well as that between hepatocellular carcinoma tissues with and without lung metastasis. As shown in Figure 3D,E, compared with normal liver tissues, 67 genes related to lung metastasis were highly expressed in hepatocellular carcinoma tissues (Table S2), and compared with hepatocellular carcinoma tissues without metastasis, most of 67 genes related to lung metastasis were also upregulated in hepatocellular carcinoma tissues with lung metastasis in TCGA dataset (Table S3).

Gene Ontology Annotations and Kyoto Encyclopedia of Genes and Genomes Pathway Analysis of 67 Hub Genes Related to Lung Metastasis
The KEGG pathway enrichment analysis was performed on 67 hub genes related to lung metastasis. As shown in Figure 4A, a total of four pathways were found to be enriched, and these were associated with the ribosome and RNA polymerase. The GO analysis results showed that the 67 hub genes were mainly distributed among items related to the ribosome, mitochondria, and RNA transcription for cellular components (CCs), such as the ribosomal subunit, mitochondrial protein complex, and RNA polymerase II-core complex ( Figure 4B). They were mainly concentrated in items related to the ribosome composition and RNA polymerase activity for molecular functions (MFs), such as the structural constituent of the ribosome, RNA polymerase II activity, and snoRNA binding ( Figure 4C). They were enriched in items related to ribosome synthesis and protein translation processes for biological processes (BP), such as SRP-dependent co-translational protein targeting of the membrane, ribosome biogenesis, and rRNA metabolic processes ( Figure 4D). This indicates that ribosome and rRNA-related pathways might play important roles in lung metastasis associated with hepatocellular carcinoma.

Protein-Protein Interaction Network Analysis of 67 Hub Genes Related to Lung Metastasis
After entering the 67 genes related to lung metastasis in the STRING database, 67 nodes and 246 edge PPI networks were obtained ( Figure 5A). In order to simplify the PPI network to find the most critical core genes, the cytoHubba plug-in from the Cytoscape software was used to obtain the degree score of each gene. The top 10 genes were ribosomal proteincoding genes or genes related to ribosome production, including NHP2L1, FBL, RPS15, RPS16, RPL18, RPL35A, RPLP2, RPL27A, RPL24, and RPS11 ( Figure 5B). In addition, we obtained the core subnetwork through the MCODE plug-in of the Cytoscape software. As shown in Figure 4C, we found that FBL was the seed gene of the core subnetworks.

FBL Is Highly Expressed in Hepatocellular Carcinoma, and Its High Expression Is Closely Related to Poor Prognosis and Lung Metastasis in Hepatocellular Carcinoma Patients
According to the results of the Protein-Protein Interaction network analysis, we selected the ribosome-related gene FBL, which might play an important role for in-depth studies. In order to determine whether FBL is highly expressed in hepatocellular carcinoma, we analyzed FBL expression in samples of hepatocellular carcinoma tissue and normal liver tissue adjacent to the cancer in the TCGA and GSE14520. As shown in Figure 6A,C, the concentration of FBL mRNA in hepatocellular carcinoma tissues was significantly higher than that in normal liver tissues (p < 0.001). To improve the comparability, we further carried out a paired analysis of the tumor tissues and adjacent normal tissues from the above datasets, and the differences were found to be statistically significant ( Figure 6B,D). Moreover, in order to analyze the correlation between FBL expression and clinical prognosis in patients with hepatocellular carcinoma, we divided the two cohorts (TCGA and GSE14520) into a high-expressing FBL group and a low-expressing FBL group according to the median of FBL expression. The results of the survival analysis showed that high FBL expression was positively correlated with a shorter overall survival (OS) (Figure 6E,G) and shorter disease-free survival (DFS) (Figure 6F,H) in hepatocellular carcinoma patients, and the differences were found to be statistically significant.  composition and RNA polymerase activity for molecular functions (MFs), such as the structural constituent of the ribosome, RNA polymerase II activity, and snoRNA binding ( Figure 4C). They were enriched in items related to ribosome synthesis and protein translation processes for biological processes (BP), such as SRP-dependent co-translational protein targeting of the membrane, ribosome biogenesis, and rRNA metabolic processes (Figure 4D). This indicates that ribosome and rRNA-related pathways might play important roles in lung metastasis associated with hepatocellular carcinoma.  In order to verify the above findings, we first detected the expression of the FBL protein in five pairs of matched hepatocellular carcinoma tissues and adjacent normal liver tissues by Western blot analysis ( Figure 7A). The expression of FBL mRNA was detected by RT-qPCR in 16 pairs of matched hepatocellular carcinoma tissues and adjacent normal liver tissues ( Figure 7B). The results showed that, compared with adjacent normal liver tissues, the expression of FBL in terms of the mRNA level and protein level in hepatocellular carcinoma tissues was upregulated, and the differences were statistically significant. Meanwhile, immunohistochemistry was used to detect FBL expression in 229 patients with hepatocellular carcinoma in our hospital. According to the immunohistochemistry results, we divided hepatocellular carcinoma patients into a low FBL expression group ( Figure 7C, Negative and Weak) and a high FBL expression group ( Figure 7C, Moderate and Strong), and the clinical data were sorted and analyzed according to the grouping (Table 1). We found that, compared with hepatocellular carcinoma patients with low FBL expression, hepatocellular carcinoma patients with high FBL expression had a higher probability of experiencing lung metastasis ( Figure 7D). In addition, consistent with the results of the survival analysis in the TCGA and GSE14520, the overall survival of patients with high FBL expression was shorter ( Figure 7E), and disease-free survival was also shorter ( Figure 7F), and the differences were statistically significant. Moreover, as shown in Table 1, the high expression of FBL in hepatocellular carcinoma patients was closely related to the serum alpha-fetoprotein level, the number of tumors, the pathological grade, vascular invasion, and TNM staging.

Protein-Protein Interaction Network Analysis of 67 Hub Genes Related to Lung Metastasis
After entering the 67 genes related to lung metastasis in the STRING database, 67 nodes and 246 edge PPI networks were obtained ( Figure 5A). In order to simplify the PPI network to find the most critical core genes, the cytoHubba plug-in from the Cytoscape software was used to obtain the degree score of each gene. The top 10 genes were ribosomal protein-coding genes or genes related to ribosome production, including NHP2L1, FBL, RPS15, RPS16, RPL18, RPL35A, RPLP2, RPL27A, RPL24, and RPS11 ( Figure 5B). In addition, we obtained the core subnetwork through the MCODE plug-in of the Cytoscape software. As shown in Figure 4C, we found that FBL was the seed gene of the core subnetworks.  ure 6B,D). Moreover, in order to analyze the correlation between FBL expression and clin-ical prognosis in patients with hepatocellular carcinoma, we divided the two cohorts (TCGA and GSE14520) into a high-expressing FBL group and a low-expressing FBL group according to the median of FBL expression. The results of the survival analysis showed that high FBL expression was positively correlated with a shorter overall survival (OS) (Figure 6E,G) and shorter disease-free survival (DFS) (Figure 6F,H) in hepatocellular carcinoma patients, and the differences were found to be statistically significant. In order to verify the above findings, we first detected the expression of the FBL protein in five pairs of matched hepatocellular carcinoma tissues and adjacent normal liver tissues by Western blot analysis ( Figure 7A). The expression of FBL mRNA was detected by RT-qPCR in 16 pairs of matched hepatocellular carcinoma tissues and adjacent normal liver tissues ( Figure 7B). The results showed that, compared with adjacent normal liver

Knockdown of FBL Suppresses Proliferation and Stemness in HCC Cell Lines
As shown in Figure 8A, the expression of FBL in HCC cell lines (Huh7 and PLC/PRF/5) were higher than that in immortalized liver cell line (LO2). The huh7 and PLC/PRF/5 cell lines were chosen to establish the FBL knockdown models through siRNA in our study. The knockdown efficiency of FBL was verified by Western blotting analysis, indicating that the silencing efficiency of two siRNA sequences was >70%. The siFBL#2 sequence was found to have a higher silencing efficiency than the siFBL#1 sequence ( Figure 8B). To determine the effect of FBL knockdown in the proliferation of cells mentioned above, Cell-counting kit-8 (CCK-8) assays were performed, and the results showed that the silencing of FBL inhibited the proliferation of HCC cells compared with the control ( Figure 8C). Furthermore, EdU assays also revealed that the silencing of FBL could inhibit the proliferation of HCC cells ( Figure 8D). As we can see in Figure 8C, DNA synthesis was significantly reduced in FBL-silenced groups (red dots represented cells undergoing DNA synthesis). It was reported that FBL played an important role in promoting pluripotency in embryonic stem cells, but the role of FBL in the stemness of HCC cells is unknown [15]. Thus, we performed sphere formation and colony formation assays to compare the level of stemness between the FBL-silenced group and the control group. The results indicate that the silencing of FBL could remarkably reduce sphere formation ( Figure 8E) as well as colony formation ( Figure 8F) in HCC cell lines.

Knockdown of FBL Suppresses Proliferation and Stemness in HCC Cell Lines
As shown in Figure 8A, the expression of FBL in HCC cell lines (Huh7 and PLC/PRF/5) were higher than that in immortalized liver cell line (LO2). The huh7 and PLC/PRF/5 cell lines were chosen to establish the FBL knockdown models through siRNA in our study. The knockdown efficiency of FBL was verified by Western blotting analysis,  HBV, hepatitis B virus; AFP, Alpha-fetoprotein; TNM, tumor-node-metastasis. * p < 0.05, ** p < 0.01, *** p < 0.001. χ 2 test is used here.

Knockdown of FBL Suppresses Migration and Invasion in HCC Cell Lines
Migration and invasion capacities are characteristics of highly malignant cells and are important causes of tumor metastasis. The results mentioned above show that FBL is associated with metastasis, so we were interested to determine whether FBL could affect the migration and invasion of hepatoma cells. To investigate the effect of FBL in migration, we conducted both wound healing assays and Transwell migration assays. The results showed that FBL silencing evidently slowed the speed of wound healing ( Figure 9A,B) and reduced the number of migratory cells in the Transwell migration assay ( Figure 9C). In addition, Transwell invasion assays were carried out to test the invasion capacity after FBL silencing, and we found that FBL knockdown greatly weakened the capacity of invasion in HCC cell lines ( Figure 9D). the proliferation of HCC cells ( Figure 8D). As we can see in Figure 8C, DNA synthesis was significantly reduced in FBL-silenced groups (red dots represented cells undergoing DNA synthesis). It was reported that FBL played an important role in promoting pluripotency in embryonic stem cells, but the role of FBL in the stemness of HCC cells is unknown [15]. Thus, we performed sphere formation and colony formation assays to compare the level of stemness between the FBL-silenced group and the control group. The results indicate that the silencing of FBL could remarkably reduce sphere formation ( Figure 8E) as well as colony formation ( Figure 8F) in HCC cell lines.

FBL Might Be Involved in Multiple Signaling Pathways in Hepatocellular Carcinoma
The hepatocellular carcinoma samples in the TCGA and GSE14520 databases were classified into FBL high expression or FBL low expression groups according to the median FBL expression. We performed a gene set enrichment analysis (GSEA) with the MSigDB hallmark gene set on the FBL high expression and low expression groups from the TCGA and GSE14520 datasets. As shown in Figure 10A, there were five gene sets in the FBL high expression group that were simultaneously enriched in the above two datasets, namely the MYC_TARGETS_V1, MYC_TARGETS_V2, G2M_CHECKPOINT, E2F_TARGETS, and DNA_REPAIR gene sets. According to previous research, all of the gene sets mentioned above are closely associated with cell proliferation, and the MYC gene and E2F gene family were also found to be closely related to the stemness and metastasis of tumors, which suggests that FBL might participate in the regulation of these pathways and affect the proliferation, stemness, migration, and invasion of hepatocellular carcinoma. In addition, as shown in Figure 10B, three gene sets in the FBL low expression group were simultaneously enriched in the above two datasets. These three gene sets were BILE_ACID_METABOLISM, FATTY_ACID_METABOLISM, and ADIPOGENESIS, which suggests that FBL might be related to bile acid metabolism and fatty acid metabolism in hepatocellular carcinoma.

Knockdown of FBL Suppresses Hepatocellular Carcinoma Cell Growth In Vivo
It has been reported that the knockdown of FBL is a negative factor for tumor growth in many cancers, such as breast cancer and prostate cancer [16,17]. However, the role of FBL in hepatocellular carcinoma cell growth in vivo is unclear. Therefore, we subcutaneously injected Huh7 cells stably transfected with the negative control lentivirus (shNC) and the FBL knockdown lentivirus (shFBL) into nude mice, respectively. As shown in Figure 11A, the knockdown efficiency of FBL was verified by Western blotting analysis, indicating that the silencing efficiency of FBL knockdown lentivirus was significant. The results showed that the knockdown of FBL significantly inhibited tumor growth and tumorigenesis compared with the control group ( Figure 11B-D).   FBL in hepatocellular carcinoma cell growth in vivo is unclear. Therefore, we subcutaneously injected Huh7 cells stably transfected with the negative control lentivirus (shNC) and the FBL knockdown lentivirus (shFBL) into nude mice, respectively. As shown in Figure 11A, the knockdown efficiency of FBL was verified by Western blotting analysis, indicating that the silencing efficiency of FBL knockdown lentivirus was significant. The results showed that the knockdown of FBL significantly inhibited tumor growth and tumorigenesis compared with the control group ( Figure 11B-D).

Discussion
Lung metastasis is an important factor leading to poor prognosis in hepatocellular carcinoma patients. In-depth research on the molecular mechanism associated with lung metastasis in hepatocellular carcinoma will help us to predict the prognosis of hepatocellular carcinoma more accurately and provide a scientific basis for its prevention and treatment.

Discussion
Lung metastasis is an important factor leading to poor prognosis in hepatocellular carcinoma patients. In-depth research on the molecular mechanism associated with lung metastasis in hepatocellular carcinoma will help us to predict the prognosis of hepatocellular carcinoma more accurately and provide a scientific basis for its prevention and treatment.
Many previous studies found that the phenomenon of the over-synthesis of ribosomes exists in a variety of tumors, and it plays an important role in the occurrence and development of tumors [18,19]. As the "factories" of protein synthesis, ribosomes are mainly composed of ribosomal proteins and ribosomal RNA. Ribosomal proteins are divided into large ribosomal proteins (RPL) and small ribosomal proteins (RPS). Previous studies showed that ribosomal proteins mainly participate in the process of protein synthesis. However, with the deepening of research, researchers have found that ribosomal proteins are also involved in the activation of proto-oncogenes or tumor suppressor genes that regulate the cell cycle through multiple signal pathways to promote the growth and proliferation of cells [20]. In addition, a variety of ribosomal proteins have been reported to be associated with tumor metastasis. For instance, the high expression of RPS6 in lung cancer is significantly related to the increased risk of early metastasis [21]. Knocking down RPL39 could significantly inhibit the occurrence of lung metastasis of breast cancer in immunodeficiency mouse models, and most breast cancer patients with lung metastasis have acquired mutations of RPL39 [22]. Knocking out RPL15 using CRISPR technology significantly weakened lung metastasis in breast cancer circulating tumor cells in mouse models, and with RPL15 overexpression, breast cancer circulating tumor cells had stronger lung and ovarian metastasis ability [23]. However, the tumor metastasis mechanism affected by the ribosomal protein has not been clearly reported in the past research. According to the clues provided in the existing research, there are two types of mechanisms: the first is that a variety of ribosomal proteins could promote the composition and activity of ribosomes, which might affect the translation efficiency of specific mRNA types [24]; the other is that it might affect the activation or inhibition of signal pathways through the binding of specific proteins, thereby affecting tumor metastasis [20]. Our study found that ribosomal protein genes account for a large proportion of the modules related to lung metastasis analyzed by WGCNA, such as RPS15, RPS16, RPL18, RPL35A, RPLP2, RPL27A, RPL24, and RPS11. These ribosomal proteins are likely to play important roles in lung metastasis associated with hepatocellular carcinoma, but the related mechanisms need further verification through in vivo and in vitro experiments.
Fibrillarin (FBL) is a 2 -O methyltransferase that is involved in ribosomal biogenesis [25]. It is the main constituent protein of box C/D snoRNP (small nucleolar ribonucleoprotein particles). It is upregulated in prostate cancer, breast cancer, and other tumors, and it can promote tumor cell growth, proliferation, and drug resistance [26][27][28]. Additionally, recent research suggested that increased Fibrillarin expression is related to poor tumor prognosis in hepatocellular carcinoma [29]. However, few studies have investigated whether FBL can promote the progression of hepatocellular carcinoma. Our research preliminarily proves that FBL could affect the proliferation, stemness, migration, and invasion of hepatocellular carcinoma cells, which might be involved in the activation of MYC and E2F signaling. In previous research, it was shown that the tumor suppressor gene TP53 and the proto-oncogene MYC are upstream regulators of FBL [16,27]. In breast cancer cells, TP53 can inhibit the expression of FBL. In prostate cancer, however, the expression of FBL can be upregulated by MYC. As for hepatocellular carcinoma, further study is required to determine which gene plays the main regulatory role in FBL expression. Regarding the downstream regulation of FBL, in addition to participating in the regulation of ribosome synthesis, the latest research shows that FBL could directly modify specific types of mRNA by 2 -O methylation, thereby inhibiting its translation efficiency and achieving the regulation of specific gene expression, which suggests that the effect of FBL on the biological function of hepatocellular carcinoma may also play a role in this regard [30]. Besides, our study demonstrated that genes related to MYC and E2F signaling pathway, such as BYSL, IPO4, PES1, SSB, HNRNPA3, RSL1D1, APEX1, NPM1, EIF3D, RUVBL2, TRIM28, IMPDH2, ILF3, HMGA1, and XRCC6, were strongly positively correlated to FBL, which means that they may be the downstream molecules responsible for FBL-induced lung metastasis.
In general, our research provides ideas for the prediction of the prognosis of hepatocellular carcinoma patients as well as a scientific basis for research on the mechanisms associated with lung metastasis occurrence in hepatocellular carcinoma.

Conclusions
Ribosome-biogenesis-related proteins may play an important role in lung metastasis associated with hepatocellular carcinoma. As a hub gene in ribosome-biogenesis-related proteins, FBL, which is closely related to the poor prognosis of HCC patients, can promote the growth and metastasis of hepatocellular carcinoma.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/bioengineering9080396/s1, Table S1: The hub genes in Turquoise module and Light cyan module and 67 hub genes related to lung metastasis coming from unit of Turquoise module and Light cyan module; Table S2: The expression of 67 genes between normal liver tissues and hepatocellular carcinoma tissues; Table S3: The expression of 67 genes between hepatocellular carcinoma tissues without metastasis and hepatocellular carcinoma tissues with lung metastasis; Table S4: 53 genes correlated to FBL with absolute coefficient > 0.4 in TCGA, GSE14520 and PDC000198.  . are really grateful to Jianyou Liao (Sun Yat-Sen Memorial Hospital, Sun Yat-Sen University) for his suggestions regarding the writing and data analysis. Additionally, we would like to express our gratitude to Jianming Zeng (University of Macau) and all members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

Conflicts of Interest:
The authors declare no conflict of interest.