Lillian I Dolapchiev,1,* Kristyn A Gonzales,1,* Lorenzo R Cruz,1 Mihai Gagea,2 Heather L Stevenson,3 Suet-Ying Kwan,1 Laura Beretta1

1Department of Molecular and Cellular Oncology, The University of Texas MD Anderson Cancer Center, Houston, Texas, USA; 2Department of Veterinary Medicine & Surgery, The University of Texas MD Anderson Cancer Center, Houston, Texas, USA; 3Department of Pathology, The University of Texas Medical Branch, Galveston, Texas, USA

Correspondence: Laura Beretta, Department of Molecular and Cellular Oncology, The University of Texas MD Anderson Cancer Center, 1515 Holcombe Boulevard, Houston, TX, 77030, USA, Tel +1 713 792 9100, Email [email protected]

Purpose: Hepatocellular carcinoma (HCC) related to metabolic dysfunction-associated steatotic liver disease (MASLD) is often diagnosed at a late stage, and its incidence is increasing. Predictive biomarkers are therefore needed to identify individuals at high risk of HCC. We aimed to characterize the gut microbiome and hepatic transcriptome associated with HCC development in female mice with hepatocyte-deletion of Pten (HepPten−). These mice present with large variations in HCC development, making them a powerful model for biomarker discovery.
Methods & Results: Sequencing of stool 16S and hepatic RNA was performed on a first set of mice. Among all liver histology parameters measured, the strongest association with microbiome composition changes was with the number of tumors detected at necropsy, followed by inflammation. The gut microbiome of mice with more than 2 tumors was enriched with Lachnospiraceae UCG and depleted of Palleniella intestinalis and Odoribacter. In contrast, hepatic transcriptomic changes were most strongly associated with tumor burden, followed by liver fibrosis. The 840 differentially expressed genes correlating with tumor burden were enriched in leukocyte extravasation and interleukin 10 receptor A (IL10RA) pathways. In addition, the abundance of Spp1-high epithelial cells is correlated with tumor burden. Association between tumor number and depletion of Palleniella intestinalis, and between tumor burden and circulating levels of C-X-C motif chemokine ligand 13 (CXCL13) and stem cell factor (SCF), was further validated in an independent set of mice.
Conclusion: We identified microbiome components contributing to liver carcinogenesis by inducing inflammation, and changes in hepatic gene expression and hepatic cells distribution that contribute to tumor growth. Such information can be highly valuable for the development of new prevention strategies as well as of new biomarkers for risk modeling in HCC.

Introduction

Approximately one in four adults worldwide is affected by non-alcoholic fatty liver disease (NAFLD), recently renamed metabolic dysfunction-associated steatotic liver disease (MASLD).1 MASLD is a progressive disease, which encompasses a spectrum of chronic liver damage. Patients initially present with simple steatosis that can progress to steatohepatitis (MASH), cirrhosis or hepatocellular carcinoma (HCC). MASH is typically characterized by hepatic fibrosis and inflammation and is a leading cause of cirrhosis in adults in the United States.2,3 Advanced liver fibrosis is the main determinant of mortality in MASH patients and a major risk factor for HCC.1 Although the incidence rate of HCC in MASH patients is lower than that of viral hepatitis patients, the growing obesity epidemic suggests that the number of MAFLD-related HCC patients will continue to increase.4,5

Sensitive prediction models are needed to identify those with MASH at risk of developing HCC.6 Intestinal dysbiosis, increased permeability of the intestinal epithelium, and translocation of microbiota have been shown to induce chronic inflammation, leading to the progression of liver disease.7 It is well recognized that the gut microbiome contributes to the development of MASH,8,9 progression to advanced liver fibrosis,10–15 and even HCC.16–18 The use of gut microbiome profiles to predict HCC outcome is therefore a potential complementary strategy to the use of host hepatic molecular factors.

Sex differences in hepatocarcinogenesis are well known, with men at higher risk for HCC development. Such gender disparity has been well documented in different etiologies of HCC including diabetes with a reduced female-to-male ratio of relative risk19 and alcohol consumption with male mice more susceptible to HCC incidence.20 At the molecular level, several drivers of gender disparity in hepatocarcinogenesis have been reported. These include MyD88-dependent interleukin-6 production,21 epidermal growth factor receptor,22 TGFβ-SMAD3 signaling,20 lipid metabolism and extracellular matrix organization,23 estrogen receptor 1 inhibition of Wnt/β-catenin signaling,24 genetic effects on gene expression,25 Hyperpolarization-activated cyclic nucleotide-gated cation channel 3.26 Gut microbiota was also reported to contribute to gender disparity in hepatocarcinogenesis.27,28 Similar gender disparity is seen in all mouse models of HCC, and therefore, female mice are underrepresented in mouse studies of hepatocarcinogenesis. Gender difference in HCC incidence in HepPten− was previously reported, with most male mice developing liver tumors, while only 50% of the female mice developed tumors.29 Female mice therefore represent a powerful model to dissect the molecular changes associated with HCC development in the context of MASH.

In this study, we aimed to identify alterations in the liver transcriptome and gut microbiome associated with HCC development and growth. To that end, female mice with hepatocyte-specific deletion of Pten (HepPten−) mice were used. These mice develop MASH and HCC.30 In this model, and in contrast to male mice, the female mice exhibit great heterogeneity with respect to the age of tumor onset.

Materials and Methods Mouse Necropsies and Histology Analysis

