Identification and characterization of sexual dimorphism‑linked gene expression profile in hepatocellular carcinoma.

Hepatocellular carcinoma (HCC) is sexually disparate in humans, with a significantly increased prevalence in males. The molecular mechanisms by which the inhibition or development of liver cancer are facilitated require further investigation with regard to sex factors affecting disease progression. In the present study, functional signatures of differentially expressed genes (DEGs) were screened in female and male tumors via bioinformatics analysis. The following gene chip expression profiles were downloaded from the Gene Expression Omnibus: GSE19665, GSE23342 and GSE9843. They comprised cancerous and non‑cancerous tissue from patients with HCC and included critical sex features. Further evaluation of selected DEGs in the two sexual groups was performed via hierarchical clustering analysis. Venn diagram and functional protein‑protein interaction (PPI) network analyses were performed. Survival analysis of patients with differences in gene expression levels was subsequently performed using the Kaplan‑Meier Plotter database. Certain identified DEGs were common in female and male tumor samples, whereas others exhibited a sexually‑biased expression profile. Gene Ontology revealed that the cell cycle module 'biological process' was enriched in tumors derived from both sexes, whereas the metabolic pathways and drug metabolism modules were only significantly enriched in cancer tissues from male subjects. A number of hub DEGs in the cell cycle and p53 signaling pathways were involved in significant protein‑protein interaction (PPI) modules, including CDK1 and CCNB1. These DEGs were upregulated in tumors derived from female subjects compared with those derived from male subjects, and could be used as markers of poor prognosis in male patients. Other genes, such as CYP3A4 and SERPINA4, were identified in metabolic pathways, and were downregulated in male compared with female subjects. These genes were associated with a decreased survival rate. The data demonstrated that sex differences in physiology may regulate the levels of gene expression and/or activity, including gene function associated with oncogenesis and the outcomes of liver cancer. Additional surveys are required to explore in detail the molecular mechanisms underlying the differences in gene expression between the two sexes during the development of liver cancer.

Introduction

Sexual disparity has long been recognized as an independent key factor that can affect susceptibility to, and the incidence, development and management of various types of diseases, including vascular and heart disease, brain disorders and cancer (1–3). The susceptibility to cancer and the risk of malignancy are generally much higher in male compared with in female subjects (4–6). Therefore, sex is considered a prognostic indicator for various cancers (7); however, the endogenous molecular signatures of sex-based differences in different types of cancer require further investigation.

Liver cancer is the 5th most common type of cancer and the second leading cause of mortality in males globally, whereas it is the 7th most common cancer and the 6th leading cause of mortality in women (3,8). Hepatocellular carcinoma (HCC) accounts for 75–90% of primary liver cancers. The 5-year survival of patients with HCC is only 7%, and the majority of the cases result from risk factors, such as hepatitis B or hepatitis C virus infection (9). A previous study reported that the incidence of HCC is affected by sex (10) and predominant in males, with the male:female ratio within populations usually ranging between 2:1 and 7:1 (11–13).

In human and rodent females, differences in the incidence of HCC have been attributed to the expression of sexual hormones and associated receptors, such as estrogen and the estrogen receptor (ER) (14–16). For example, estrogen exerted protective effects against chemically-induced HCC tumorigenesis by inhibiting IL-6 production in macrophage-Kupffer cells, which in turn reduced tumor proliferation (16). However, the expression levels of IL-6 did not necessarily correlate with tumor load in female and male mice when forkhead box protein A (Foxa)1 and 2 is ablated in hepatocytes (6). The steroid hormone receptors ER and androgen receptor co-regulate Foxa factors, Foxa1 and Foxa2, which are transcription factors involved in carcinogen-induced HCC development or proliferation (6,17); however, no mutations in Foxa1/2 have been reported in humans. Comprehensive signatures of significantly sex-affected pathways, including immune response, cell cycle, metabolism-related, and DNA repair and p53 signaling pathways have been identified in numerous cancer types (10); however, sex-dependent gene expression signatures in liver cancer remain unclear (18).