All animal work was performed in accordance with the approval of the Institutional Animal Care and Use Committee of The University of Texas MD Anderson Cancer Center (IACUC protocol #00000897-RN03, approved 04/22/2021) and conducted in accordance with the Animal Welfare Act, the Guide for the Care and Use of Laboratory Animals and the Public Health Service Policy. Experimental mice were female PtenloxP/loxP;Alb-Cre+ (HepPten−) and were obtained through the crossing of PtenloxP/loxP;Alb-Cre+male mice and PtenloxP/+;Alb-Cre+ female mice (C57BL/6 background). Stool samples were collected the day prior to necropsy, in a sterile tube that was snap frozen and stored at −80°C. At the time of necropsy, body, liver and cecum weight were collected together with blood by cardiac puncture, liver tissue, liver tumors, and colon. Tumor diameters and colon length were measured. Tumor and liver sections were also snap-frozen in liquid nitrogen. Formalin-fixed tissue sections were sectioned and stained with hematoxylin and eosin (H&E) and/or Masson’s Trichrome. Histopathologic analysis was performed blindly by pathologists. We adopted the NAFLD Activity Score (NAS)31 to include both macrovesicular and microvesicular steatosis grades and calculated as an unweighted sum of the grading scores of macrosteatosis (0–3), microsteatosis (0–3), lobular inflammation (0–3) and ballooning (0–2). Bile ductular reaction and liver fibrosis were also assessed on a 0 to 3 and 0 to 4 scale, respectively. Aspartate aminotransferase (AST) and alanine aminotransferase (ALT) measurements in serum (30μL of 6-fold dilution) were determined using the ACE Axcel clinical chemistry system (Diagnostic Technologies).

Stool DNA Extraction, 16S rRNA Amplicon Sequencing and Bioinformatic Analysis

Stool samples were analyzed for 16S sequencing at the MD Anderson Cancer Center Microbiome Core Facility. Bacterial genomic DNA was extracted using the QIAamp Fast DNA stool mini Kit (Qiagen). The V4 region of the bacterial 16S rDNA gene was amplified by PCR (forward primer: 5’-AATGATACGGCGACCACCGAGATCTACACGCTXXXXXXXXXXXXTATGGTAATTGTGTGYCAGCMGCCGCGGTAA-3’, where XXXXXXXXXXXX is an index sequence for multiplexing libraries; reverse primer: 5’-CAAGCAGAAGACGGCATACGAGATAGTCAGCCAGCCGGACTACNVGGGTWTCTAAT-3’). Libraries were purified using Zymo I-96 column purification and analyzed on the Agilent 4200 Tapestation system (Agilent). Barcoded amplicons were pooled at equal concentrations. Pooled libraries were quantified by Qubit fluorometer, and the molarity was calculated based on amplicon size. Sequencing was performed using 250 bp paired-end on the Illumina MiSeq platform (Read1 seq primer: 5’-TATGGTAATTGTGTGYCAGCMGCCGCGGTAA-3’; Read2 seq primer: 5’-AGTCAGCCAGCCGGACTACNVGGGTWTCTAAT-3’; and index primer: AATGATACGGCGACCACCGAGATCTACACGCT. Paired-end reads were de-multiplexed and split in QIIME. Merging of paired-end reads to create consensus sequences was done by VSEARCH v7, allowing up to a maximum of 10 mismatched. The cluster_otus command, an implementation of UPARSE algorithm, was used to perform 97% related operational taxonomic units (OTU) clustering. Denoising was done by the unoise3 command. OTUs were then subject to taxonomy assignment using the Mothur with Silva database.

RNA Sequencing and Cell Deconvolution

RNA sequencing was performed at the University of Texas MD Anderson Advanced Technology Genomics Core. Following cDNA synthesis, end repair, adaptor ligation and PCR amplification with the KAPA Stranded mRNA-Seq kit, the enriched cDNA libraries were subjected to strand-specific mRNA sequencing on the Illumina NovaSeq 6000 system, using the SP flow cell and 100 base-pair paired-end reads. Sequence files were generated in FASTQ format (bcl2fastq, v2.20.0). Using the STAR program (v2.7.0f), raw reads were aligned to the mouse reference genome GRCm39, and the number of aligned read per gene using the annotations GRCm39.103 was counted. Genes with maximum abundance below 0.5 counts per million (CPM) across all samples were excluded, leaving 17,315 genes for differential expression analysis. Criteria for differentially expressed genes (DEGs) were DESeq2 adjusted p values ≤0.05 and FC ≥1.5 or ≤-1.5. Furthermore, the abundance of DEGs had to significantly correlate with tumor burden (Spearman correlation Benjamini–Hochberg (BH) adjusted p values ≤0.05). To impute relative proportions of cell populations in each liver sample, CIBERSORTx32 was implemented using total liver RNA sequencing data. Sequencing data normalized to CPM was used for analysis. The reference signature matrix for adult mouse liver was obtained from a library of cell profile matrices at: https://github.com/Nanostring-Biostats/CellProfileLibrary.33

Enzyme-Linked Immunoassays (ELISA)

CXCL13 concentration was measured in serum samples using the Quantikine ELISA Mouse CXCL13/BLC/BCA-1 Immunoassay (catalog #MCX130, R&D Systems (bio-techne), Minneapolis) according to manufacturer’s instructions. Kitl (SCF) concentration was measured in EDTA plasma samples using the Quantikine ELISA Mouse SCF Immunoassay kit (catalog #MCK00, R&D Systems (bio-techne), Minneapolis) following manufacturer’s instructions. Optical density was read using a Molecular Devices SpectraMax ABS Plus Microplate Reader set at 450 nm with the correction wavelength set at 570 nm.

Stool DNA Quantitative PCR (qPCR)

DNA was extracted from stool samples using the NucleoSpin DNA Stool kit (#740472.50, Macherey-Nagel). DNA concentration was measured using a NanoDrop 1000 Spectrophotometer (Thermo Scientific). Extracted DNA samples were diluted to 1 ng/µL aliquots. The species-specific forward primer (5’-GAACTGGCGTGCTTGAGTCC-3’) and reverse primer (with 5’-TTCGAGCATCAGCGTCAGTT-3’) were designed by NCBI primer-BLAST based on the most abundant Prevotellaceae OTU observed by 16S sequencing. That OTU sequence perfectly aligned to Palleniella intestinalis by NCBI nucleotide BLAST. The absolute concentration of Palleniella intestinalis was determined by qPCR, in each DNA stool sample using 6ng of DNA and iTaq Universal SYBR Green SuperMix (#172-5124, Bio-Rad). Each qPCR plate was run in a CFX Connect Real-Time System (Bio-Rad) with the following cycling conditions: denaturation at 95°C for 5 minutes for 1 cycle, and for 40 cycles, denaturation at 95°C for 15 seconds, annealing at 60°C for 30 seconds and extension at 72°C for 40 seconds. The number of species-specific 16S gene copies was calculated from a standard curve, generated using Palleniella intestinalis genomic DNA from the Leibniz-Institute DSMZ-German Collection of Microorganisms and Cell Cultures.

Statistical Analyses

Unless specified otherwise, the statistical analyses were performed in R (version 4.1.2; R Foundation for Statistical Computing, Vienna, Austria). Principal coordinates analysis (PCoA) was performed using the “cmdscale” function and the weighted UniFrac distances of the OTU tables. Beta dispersion and permutational multivariate analysis of variance (PERMANOVA) tests, using weighted UniFrac distances, were performed with the Vegan package. Partitioning Around Medoids (PAM) clustering based on weighted UniFrac distances for microbiome data and Euclidean distances for gene expression data was performed using the “pam” function in the cluster package. The “fviz_nbclust” function in the factoextra package suggested four clusters. Kruskal–Wallis test was performed to identify bacterial taxa with significantly different abundance across clusters. To evaluate the relationship between microbiome composition and tumor and histology parameters, distance-based redundancy analysis (RDA) was performed with the “capscale” function in the Vegan package, using weighted UniFrac distances as the response variable, and tumor burden, tumor number, inflammation, and BDR severity as explanatory variables. To evaluate the relationship between hepatic gene expression profiles and liver tumor development and fibrosis severity, RDA was again performed with the “capscale” function in the Vegan package, using gene expression as the response variable and tumor number, tumor burden, and fibrosis stage as explanatory variables. Analysis of variance (ANOVA)-like permutation and marginal tests were performed using the “anova.cca” function in the Vegan package, to determine the significance of the model and of each explanatory variable. Fitting variables onto RDA plots were performed using the “envfit” function in the Vegan package, using linear constraints. Differences in bacterial abundance were assessed using the linear discriminant analysis (LDA) effect size (LefSe) tool,34 with p<0.05 and log10 LDA score >2 considered significant. To determine associations between bacterial taxa and histology parameters, mice were divided into low and high groups based on median values. Pairwise correlations were performed using Spearman correlation. Correlation p-values were adjusted for multiple testing by the Benjamini–Hochberg method. A Shapiro–Wilk test was used to test for the normal distribution of data prior to statistical analyses. Mann–Whitney test was used to compare significant (p < 0.05) associations between bacterial taxa in mice with high or low tumor number. To determine whether a relationship exists between tumor burden and serum CXCL13 or plasma SCF concentrations, Spearman correlation was used. Spearman correlation was also used to determine the relationship between CXCL13 or SCF gene expression and CXCL13 serum or SCF plasma protein concentrations, respectively. Significance was defined as a p-value <0.05. To determine whether a relationship exists between Palleniella intestinalis abundance in stool and tumor number, Spearman correlation was used.

Results Liver Tumor Development in HepPten− Female Mice

The presence of liver tumors was analyzed in a first set of 25 HepPten− female mice at the time of necropsy, when mice were around 13 months old. Six among the 25 mice (24%) did not have any tumor. For the remaining 19 mice, the number of tumors per mouse ranged from 1 to 16, with an average of 4.5 tumors and a median of 3 tumors (Figure 1A). For these mice with tumors, the tumor burden ranged from 8.4 to 4774.1 mm3 (mean = 963.3 mm3; median = 388.4 mm3) (Figure 1A). Liver histological parameters pertaining to MASH and fibrosis were also analyzed in these 25 HepPten− female mice. Liver fibrosis scores and bile ductular reaction (BDR) scores ranged from 0 to 3 and 0 to 1.5, respectively, with medians for both of 0.5. Ten mice (40%) presented with no liver fibrosis and 7 mice (28%) presented with no BDR (Figure 1B). The extent of NAFLD severity also showed large variation with modified NAFLD scores (mNAS) ranging from 0.5 to 6.5 (median = 3.5) (Figure 1C). Fifteen mice (60%) had no macrovesicular steatosis, while 1 mouse (4%) had no detectable microvesicular steatosis. Inflammation scores ranged from 0 to 1.5 (median = 0.50), while ballooning scores ranged from 0 to 2 with a low median of <0.5 (Figure 1C).

Figure 1 Liver tumor and liver histology characteristics for the 25 HepPten− female mice. (A) Scatterplot showing number of tumors and tumor burden per mouse observed in 25 HepPten− female mice. Liver histological parameter scores for (B) liver fibrosis, bile ductular reaction and for (C) MASH parameters in these 25 mice. Median and range bars are shown.

Gut Microbiome Changes Associated with Liver Tumor Development

Stool samples collected from 23 of the 25 mice, the day prior to necropsy were analyzed using 16S sequencing. Relationship between gut microbiome composition and liver tumor development was revealed by PAM clustering. PCoA plots of weighted UniFrac distances resulted in three clusters: cluster A (n = 4), cluster B (n=11), and cluster C (n = 8) (Figure 2A). Among all collected necropsy and liver histological parameters, tumor burden, number of tumors, inflammation and BDR showed significant differences between clusters (Supplementary Table S1). Mice in cluster C had the highest number of tumors (mean = 6.9; p = 0.005) as well as the highest average tumor burden (1418 mm3; p = 0.040). Mice in cluster C also had high liver inflammation scores (mean = 0.6; p = 0.003) and BDR scores (mean = 0.6, p = 0.003). In contrast, mice in cluster B had the lowest average number of tumors (1.5; p = 0.015), inflammation scores (0.2; p=0.017), and BDR scores (0.2; p = 0.017). To further evaluate the relationship between microbiome composition and liver tumor development and growth, and inflammation and BDR severity, redundancy analysis (RDA) was performed based on weighted UniFrac distances and using these tumor and histology parameters as explanatory variables (Figure 2B). The overall model was highly significant with the following four parameters (number of tumors, tumor burden, inflammation scores and BDR scores) explaining 30.85% of the variation in microbiome composition (p=0.016). A significant association was, however, only observed with the number of tumors (9.94% variance explained, permutation test p = 0.034).

Figure 2 Relationship between stool microbiome profiles and liver parameters. (A) PCoA plot of stool microbiome profiles from 23 female HepPten- mice, based on weighted UniFrac distances. Samples are grouped by PAM clusters. (B) Triplot of the weighted UniFrac distance-based RDA, which was performed to assess the relationship between selected liver parameters (explanatory variables) and the stool microbiome profiles. p-values and % variance explained with the ANOVA-like significance test of the model and the marginal test of the explanatory variables, are shown.

To identify the specific taxa associated with number of tumors, LEfSe analysis was performed from kingdom to species. A total of 28 taxa were identified to be significantly associated with the number of tumors >2 (Figure 3A). Of these 28 taxa, 13 taxa were significantly depleted in mice with high tumor number, while the remaining 15 taxa were significantly enriched in mice with high tumor number (Figure 3B). The abundance of 27 of these bacteria showed significant positive or negative correlation with the number of tumors detected in these mice, as measured by Spearman correlation coefficient (Figure 3C). Strong enrichment in mice with more than 2 tumors and positive association with number of tumors were observed for Lachnospiraceae UCG-006 (fold change [FC]=5.6, p = 0.021; rs=0.44, p = 0.035), Monoglobus (FC = 30.8, p = 0.001; rs = 0.67, p < 0.001) and Gordonibacter (FC = 4.9, p = 0.004; rs = 0.56, p = 0.005). Strong depletion and negative association were observed for Odoribacter (FC = −8.3, p = 0.002, rs = −0.65, p = 0.001) and Prevotellaceae UCG-001 (FC = −3.8, p < 0.001, rs=−0.68, p < 0.001).

Figure 3 Bacterial taxa with altered abundance in mice with high tumor numbers. (A) Cladogram showing bacterial taxa with significantly different abundances between mice with high (>2) and low tumor (≤2) numbers, as assessed by the LEfSe algorithm. Classifications at the phylum (p_), class (c_), order (o_), family (f_), genus (g_), and species (s_) levels are shown. (B) Corresponding volcano plots showing all taxa that are significantly depleted or enriched in mice with a high tumor number, with taxa labelled as in the cladogram. Significance (p<0.05) was assessed using a Mann–Whitney test. The x-axis represents the base 2 log of the fold change (high/low); the y-axis represents the negative base 10 log of the p-value for each taxon. (C) Spearman correlation matrix between tumor number and significant taxa from the LEfSe analyses. + indicates a significant correlation (p<0.05).

Association Between Prevotellaceae UCG-001/Palleniella Intestinalis and Tumor Number: Validation in an Independent Set of Mice

At the genus level, 94 OTUs were assigned to Prevotellaceae UCG-001 in the 16S sequencing data. The most abundant OTU accounted for 83.6% of all reads and aligned perfectly to Palleniella intestinalis, a member of the Prevotellaceae family. The second most abundant OTU accounted for only 4.8% of all reads assigned to Prevotellaceae UCG-001. To validate the negative correlation between stool content of Prevotellaceae UCG-001 and tumor number, we designed a set of primers based on the most abundant OTU and therefore quantified Palleniella intestinalis in the stool of the same set of 25 mice, by qPCR. For this set of mice, Palleniella intestinalis concentration ranged from 633 to 132,323,834 16S gene copies per mg stool (median = 24,493). The mouse with the lowest Palleniella intestinalis abundance had 3 tumors, while the mouse with the highest abundance had 1 tumor. A strong negative correlation was observed between Palleniella intestinalis abundance and tumor number (rs = −0.72, p < 0.001) (Figure 4A), confirming the inverse relationship between level of this member of the Prevotellaceae family and tumor number.

Figure 4 Relationship between Palleniella intestinalis concentration and tumor number in mice in discovery and validation groups. Scatterplots showing Spearman correlation between the concentrations of Palleniella intestinalis and tumor number in (A) the HepPten− discovery female mice group, and (B) the HepPten− validation female mice group. Spearman correlation p-values and rank correlation coefficients are shown.

The association between stool Palleniella intestinalis abundance and tumor numbers was further validated in an independent set of 71 HepPten− female mice. At necropsy, ten of the 71 mice (14%) did not have any tumors. For the remaining 61 mice, the number of tumors per mouse ranged from 1 to 18, with an average of 3.6 tumors and a median of 3 tumors (Supplementary Figure S1A). For the mice with tumors, the tumor burden ranged from 11.2 to 5544.8 mm3 (mean = 865.9 mm3; median = 252.2 mm3) (Supplementary Figure S1A). Liver histological parameters pertaining to MASH and fibrosis were also analyzed in these 71 HepPten− female mice. Liver fibrosis scores and BDR scores ranged from 0 to 2 and 0 to 3, respectively, and with medians of 0 and 0.5, respectively. Thirty-nine mice (55%) presented with no liver fibrosis and 29 mice (41%) presented with no BDR (Supplementary Figure S1B). The extent of NAFLD severity also showed large variation with mNAS ranging from 0 to 7.5 (median = 4.5). Macrovesicular scores ranged from 0 to 3 (median = 0), and microvesicular scores ranged from 0 to 3 (median = 3). Thirty-six mice (51%) had no macrovesicular steatosis, while 1 mouse (1.4%) had no detectable microvesicular steatosis. Inflammation scores ranged from 0 to 3 (median = 0.50), while ballooning scores ranged from 0 to 2 with a median of 0.5 (Supplementary Figure S1C). Palleniella intestinalis concentrations ranged from 6,270,758 to 122,607,493 16S gene copies per mg of stool (median=34,261,633). The mouse with the lowest Palleniella intestinalis abundance had 8 tumors, while the mouse with the highest abundance had only 1 tumor. A significant negative correlation was observed between Palleniella intestinalis abundance and tumor number (rs = −0.42, p < 0.001) (Figure 4B), confirming an inverse relationship between Palleniella intestinalis and tumor number.

Liver Transcriptome Changes Associated with Liver Tumor Growth

RNAseq was performed on liver collected at necropsy from the 25 mice. Relationship between hepatic gene expression and liver tumor development was revealed by PAM clustering. PCoA plots resulted in three clusters: cluster A (n = 8), cluster B (n=7) and cluster C (n = 10) (Figure 5A). Among all collected necropsy and liver histological parameters, liver weight, tumor burden, number of tumors, fibrosis and ALT levels, showed significant differences between clusters (Supplementary Table S2). Mice in cluster B had the highest average liver weight (5.19 g; p < 0.001), the highest average number of tumors (7.0; p=0.008), as well as the highest average tumor burden (2264 mm3; p < 0.001). Mice in cluster B also had the highest average plasma levels of ALT (279.4 U/L, p = 0.034). Mice in cluster A had the lowest average fibrosis scores (0.3; p = 0.046). To further evaluate the relationship between hepatic gene expression in these 25 mice and liver tumor development and fibrosis severity, RDA was performed using these tumor and histology parameters as explanatory variables (Figure 5B). The overall model was highly significant with the following three parameters (number of tumors, tumor burden, fibrosis scores) explaining 48.3% of the variation in gene expression (p = 0.001). A significant association was observed with tumor burden (17.51% variance explained, permutation test p = 0.002) and fibrosis (7.41% variance explained, permutation test p = 0.047).

Figure 5 Relationship between liver transcriptome profiles and tumor growth. (A) PCoA plot of liver RNAseq gene expression from the same 25 female HepPten− mice. Samples are grouped by PAM clusters. (B) Triplot of the gene expression RDA, which was performed to assess the relationship between selected liver parameters (explanatory variables) and the liver RNAseq data. p-values and % variance explained with the ANOVA-like significance test of the model and the marginal test of the explanatory variables, are shown. (C) Volcano plot showing the 840 DEGs identified in the liver of mice by DESeq2, all significantly correlated with tumor burden. 706 DEGs were enriched (colored in yellow) in mice with high tumor burden while the remaining 134 DEGs were depleted (colored in blue) in mice with high tumor burden. The base 2 log of the fold change (high/low) was graphed on the x-axis while the negative base 10 log of the DESeq2 adjusted p-values were graphed on the y-axis. Genes that were labeled represent the genes that were most depleted or most enriched. (D) Bar graph (left) showing the significant p-values for 5 canonical pathways enriched by DEGs with tumor burden as determined by IPA analysis and a bar graph (right) showing the significant p-values for 5 upstream regulators.

A total of 840 DEGs, 706 upregulated and 134 downregulated, were identified in liver of mice with high tumor burden (>89 mm3) (Figure 5C). All DEGs were selected for having significant correlation with tumor burden (Supplementary Table S3). Genes with strong abundance changes and positive correlation with tumor burden included solute carrier family 4 member 9 (Slc4a9) (FC = 13.1, p < 0.001; rs = 0.73, p = 0.003), RP1 axonemal microtubule associated (Rp1) (FC = 12.2, p = 0.010; rs = 0.66, p = 0.009), protein phosphatase 2 regulatory subunit B-gamma (Ppp2r2c) (FC = 9.3, p < 0.001; rs = 0.74, p = 0.003), tumor-associated calcium signal transducer 2 (Tacstd2) (FC = 6.5, p = 0.001; rs = 0.64, p = 0.011), oncoprotein-induced transcript 1 (Oit1) (FC = 4.2, p = 0.001; rs=0.70, p = 0.005), galactose-3-O-sulfotransferase 2 (Gal3st2) (FC = 3.7, p = 0.001; rs = 0.80, p = 0.002), chemokine (C-C motif) ligand 8 (Ccl8) (FC = 3.6, p = 0.004; rs= 0.71, p = 0.004), CXCL13 (FC = 3.4, p < 0.001; rs=0.71, p = 0.005), follistatin (Fst) (FC = 3.0, p < 0.001; rs = 0.75, p = 0.003), and thrombospondin 4 (Thbs4) (FC = 2.7, p = 0.034; rs=0.73, p = 0.003). The strongest depletion and negative associations were observed for several genes encoding for major urinary proteins: Mup9 (FC=−5.7, p = 0.001; rs=−0.73, p = 0.003), Mup17 (FC = −3.8, p = 0.016; rs=−0.75, p = 0.003), Mup1 (FC = - 3.7, p = 0.009; rs= −0.69, p = 0.006), Mup15 (FC = −3.4, p = 0.043; rs= −0.74, p = 0.003), Mup11 (FC = −3.3, p = 0.005; rs= −0.73, p = 0.003), Mup19 (FC = −3.0, p = 0.007; rs=−0.72, p = 0.004), Mup13 (FC = −2.8, p = 0.026; rs=−0.72, p = 0.004), Mup22 (FC = −2.6, p = 0.018; rs= −0.71, p = 0.005) and Mup2 (FC = −2.1, p = 0.001; rs= −0.74, p = 0.003). Other genes included aquaporin 4 (Aqp4) (FC = −2.8, p < 0.001; rs= −0.75, p = 0.003), glycerol-3-phosphate acyltransferase, mitochondrial (Gpam) (FC = −2.1, p < 0.001; rs= −0.72, p = 0.004), thyrotropin releasing hormone degrading enzyme (Trhde) (FC = −2.0, p < 0.001; rs= −0.76, p = 0.002), the retinoic acid early transcripts 1E (Raet1e) (FC=−1.8, p = 0.001; rs=−0.78, p=0.002), and 1Delta (Raet1d) (FC = −1.8, p = 0.001; rs= −0.74, p = 0.003). Overall, DEGs associated with tumor burden were strongly enriched in leukocyte extravasation signaling (p = 1.6x10−7) and Th1 pathway (p = 2.3x10−7), as determined by IPA analysis. Upstream regulators included interleukin 10 receptor A (IL10RA) (p = 3.8x10−12) and ubiquitin-specific peptidase 22 (USP22) (p = 3.2x10−11) (Figure 5D).

To determine whether increased tumor burden was associated with changes in cellular distribution in the liver, cell deconvolution analysis using a mouse liver matrix was performed on the hepatic transcriptomic profiles obtained from the 25 mice. As anticipated, hepatocytes were the dominant cell type with a significant reduction of Fabp1-high hepatocytes in mice with high tumor burden (28.0% vs 52.6%, FC=−1.9, p < 0.001) and concomitant increase in progenitor cells labeled as Spp1-high epithelial cells (3.3% vs 1.0%, FC = 3.4, p = 0.010). Endothelial cells, Kupffer cells, and Chil3-high macrophages were also significantly enriched in liver with high tumor burden (1.7% vs 0.3%, FC = 6.2, p = 0.007; 8.5% vs 3.5%, FC = 2.5, p < 0.001; 1.4% vs 0.6%, FC = 2.2, p = 0.028, respectively) (Figure 6A). Among these cell types, strong positive correlations were observed between tumor burden and presence of Kupffer cells (rs = 0.65, p < 0.001) or Spp1-high epithelial cells (rs=0.60, p = 0.001). In contrast, a strong negative correlation was observed between tumor burden and depletion of Fabp1 hepatocytes (rs = −0.73, p < 0.001) (Figure 6B).

Figure 6 Relationship between hepatic cell types and tumor burden in discovery mice. (A) Cellular distribution of cell-types (relative percent) in the liver of the 25 discovery mice with high tumor burden (11 left-most vertical bars) versus low tumor burden (14 right-most vertical bars). Cell deconvolution analysis with a mouse liver matrix was used. (B) Scatterplots showing Spearman correlation between tumor burden and the presence of Kupffer cells (left), between tumor burden and the presence of Spp1-high epithelial cells (center), and between tumor burden and the presence of Fabp1 hepatocytes (right). Spearman correlation p-values and rank correlation coefficients are shown.

Circulating CXCL13 and SCF as Biomarkers for HCC Prediction

Hepatic expression of the chemokine CXCL13 and the stem cell factor SCF genes was found to strongly correlate with tumor burden (Supplementary Table S3). Because both CXCL13 and SCF are detectable in blood, we measured circulating concentrations of CXCL13 and SCF in the same set of 25 mice. Serum concentration of CXCL13 ranged from 0.58 to 4.32 ng/mL (median = 1.45 ng/mL). The mouse with the highest CXCL13 level had a large tumor burden of 398.57 mm3. A significant positive correlation was observed between tumor burden and CXCL13 serum concentrations (r = 0.47, p = 0.025) (Figure 7A). We also observed a positive correlation between CXCL13 hepatic gene expression and CXCL13 protein serum concentration (r = 0.64, p < 0.001), suggesting a direct read of CXCL13 hepatic gene abundance in blood (Figure 7B). Similar trends were observed for SCF plasma concentration. SCF concentrations ranged from 101 to 217 pg/mL (median of 135 pg/mL). The mouse with the highest SCF level had a high tumor burden of 769.62 mm3. A strong positive correlation was observed between tumor burden and SCF plasma concentrations (r = 0.60, p = 0.003) (Figure 7C), although SCF hepatic gene expression levels and SCF plasma protein concentrations had only a weak, non-significant correlation (r = 0.29, p = 0.176) (Figure 7D).

Figure 7 Significant association between tumor burden and protein biomarkers CXCL13 and SCF. (A) Scatterplot showing correlation between tumor burden (mm3) and serum CXCL13 protein concentration (ng/mL) in the HepPten−discovery female mice set. (B) Correlation between CXCL13 hepatic expression using RNA sequencing measured in adjusted counts per million (CPM) and serum CXCL13 protein concentration (ng/mL) from the HepPten− discovery female mice set. (C) Scatterplot showing correlation between tumor burden (mm3) and SCF plasma protein concentration (pg/mL) in the HepPten− discovery female mice set. (D) Correlation between SCF hepatic expression using RNA sequencing measured in adjusted counts per million (CPM) and SCF plasma protein concentration (pg/mL) from the HepPten− discovery female mice set. (E) Correlation between tumor burden and serum CXCL13 protein concentration (ng/mL) in the HepPten− validation female mice set. (F) Correlation between tumor burden and SCF plasma protein concentration (pg/mL) in the HepPten− validation female mice set. Spearman correlation was used for all six plots. Spearman correlation p-values and rank correlation coefficients are shown. p<0.05 was considered significant.

Validation of these two protein biomarkers was performed on the independent set of 71 HepPten− female mice described above and in Supplementary Figure S1. In this set of mice, CXCL13 serum concentrations ranged from 0.49 to 7.13 ng/mL (median = 1.16 ng/mL). The mouse with the highest CXCL13 concentration had a very high tumor burden of 1604.1 mm3. Validating the data obtained in the first set of mice, a significant positive correlation was observed between tumor burden and CXCL13 serum concentrations (r = 0.39, p < 0.001) (Figure 7E). SCF plasma concentrations ranged from 54 to 206 pg/mL (median = 81 pg/mL). The mouse with the highest SCF concentration had a tumor burden of 551.7 mm3. In addition to validating the data obtained in the first set of mice, a significant positive correlation was observed between tumor burden and SCF plasma concentrations (r = 0.24, p = 0.049) (Figure 7F).

Discussion

The increasing number of patients with MASLD-associated HCC, combined with the late diagnosis of this disease, warrants the need for preventative measures and earlier identification of those at risk for MASLD-associated HCC. In this study, we aimed to identify molecular markers that can be used to predict HCC risk in mice with MASLD, by correlating tumor development and growth with changes in gut microbiome composition and hepatic gene expression, in HepPten− female mice. Female HepPten− mice exhibit great heterogeneity with respect to the age of tumor onset, representing therefore a relevant model to dissect the molecular events associated with development and growth of MASLD-associated HCC. The most remarkable finding of the study was that gut microbiome changes were strongly associated with number of tumors, indicative of tumor development, while hepatic transcriptomic changes were strongly associated with tumor burden, indicative of tumor growth. Interestingly, gut microbiome changes were also associated with hepatic inflammation, representing therefore a likely intermediate between the gut and liver tumor development. In contrast, hepatic changes correlating with tumor burden were also associated with liver fibrosis severity.

The gut microbiome has been proposed as a tool for non-invasive biomarkers for early HCC.35,36 Overall, gut microbial diversity increased with progression from cirrhosis to HCC, with a decrease of butyrate-producing genera and a concomitant increase of lipopolysaccharide-producing genera.36 Furthermore, probiotics can be effective in reversing steatohepatitis that occurs due to dysregulated bile acid synthesis and dysbiosis in mice.37 In our study, high tumor number was associated with depletion of Prevotellaceae_UCG-001, representing a cluster of uncultured genera. Inulin, a dietary fiber used as a prebiotic to alleviate glucose and lipid metabolism disorders in mice and humans, induced an enrichment of Prevotellaceae UCG-001 that negatively correlated with markers of glycolipid metabolism disorders.38 Alignment of the representative sequences of Prevotellaceae_UCG-001 in our 16S sequencing data indicated that the cluster was primarily comprised of Palleniella intestinalis. Palleniella intestinalis is a novel species first described in the mouse gut as “Prevotella intestinalis”,39 and later reclassified from Prevotella into the novel Palleniella genus in the Prevotellaceae family.40 Genomic interrogation found that Palleniella intestinalis consistently formed a distinct phylogenetic branch outside of the Prevotella genus, as an intermediate between Prevotella, Xylanibacter and Alloprevotella species. Analysis of its genome predicted the production of acetate, propionate and L-glutamate. Analytical profile index testing further indicated that the species possessed beta-galactosidase, alpha-glucosidase, alpha-arabinosidase, leucyl glycin arylamidase, alanine arylamidase and glutamyl-glutamate arylamidase activity.40 Strong depletion in mice with more than two tumors and negative association with tumor number were also observed for Odoribacter. Depletion of Odoribacter in MAFLD patients and of Odoribacter splanchnicus in HCC patients was previously reported.35,41

Finally, strong enrichment in mice with more than 2 tumors and positive association with the number of tumors were observed for Lachnospiraceae UCG-006. Lachnospiraceae colonizes the intestinal lumen of the host. Their species richness and relative abundances increase during the host’s life. Members of Lachnospiraceae are among the main producers of short-chain fatty acids, and some taxa of Lachnospiraceae are associated with extraintestinal diseases.42 Most importantly, Lachnospiraceae has been associated with human liver cancer development. Enrichment of Lachnospiraceae was found in the gut microbiome of HCC patients,43 as well as in human liver tumor tissues.44 In addition, diet can influence the gut microbiome as reviewed in the context of high fat diet and gastrointestinal cancers.45 One particular study found that feeding mice a high-fat diet resulted in an increase in Lachnospiraceae.46 Another study found that a red meat diet in rats promoted an enrichment of Lachnospiraceae_UCG-010, together with increased liver inflammation and fibrosis.47 Future studies on the effect of diet on the microbiome signatures, which we found to be associated with HCC development, would be highly valuable.

Gene expression changes associated with tumor burden in HepPten− female mice were strongly enriched in leukocyte extravasation signaling and Th1 pathway, with IL10RA and USP22 predicted as major upstream regulators. Leukocyte infiltration is a hallmark of hepatic inflammation. L-selectins on leukocytes initiate the leukocyte extravasation process and adhesion to liver sinusoidal endothelial cells.48 Leukocyte extravasation plays an essential role in many pathophysiological processes, including fibrosis and cancer progression.49 Recruitment of activated leukocytes from peripheral blood into the tumor tissue is a crucial step of the immune response. In the context of head and neck cancer, the level of leukocyte extravasation-related genes was higher in cancer-associated fibroblasts than in normal fibroblasts.50 The IL-10 network is among the most important pathways linking cancer and inflammation. Disruption of the interactions of IL-10 with its receptors (IL-10RA and IL-10RB) may enhance inflammation and modulate anti-tumor immunity.51 Interestingly, innate sensing of translocated microbiota induces IL-10R signaling, driving the loss of systemic T-cell immunity during chronic liver injury.52 USP22 is abnormally expressed in various malignant tumors.53 In HCC, high USP22 expression is associated with poor prognosis.54,55 USP22 may stabilize programmed death ligand 1 by deubiquitination while also regulating T-cell infiltration into tumors. Ablation of USP22 in immunosuppressive regulatory T cells leads to reduced tumor burden in several cancers.56 Importantly, USP22 promotes tumorigenesis and progression in HCC, by promoting stemness.53,57

As hepatic expression of CXCL13 and SCF genes strongly correlated with tumor burden, we tested their potential as circulating HCC biomarkers. Indeed, a positive correlation was observed between tumor burden and CXCL13 serum concentrations or SCF plasma concentrations. CXCL13, originally identified as a B-cell chemokine, plays an important role in the immune system. The interaction between CXCL13 and its receptor CXCR5 builds a signaling network that regulates the development of inflammatory diseases and cancers.58 In the liver, CXCL13 mediates the recruitment of intrahepatic CXCR5+CD8+ T cells and promotes aggregation of CD19+ B cells and CXCR5+CD4+ T cells.59 CXCL13 was suggested to be involved in lymphocyte hyperresponsiveness associated with fatty liver inflammation in obese mice.60 SCF is a multifunctional cytokine playing major roles in liver embryogenesis and in the pathogenesis of liver diseases, including HCC and MAFLD.61 Increased expression of SCF in HCC patients is associated with poor prognosis.62 In the developing liver, hepatic stellate and endothelial cells create a niche for perisinusoidal vascular hematopoietic stem cells by producing SCF.63 In the adult liver, cells expressing SCF include progenitor cells, bile epithelial cells, and some hepatocytes, participating in liver tissue repair.64 Serum levels of SCF were reported to be associated with HCC development in patients with hepatitis C.65

Overall, the transcriptomic changes associated with tumor burden suggested a major role of progenitor cells in HCC development and growth, as confirmed by cell deconvolution analysis. Progenitor cells, labeled as Spp1-high epithelial cells in the matrix, were significantly enriched in the livers of mice with high tumor burden and their numbers strongly correlated with tumor burden. A major role of hepatic progenitor cells in HCC development and growth has been extensively reported (reviewed in66). We have previously demonstrated the role of Spp1 in predicting HCC in human cohorts.67,68

The study has some limitations including the use of only female mice. Follow-up studies to confirm the universality of the identified biomarkers, with a second mouse model of MASH-related HCC and with human cohorts would also be valuable. In addition, the research focuses on correlation rather than causation. Follow-up pre-clinical studies testing specific intervention approaches would be necessary to validate potential prevention strategies.

Altogether, these results suggest that specific microbiome components contribute to liver carcinogenesis by inducing inflammation, while changes in hepatic gene expression contribute to tumor progression. The molecular events identified in this study suggest a step-wise influence of the gut microbiome on the liver and HCC tumorigenesis, and a role played by alterations in liver transcriptome and cellular diversity in HCC growth. Such information can be highly valuable for the development of new prevention strategies such as probiotics as well as of risk modeling using non-invasive biomarkers for HCC.

Abbreviations

ALT, alanine aminotransferase; AST, aspartate aminotransferase; BDR, bile ductular reaction; BH, Benjamini-Hochberg; CPM, counts per million; CXCL13, chemokine (C-X-C motif) ligand 13; DEGs, differentially expressed genes; HepPten−, hepatocyte-deletion of Pten−; IL10RA, interleukin 10 receptor A; IPA, Ingenuity Pathway Analysis; LefSe, linear discriminant analysis effect-size method; mNAS, modified NAFLD scores; NAS, NAFLD Activity Score; OTU, operational taxonomic units; PAM, partitioning around medoids; PCoA, principal coordinates analysis; qPCR, quantitative PCR; RDA, redundancy analysis; RNAseq, RNA sequencing; SCF, stem cell factor.

Acknowledgments

We would like to thank Dr Tina Chang, Lauren K. McDaniel, and Ivonne I. Flores from the MD Anderson Cancer Center Microbiome Core Facility for their help with stool specimen processing and sequencing.

Funding

This study was funded in part by the MD Anderson Cancer Center SPORE in Hepatocellular Carcinoma Grant P50 CA217674 from the National Cancer Institute (NCI) and the NCI Grant R01 CA204665, to Laura Beretta. The MD Anderson Microbiome Core Facility is supported by the MD Anderson Cancer Center Support Grant (CCSG) P30 CA016672 from NCI. The MD Anderson Cancer Center Advanced Technology Genomics Core was supported in part by P30 CA016672 and NIH 1S10OD024977. The MD Anderson Research Animal Support Facility was supported by the NIH/NCI under award number P20CA016672. Lorenzo Cruz was supported by the 1R25CA240137 UPWARDS Training Program and the CPRIT Research Training Award CPRIT Training Program (RP210028).

Disclosure

The authors have declared that no conflict of interests exists in this work.

References

1. Younossi ZM, Marchesini G, Pinto-Cortez H, Petta S. Epidemiology of nonalcoholic fatty liver disease and nonalcoholic steatohepatitis: implications for liver transplantation. Transplantation. 2019;103(1):22–27. doi:10.1097/TP.0000000000002484

2. Younossi ZM, Blissett D, Blissett R, et al. The economic and clinical burden of nonalcoholic fatty liver disease in the United States and Europe. Hepatol. 2016;64(5):1577–1586. doi:10.1002/hep.28785

3. Younossi ZM, Koenig AB, Abdelatif D, Fazel Y, Henry L, Wymer M. Global epidemiology of nonalcoholic fatty liver disease-Meta-analytic assessment of prevalence, incidence, and outcomes. Hepatol. 2016;64(1):73–84. doi:10.1002/hep.28431

4. Estes C, Anstee QM, Arias-Loste MT, et al. Modeling NAFLD disease burden in China, France, Germany, Italy, Japan, Spain, United Kingdom, and United States for the period 2016-2030. J Hepatol. 2018;69(4):896–904. doi:10.1016/j.jhep.2018.05.036

5. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2019. Ca a Cancer J Clinicians. 2019;69(1):7–34.

6. Peng C, Stewart AG, Woodman OL, Ritchie RH, Qin CX. Non-alcoholic steatohepatitis: a review of its mechanism, models and medical treatments. Front Pharmacol. 2020;11:603926. doi:10.3389/fphar.2020.603926

7. Schnabl B, Brenner DA. Interactions between the intestinal microbiome and liver diseases. Gastroenterology. 2014;146(6):1513–1524. doi:10.1053/j.gastro.2014.01.020

8. Puri P, Sanyal AJ. The intestinal microbiome in nonalcoholic fatty liver disease. Clin Liver Dis. 2018;22(1):121–132. doi:10.1016/j.cld.2017.08.009

9. Schwimmer JB, Johnson JS, Angeles JE, et al. Microbiome signatures associated with steatohepatitis and moderate to severe fibrosis in children with nonalcoholic fatty liver disease. Gastroenterology. 2019;157(4):1109–1122. doi:10.1053/j.gastro.2019.06.028

10. Adams LA, Wang Z, Liddle C, et al. Bile acids associate with specific gut microbiota, low-level alcohol consumption and liver fibrosis in patients with non-alcoholic fatty liver disease. Liver Int. 2020;40(6):1356–1365. doi:10.1111/liv.14453

11. Caussy C, Tripathi A, Humphrey G, et al. A gut microbiome signature for cirrhosis due to nonalcoholic fatty liver disease. Nat Commun. 2019;10(1):1406. doi:10.1038/s41467-019-09455-9

12. Dong TS, Katzka W, Lagishetty V, et al. A microbial signature identifies advanced fibrosis in patients with chronic liver disease mainly due to NAFLD. Sci Rep. 2020;10(1):2771. doi:10.1038/s41598-020-59535-w

13. Loomba R, Seguritan V, Li W, et al. Gut microbiome-based metagenomic signature for non-invasive detection of advanced fibrosis in human nonalcoholic fatty liver disease. Cell Metab. 2017;25(5):1054–1062e1055. doi:10.1016/j.cmet.2017.04.001

14. Oh TG, Kim SM, Caussy C, et al. A universal gut-microbiome-derived signature predicts cirrhosis. Cell Metab. 2020;32(5):878–888.e6. doi:10.1016/j.cmet.2020.06.005

15. Sharpton SR, Schnabl B, Knight R, Loomba R. Current concepts, opportunities, and challenges of gut microbiome-based personalized medicine in nonalcoholic fatty liver disease. Cell Metab. 2021;33(1):21–32. doi:10.1016/j.cmet.2020.11.010

16. Ma C, Han M, Heinrich B, et al. Gut microbiome-mediated bile acid metabolism regulates liver cancer via NKT cells. Science. 2018;360(6391). doi:10.1126/science.aan5931.

17. Ponziani FR, Bhoori S, Castelli C, et al. Hepatocellular carcinoma is associated with gut microbiota profile and inflammation in nonalcoholic fatty liver disease. Hepatol. 2019;69(1):107–120. doi:10.1002/hep.30036

18. Sydor S, Best J, Messerschmidt I, et al. Altered microbiota diversity and bile acid signaling in cirrhotic and noncirrhotic NASH-HCC. Clin Transl Gastroenterol. 2020;11(3):e00131. doi:10.14309/ctg.0000000000000131

19. Fang HJ, Shan SB, Zhou YH, Zhong LY. Diabetes mellitus and the risk of gastrointestinal cancer in women compared with men: a meta-analysis of cohort studies. BMC Cancer. 2018;18(1):422. doi:10.1186/s12885-018-4351-4

20. Brandon-Warner E, Walling TL, Schrum LW, McKillop IH. Chronic ethanol feeding accelerates hepatocellular carcinoma progression in a sex-dependent manner in a mouse model of hepatocarcinogenesis. Alcohol Clin Exp Res. 2012;36(4):641–653. doi:10.1111/j.1530-0277.2011.01660.x

21. Naugler WE, Sakurai T, Kim S, et al. Gender disparity in liver cancer due to sex differences in MyD88-dependent IL-6 production. Science. 2007;317(5834):121–124. doi:10.1126/science.1140485

22. Keng VW, Sia D, Sarver AL, et al. Sex bias occurrence of hepatocellular carcinoma in Poly7 molecular subclass is associated with EGFR. Hepatol. 2013;57(1):120–130. doi:10.1002/hep.26004

23. Liu L, Yu K, Huang C, et al. Sex differences in hepatocellular carcinoma indicated BEX4 as a potential target to improve efficacy of lenvatinib plus immune checkpoint inhibitors. J Cancer. 2022;13(11):3221–3233. doi:10.7150/jca.73051

24. Bhat M, Pasini E, Pastrello C, et al. Estrogen receptor 1 inhibition of Wnt/beta-Catenin signaling contributes to sex differences in hepatocarcinogenesis. Front Oncol. 2021;11:777834. doi:10.3389/fonc.2021.777834

25. Natri HM, Wilson MA, Buetow KH. Distinct molecul