The goal of the present study was to identify differences in the transcription of specific genes identified in the tumors of male and female subjects using bioinformatic tools (Fig. 1) (19). Microarray technology and bioinformatic analysis have been increasing regarded as a useful approach to detect differentially expressed genes (DEGs) and underlying molecular functional pathways involved in tumor carcinogenesis and progression (20). At present, no detailed bioinformatic analysis has been performed on number of microarrays of liver cancer and non-cancerous tissues which have been uploaded to the Gene Expression Omnibus (GEO) (20,21). Therefore, in the present study, Gene Ontology (GO) analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment, and protein-protein interaction (PPI) network and survival analyses were performed systematically to characterize gene expression profiles downloaded from GEO, identify genes associated with the effects of sex on carcinogenesis and progression, and emphasize the importance of sex differences with respect to the available animal or clinical trial data for HCC. Additionally, expression data were obtained from The Cancer Genome Atlas (TCGA) to validate the sex-biased gene expression profiles.

Figure 1.

Work flow chart of the bioinformatics analysis performed on the gene chip HCC datasets. The chart contains three main steps: Identification of DEGs in HCC according to patient sex; identification of the gene functions and pathways differentially expressed between normal and tumor tissues according to sex; and survival analysis of sex-associated genes. DEG, differentially expressed gene; GO, Gene Ontology; HCC, hepatocellular carcinoma.

Materials and methods

Microarray data preparation

Gene expression profiles that provided information concerning the sex features of subjects, including GSE19665 (22), GSE23343 (23) and GSE9843 (24) were downloaded from the GEO (http://www.ncbi.nlm.nih.gov/geo), a public functional genomics dataset with high-throughput gene expression data chips and microarrays (21). The selection was based on a GeneChip Human Genome U133 Plus 2.0 Array platform (Affymetrix; Thermo Fisher Scientific, Inc.). The datasets comprised 108 tissue samples, including normal liver samples from 4 female subjects and 13 male subjects, in addition to 28 and 63 HCC samples from female and male subjects, respectively.

Furthermore, raw tissue mRNA expression data and corresponding clinical information from patients with HCC were downloaded from The Cancer Genome Atlas (TCGA; http://cancergenome.nih.gov). The selection included 22 and 28 paired samples of normal liver tissues and corresponding HCC samples from 117 male and 245 female subjects with HCC, respectively.

Data process in analysis and identification of sex-biased DEGs

The raw data were preprocessed and analyzed by utilizing the affy (version 1.54.0; http://www.bioconductor.org/packages/release/bioc/html/affy.html) package in R (version 3.5.1). Initially, a robust multi-array average algorithm was used to preprocess data, including background correction, quantile normalization and calculation of gene expression values, following conversion of raw data into an expression matrix. Subsequently, the t-test method in the Linear Models for Microarray Data (limma; version 3.32.10) package in R was used to calculate the P-values and identify the DEGs between the HCC and normal liver tissue samples (25–27). Quantified gene expression data obtained from TCGA datasets were normalized and analyzed via the same process using the limma package.

The Benjamini-Hochberg method (28) was applied to calculate adjusted P-values and the false discovery rate (FDR). Probe sets without corresponding gene symbols or genes with >1 probe set were removed or averaged, respectively. Only the genes with |log2 fold change| ≥1.5 and adjusted P-value ≤0.05 were selected as DEGs. Hierarchical clustering analysis and heatmap analysis of DEGs were performed using the heatmap (version 1.0.0) package in R (29). VennDiagram (version 1.6.20, http://CRAN.R-project.org/package=VennDiagram), a package to generate high-resolution Venn and Euler plots in R, was used to construct visual analyses of the common and exclusive DEGs in the subjects of the two different sexes.

GO and pathway enrichment analysis of DEGs

The Database for Annotation, Visualization and Integrated Discovery (DAVID, http://david.abcc.ncifcrf.gov/) is a public biological resource regularly utilized for functional annotation and pathway analysis (21). In order to fully understand the biological themes of DEGs, GO (30,31) and KEGG pathway enrichment analyses (32–34) were performed using DAVID in order to examine the functions of these DEGs. A P<0.05 and a count >2 were considered as the threshold for significant differences. The top 5 genes derived from functional enrichment of each GO subset and KEGG pathway in female and male groups were depicted in a bubble diagram by the ggplot2 (version 3.1.1) package in R (35).

PPI network construction of DEGs

The STRING database (version 10.0) is an online biological resource used to detect PPIs and identify precise functions of proteins (36). Protein interactions were detected using a cut-off combination threshold of 0.4 using STRING. The data were further downloaded to construct a PPI network using Cytoscape (version 3.5.0) (37). Significant modules with stronger PPIs were selected in the plugin Molecular Complex Detection (MCODE) (38) using the following default parameters: Degree cut-off ≥2, node score cut-off ≥2, K-core ≥2 and maximum depth=100. P<0.05 was considered to indicate a statistically significant difference.

Survival analysis of sex-biased genes

The Kaplan-Meier Plotter database (http://kmplot.com/analysis) is handled by a PostgreSQL server, which integrates gene expression and clinical data simultaneously, consisting of 364 patients from GEO, European Genome-Phenome Archive (EGA; http://ega-archive.org/) and TCGA datasets (39). This database was used to perform survival analysis for identified DEGs based on the hazard ratio index, the 95% confidence interval and a log-rank P<0.05.

Results

Identification of DEGs

A total of 54,675 genes were obtained from the integrated gene chips from the GEO datasets. A total of 856 and 496 DEGs were identified in tumors from female and male subjects, respectively. The DEGs identified in the female patients included 616 upregulated and 240 downregulated genes, whereas the DEGs identified in the male subjects comprised 173 upregulated and 323 downregulated genes. The heatmaps of the bidirectional clustering suggested significant differences between the two sexes (Fig. 2A and B). Furthermore, Venn diagram analysis of DEG characterization between female and male subjects (Fig. 2C) revealed that 600 DEGs were identified only in female subjects, and 240 DEGs were found only in male subjects. A total of 505 and 95 genes were upregulated and downregulated, respectively, in female subjects, whereas 179 and 61 genes were downregulated and upregulated, respectively, in male subjects. A total of 255 genes, containing 111 upregulated genes and 144 downregulated genes, were common between female and male patients with HCC. In addition, DEGs obtained from GEO and TCGA samples intersected with each other according to sex, validating the aforementioned results (Fig. S1). The intersection contained almost all of the DEGs identified in males, and ~40% of the DEGs identified in females from GEO analysis, respectively.

Figure 2.

Analysis of sex-dependent differential gene expression in hepatocellular carcinoma. Heatmaps of DEG clustering in (A) female and (B) male patients. The top red bars indicate normal liver tissues, whereas blue bars indicate tumor samples. Green boxes represent downregulated genes, and red boxes represent upregulated genes. (C) Venn diagram of DEGs in the two sexes. Red and blue represent upregulated and downregulated DEGs in female subjects, respectively; green and yellow represents upregulated and downregulated DEGs in male subjects, respectively. DEG, differentially expressed gene.

Functional and pathway enrichment analysis of DEGs

All the identified DEGs from the female and male groups in the GEO datasets were uploaded separately into DAVID to identify the functional enrichment of the representative DEGs. A total of 220 GO terms for the female DEGs and 188 GO terms for the male DEGs were obtained. The top 10 terms in the three GO categories, including ‘biological process’ (BP), ‘cellular component’ (CC) and ‘molecular function’ (MF), were presented in Fig. 3.

Figure 3.

Gene Ontology terms and KEGG pathway enrichment of differentially expressed genes separated by sex. (A and D) Biological processes, (B and E) cellular components and (C and F) molecular functions, and (G and H) KEGG pathway analysis in female and male cohorts, respectively. KEGG, Kyoto Encyclopedia of Genes and Genomes; Rich factor, enrichment factor. Gene Ontology terms and KEGG pathway enrichment of differentially expressed genes separated by sex. (A and D) Biological processes, (B and E) cellular components and (C and F) molecular functions, and (G and H) KEGG pathway analysis in female and male cohorts, respectively. KEGG, Kyoto Encyclopedia of Genes and Genomes; Rich factor, enrichment factor. Gene Ontology terms and KEGG pathway enrichment of differentially expressed genes separated by sex. (A and D) Biological processes, (B and E) cellular components and (C and F) molecular functions, and (G and H) KEGG pathway analysis in female and male cohorts, respectively. KEGG, Kyoto Encyclopedia of Genes and Genomes; Rich factor, enrichment factor. Gene Ontology terms and KEGG pathway enrichment of differentially expressed genes separated by sex. (A and D) Biological processes, (B and E) cellular components and (C and F) molecular functions, and (G and H) KEGG pathway analysis in female and male cohorts, respectively. KEGG, Kyoto Encyclopedia of Genes and Genomes; Rich factor, enrichment factor.

The BPs of the DEGs from the female subjects were primarily enriched in mitotic nuclear division, sister chromatid cohesion, cell division, mitotic cytokinesis, DNA replication and mitotic cell cycle checkpoint. The BPs of the DEGs from the male subjects were mainly enriched in the epoxygenase P450 pathway, the exogenous drug catabolic process, the oxidation-reduction process, the drug metabolism process and the negative regulation of cell growth (Fig. 3A and D).

The CC results revealed that DEGs in the female cohorts were mainly enriched in nucleoplasm, cytosol, nucleus, cytoplasm and membrane components. The DEGs in the tumor samples from the male cohorts were mainly enriched in the following modules: Extracellular region, extracellular exosome, extracellular space, organelle membrane and blood microparticle (Fig. 3B and E).

The MF results indicated that DEGs identified in the female cohorts were primarily enriched in protein binding, ATP binding, poly(A) RNA binding, chromatin binding, protein homodimerization activity and microtubule binding. The results indicated that the DEGs in the samples from males were mainly enriched in oxidoreductase activity, oxygen binding, heme binding, iron ion binding and monooxygenase activity (Fig. 3C and F).

KEGG signaling pathway analysis was used to identify the top five enriched pathways in female subjects, which included ‘Pathways in cancer’, ‘Cell cycle’, ‘Viral carcinogenesis’, ‘Biosynthesis of antibiotics’ and ‘Oocyte meiosis’, whereas in the male subjects, the top enriched pathways included ‘Metabolic pathways’, ‘PI3K-Akt signaling pathway’, ‘Cell cycle’, ‘Retinol metabolism’, ‘Transcriptional dysregulation in cancer’ and ‘Drug metabolism-cytochrome P450’ (Fig. 3G and H).

PPI network analysis of DEGs-significant modules and core genes in network

DEGs were imported into STRING to analyze the PPIs between the various gene targets. PPI networks were constructed in Cytoscape (Fig. 4). Common hub genes between the two sexes were selected with a cut-off degree >10, and the ranked top 15 female DEGs included TOP2A, CDK1, GAPDH, CCNB1, ACLY, CCNB2, BIRC5, NDC80, CCNA2, CDC20, CDKN3, MAD2L1, BUB1, AURKA and KIF11. The top ranked male DEGs comprised TOP2A, CDK1, CCNB1, CDKN3, BIRC5, AURKA, CCNB2, MAD2L1, HMMR, TTK, EZH2, CDC20, PTTG1, RACGAP1 and NCAPG (Table I). The hub genes in each group were further analyzed for pathway enrichment, which exhibited a certain degree of similarity with the top five pathways; of note, the pathways ‘cell cycle’, ‘p53 signaling pathway’ and ‘oocyte meiosis’ were included in the two sex groups (Table II).

Figure 4.

Significant submodules of protein-protein interactions. (A) Module A, (B) module B, (C) module C and (D) module D were selected due to them exhibiting Molecular Complex Detection scores >4. Red represents upregulated nodes and green represents downregulated nodes.

Table I.

Top hub genes in the protein-protein interaction network with regard to their node degree between female and male.

A, Female
Gene name	Node degree	Clustering coefficient
TOP2A	156	0.20239
CDK1	114	0.36625
GAPDH1	111	0.11433
CCNB1	107	0.41474
ACLY	96	0.08794
CCNB2	91	0.52454
BIRC5	89	0.49413
NDC80	88	0.54075
CCNA2	87	0.53034
CDC20	86	0.53844
CDKN3	86	0.54118
MAD2L1	85	0.57843
BUB1	84	0.60614
AURKA	82	0.56218
KIF11	79	0.63746
B, Male
Gene name	Node degree	Clustering coefficient
TOP2A	78	0.26607
CDK1	54	0.53948
CCNB1	52	0.57919
CDKN3	47	0.66605
BIRC5	46	0.68213
AURKA	46	0.69952
CCNB2	44	0.75159
MAD2L1	44	0.75370
HMMR	44	0.72093
TTK	43	0.78295
EZH2	42	0.60511
CDC20	41	0.80610
PTTG1	41	0.80488
RACGAP1	41	0.83171
NCAPG	41	0.85854

Table II.

Top five enriched pathways of the hub genes in female and male subjects.

A, Female
Term	P-value	Genes
cfa04110:Cell cycle	6.24×10−9	CCNB1, CDK1, CCNB2, MAD2L1, BUB1, CDC20, CCNA2
cfa04914:Progesterone-mediated oocyte maturation	7.00×10−8	CCNB1, CDK1, CCNB2, MAD2L1, BUB1, CCNA2
cfa04114:Oocyte meiosis	1.14×10−5	CDK1, MAD2L1, BUB1, CDC20, AURKA
cfa04115:p53 signaling pathway	3.87×10−3	CCNB1, CDK1, CCNB2
cfa05203:Viral carcinogenesis	3.09×10−2	CDK1, CDC20, CCNA2
B, Male
Term	P-value	Genes
bta04110:Cell cycle	8.73×10−11	CCNB1, CDK1, MAD2L1, CCNB2, TTK, CDC20, PTTG1
bta04114:Oocyte meiosis	2.77×10−7	CDK1, MAD2L1, AURKA, CDC20, PTTG1
bta04914:Progesterone-mediated oocyte maturation	2.33×10−4	CCNB1, CDK1, MAD2L1, CCNB2
bta04115:p53 signaling pathway	4.50×10−3	CCNB1, CDK1, CCNB2
bta05166:HTLV–I infection	5.89×10−3	MAD2L1, CDC20, PTTG1

HTLV-1, human T-cell leukemia virus type 1.

GO and KEGG enrichment analysis of DEGs from subgroup analysis of the PPI network

According to the cut-off criteria, several significant modules were extracted from the male and female PPI network using the MCODE plugin for Cytoscape (Fig. 4). Module A from the female PPI network was constructed with 58 nodes and 1,564 edges, and contained a high number of hub genes (Fig. 4A). Module B in the male PPI network was constructed with 37 nodes and 651 edges (Fig. 4B). Module C in the male PPI network was constructed with 14 nodes and 59 edges (Fig. 4C). Module D in the male PPI network was constructed with 7 nodes and 14 edges (Fig. 4D). Notably, a large number of DEGs that were upregulated in the female module A were also found in the male module B. Additionally, the genes CYP3A4 and SERPINA4 were included in the male downregulated gene modules C and D (Fig. 4C and D).

To further explore the detailed functions of the genes in each module, the DEGs in the PPI networks were analyzed using GO, KEGG and STRING. Cell cycle, p53 signaling and oocyte meiosis pathways were common in the two sexes. Additional metabolic processes that were primarily enriched in the male sex, including retinol metabolism, drug metabolism and linoleic acid metabolism, were identified (Table III).

Table III.

GO and KEGG enrichment analyses of DEGs as determined by survival analysis.

A, Module A
Pathway description	P-value	Matching proteins
Cell division	2.38×10−41	ANLN, ASPM, AURKA, BIRC5, BUB1, BUB1B, CCNA2, CCNB1, CCNB2, CDC20, CDC6, CDC7, CDK1, CENPF, CEP55, CKAP2, CKS2, ECT2, HELLS, KIF11
Mitotic nuclear division	4.14×10−30	ANLN, ASPM, BIRC5, BUB1, CCNA2, CCNB1, CCNB2, CDC20, CDC6, CDK1, CENPF, CEP55, DLGAP5, HELLS, KIF14, KIF4A, MAD2L1
Oocyte meiosis	2.67×10−6	AURKA, BUB1, CCNB1, CCNB2, CDC20, MAD2L1, PTTG1, CDK1, AURKA
p53 signaling pathway	4.17×10−2	CCNB1, CDK1, CCNB2, RRM2
B, Module B
Pathway description	P-value	Matching proteins
Cell cycle	5.42×10−7	CCNB1, CDK1, BUB1B, CCNB2, CDC20, MAD2L1, PTTG1, TTK
Oocyte meiosis	4.64×10−6	CDK1, AURKA, CCNB1, CCNB2, CDC20, MAD2L1, PTTG1
Progesterone-mediated oocyte maturation	2.25×10−2	CCNB1, CDK1, MAD2L1, CCNB2
p53 signaling pathway	1.65×10−2	CCNB1, CDK1, CCNB2, RRM2
HTLV–I infection	7.93×10−3	BUB1B, CCNB2, CDC20, MAD2L1, PTTG1
C, Module C
Pathway description	P-value	Matching proteins
Retinol metabolism	2.15×10−15	CYP1A2, CYP26A1, CYP2A6, CYP2B6, CYP2C8, CYP2C9, CYP3A4, CYP4A11
Drug metabolism - cytochrome P450	9.02×10−13	CYP1A2, CYP2A6, CYP2B6, CYP2C19, CYP2C8, CYP2C9, CYP3A4
Chemical carcinogenesis	3.96×10−10	CYP1A2, CYP2A6, CYP2C19, CYP2C8, CYP2C9, CYP3A4
Linoleic acid metabolism	4.69×10−10	CYP1A2, CYP2C19, CYP2C8, CYP2C9, CYP3A4
Metabolic pathways	1.56×10−8	ACSL4, ACSL6, CYP1A2, CYP2A6, CYP2B6, CYP2C19, CYP2C8, CYP2C9, CYP3A4, CYP4A11
D, Module D
Pathway description	P-value	Matching proteins
Negative regulation of peptidase activity	3.92×10−3	ECM1, HGF, SERPINA4, SPP2
Negative regulation of endopeptidase activity	2.06×10−2	ECM1, HGF, IGF1, SERPINA4, SPP2
Negative regulation of macromolecule metabolic process	4.80×10−3	ECM1, FGFR2, HGF, IGF1, SERPINA4, SPP2
Negative regulation of cellular metabolic process	8.95×10−3	ECM1, FGFR2, HGF, IGF1, SERPINA4, SPP2
Regulation of molecular function	4.80×10−3	ECM1, FGFR2, HGF, IGF1, SERPINA4, SPP2

DEG, differentially expressed gene; GO, Gene Ontology; HTLV-1, human T-cell leukemia virus type 1; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Survival analysis of candidate genes in PPI modules

To explore the effects of each module on patient survival, the prognostic value of DEGs for each sex was analyzed using the Kaplan-Meier Plotter and data derived from the combination of the GEO, EGA and TCGA databases. First, common hub genes in modules A and B that were shared in specific common pathways for the two sexes, such as CDK1 and CCNB1, were analyzed. The expression of these two genes and other shared genes were not significantly associated with the survival of female patients (Fig. 5A and C; other genes not shown). Conversely, significantly decreased survival was observed for male patients with overexpression of these hub genes (Fig. 5B and D; other genes not shown).

Figure 5.

Sex-specific survival analysis for differentially expressed genes in patients with hepatocellular carcinoma. Survival analysis was performed for (A and B) CDK1, (C and D) CCNB1, (E and F) CYP3A4 and (G and H) SERPINA4 in females and males, respectively, using the Kaplan-Meier Plotter database. P<0.05 was considered statistically significant. CCNB1, cyclin B1; CDK1, cyclin-dependent kinase 1; CYP3A4, cytochrome P450 3A4; SERPINA4, serpin family A member 4.

Furthermore, genes identified in the subgroup PPI analysis that were specific to the male modules C and D were analyzed; it was revealed that CYP3A4 and SERPINA4, which were downregulated in male patients, were significantly associated with poor patient survival (Fig. 5F and H); however, no statistical significance was observed for the association between the expression of these genes and survival in female patients (Fig. 5E and G).

Discussion

HCC is a high-mortality malignancy with inadequate current medical treatment. Males usually exhibit an increased risk of developing HCC and poorer prognosis than females (3). The sex differences in physiological liver gene expression have been widely studied in rodents and humans (40–42); however, knowledge of sex-dependent gene expression in human HCC remains scarce, and may be of substantial biological and medical importance. In the present study, several mRNA profile datasets from GEO and TCGA were integrated using bioinformatic analysis methods to further identify sex-associated characteristics of gene expression between HCC tissues from female and male patients.

In the present study, liver samples from public datasets were obtained, and the following sex-based gene expression signatures were identified: A total of 858 DEGs from the female patient group, of which 111 were upregulated, and 144 were downregulated genes shared with the male cohort. A high number of DEGs were identified in only one sex. The study highlighted common molecular features for both sexes that were present in HCC tissue. Subsequent analysis into differences in the gene function, GO ontology and pathway analysis of DEGs between the two sexes were performed. DEGs in the female cohort were mainly enriched in the intracellular components, cell cycle, oocyte, mitotic nuclear division and cell division. DEGs in the male cohort were primarily enriched in extracellular components, metabolic pathways, epoxygenase P450 pathway, exogenous drug catabolic processes, oxidation-reduction processes and drug metabolism processes.

A majority of functional pathways were similar between the two sexual groups, as demonstrated by analysis of the common hub genes in the significant PPI networks and modules. It is well known that the cell cycle serves an important role in controlling cell proliferation, and that the dysregulation of the mitotic cycle accelerates carcinogenesis and tumor progression (43–45). Similar findings were noted in the present study; the cell cycle module was enriched in both sexes. Highly expressed hub genes, such as CDK1 and CCNB1, were shared between both sexes. These common genes were enriched in pathways such as the cell cycle pathway. Previous studies have reported that the cell cycle is a sex-affected process (10,46,47). The present study indicated that upregulated expression of CDK1 and CCNB1 was associated with poorer patient survival in males, whereas no significant differences were observed in the expression levels of these genes in female patients with HCC, suggesting that sex may affect the functional activity of genes involved in the cell cycle in HCC. Cyclin-dependent kinase 1 (CDK1) is required for regulating cell cycle progression, DNA replication and segregation, cell maturation and proliferation (48,49). The aberrant activation of CDK1 leads to a significant contribution to tumorigenesis by promoting cell proliferation (49). Cyclin B1 (CCNB1) is a member of the highly conserved cyclin family of proteins that is required for appropriate cell cycle progression (50). CCNB1 binds to CDK and forms an active complex known as maturation-promoting factor, which facilitates cell entry to mitosis (48,51). Overexpression of either CDK1 or CCNB1 is associated with advanced stage, poor differentiation, and poor prognosis in HCC (52,53). A previous study suggested that knockdown of CDK1 is more effective in males with HCC (19), consistent with the results of the present study. Previous studies demonstrated that 17 beta-estradiol exerts inhibitory and apoptotic effects on SW-13 human adrenal carcinoma cells in vitro by decreasing CCNB1 expression and blocking G2/M entry, highlighting that the expression of CCNB1 can be significantly decreased by estrogen (54). Therefore, high concentrations of estradiol may aid the protection of female patients against the development or progression of HCC by inhibiting activation of the cell cycle and subsequently carcinogenesis. Similarly upregulated expression levels of the cell cycle genes were detected between the two sexes; different levels of estrogen or ER expression may affect the activity of target genes and subsequently the prognosis of HCC.

The liver is an important organ in metabolizing and detoxifying dietary compounds and drugs. Sex-associated differences in the metabolism of drugs may account for liver toxicity; it has been shown that females exhibit increased susceptibility to drugs and associated adverse effects (55,56). In the present study, the findings indicated that DEGs in drug and metabolic pathways were enriched in the tumors of the male subjects. The majority of these genes were downregulated, suggesting that metabolic pathways, the epoxygenase P450 pathway, the exogenous drug catabolic process, the oxidation-reduction process, the oxidoreductase activity and drug metabolism processes were primarily downregulated in the tumors of male patients with HCC. Previous studies reported that oxidoreductase activity is an important component of antioxidant defenses, and that it acts as a tumor suppressor that is frequently altered in tumors (57,58). The relationship between the expression of candidate genes in significant modules and patient survival was investigated. Downregulation of the CYP3A4 gene was a significant predictor of male patient survival. Cytochrome P450 3A4 (CYP3A4) is mainly expressed in the liver and intestine, and is involved in the metabolic activation and metabolism of several pro-carcinogens (59,60). Sex-associated differences were observed in CYP3A4 enzyme activity towards different substrates with respect to their metabolic clearance (55). This enzyme exhibits higher activity in women than in men, and may be regulated by steroid hormones (55). A previous study reported that the downregulation of CYP3A4 may be associated with lower metabolic clearance of various procarcinogens, leading to higher accumulation of these compounds in the body (55); as a result, carcinogenesis is facilitated in patients with low levels of the CYP3A4 gene (55,59). Serpin family A member 4 (SERPINA4), also known as kallistatin, is a novel antiangiogenesis protein, which exerts various effects on inflammation, angiogenesis and tumor growth (61,62). Previous studies reported that SERPINA4 suppressed tumor growth by directly inhibiting cancer cell proliferation, migration and invasion, and by adjusting cancer cell signaling in various malignant tumors, including lung and breast cancers, and HCC (61–63). Kallistatin levels were positively correlated with free androgen index in women with polycystic ovary syndrome (64). In the present study, SERPINA4 was included in the analysis of the enrichment pathways and was associated with negative regulation in the metabolic process component. These findings may explain how low SERPINA4 gene expression in male patients with HCC may lead to poor survival.

Comparison of the gene expression levels between the two sexes raises questions regarding the mechanisms that confer additional protection in one sex compared with the other. The differences in physiology between the two sexes result primarily from the combination of the cell-autonomous effects of sex chromosomes, and from the activation effects of male and female sex hormones on their receptors (65). Previous studies revealed that the cell cycle and metabolic pathways are partially affected by steroid hormones and their receptors (55,66). The candidate genes or pathways identified in the present study require further investigate to assess their potential associations with estrogen or androgen and their receptors in an animal model of liver cancer.

Human sex-biased differences in hepatic genes at the genetic, regulatory and functional level can be associated with the hepatic disease state. Zhang et al (67) evaluated microarray data from liver samples from Western Europeans following the removal of primary liver tumors. They reported that significant sex-biased genes could affect diverse physiological functions, such as zinc finger clusters and lysine demethylase for epigenetic modification, and CYP3A for hepatic drug metabolism. Yuan et al (10) performed a comprehensive analysis of molecular differences with respect to sex across 13 cancer types from TCGA. Numerous sex-biased genes were enriched in sex chromosomes and accompanied with sex-biased DNA methylation. These genes were associated with the immune response, cell cycle, metabolism-pathways, DNA repair and p53 pathways. Ma et al (19) analyzed DEGs between neoplastic and normal tissue for each sex in HCC samples from TCGA. They reported that pathways associated with lipid metabolism were only significantly dysregulated in male subjects. The present study contained not only TCGA samples, but also GEO expression profiles in order to analyze different gene expression between male and female patients with respect to tissue (neoplastic vs. normal). The analyzed data included expression profiles from Asians and Europeans. Sex-biased genes that were, consistent with previous findings, associated with drug and metabolic pathways, as well as shared genes enriched in cell cycle and p53 pathways, were associated with patient prognosis in a sex-dependent manner.

The present study exhibited several limitations. First, the lack of clinical information for patients downloaded from GEO prevents adjustments for confounding factors, such as age, disease state, tumor purity or survival time, and may reduce the observation power. It should be noted that the bias of tumor and patient characteristics between male and female patients for which samples were downloaded from TCGA in the present study were not significantly different to those observed in a previous study (10). Second, the normal liver tissues utilized in the present study were primarily obtained from either non-tumor-bearing donors or liver samples adjacent to tumors. The tissues may consist of distinct cell types to corresponding tumor tissues, potentially confound comparison between normal and tumor samples. Additional patients were enrolled from TCGA to confirm the results from the GEO analysis; however, the majority of the tissues were obtained from patients with primary HCC tumors removed. Furthermore, the small number of normal tissues used in the study further sensitizes the study to potential confounds, limiting the detection power. Therefore, future studies into this topic will require rigorous analyses of large patient cohorts with greater control of clinical confounding variables. Third, the presented results were entirely based on microarray data. It remains to be established to what extent the observed sex differences in mRNA expression would contribute to sex differences at the protein level and the levels of biological activity in vivo. Finally, our study did not explore the sex differences in gene expression signatures and function derived from non-RNA-based mechanisms, such as translational regulation, protein function and post-translation modifications. Associations of sex-determined factors, such as sex chromosomes, sex hormones and their activation receptors, with sex-biased gene expression characteristics require further elucidation. In addition, the impact of sex-biased somatic alterations on the etiology, progression, treatment and personalized therapy of cancer should be considered. In summary, the present findings in sex-dependent DEGs and enriched functional pathways exhibited a certain degree of consistency with previous studies; however, the aforementioned limitations prevented a full characterization of sex-biased gene expression.

In conclusion, the present study screened sex dimorphism signatures of gene expression, gene functions and pathways in HCC using computational bioinformatic methods. Genes associated with the cell cycle and/or with the metabolic pathways were identified as dysregulated in both sexes; however, these genes exhibited distinct prognostic values for male and female patients with HCC. The differences in physiology noted between the two different sexes may account for the differences in the expression of genes, including those involved in oncogenesis and in liver cancer progression. Further studies into sex-biased DEGs in liver cancer at the genetic, regulatory and functional level should be conducted in animal models of HCC, or in larger independent cohorts with increased controlling for potential confounding factors, in order to further increase understanding of the role of these genes in hepatic disease states and therapeutic strategies.