The microarray datasets used in this study were obtained from NCBI Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) (8,9) and ArrayExpress (10). MESH terms ʽliver fibrosis', ʽmyofibroblasts' and ʽhepatic stellate cells' were used to identify potential datasets of interest, and GSE34949 and E-MTAB-2081 were selected for the analysis. Both datasets included whole-genome transcriptional data from rodent HSCs in the quiescent and activated state. For GSE34949, total RNA was extracted from primary C57BL/6 mouse HSCs on days 0 and 7 of culture. For the E-MTAB-2081 dataset, RNA samples were obtained from quiescent and activated HSCs derived from Sprague-Dawley rats. For each dataset, we performed moderate t-statistics to generate a list of upregulated genes. We used a threshold value of p<0.05 and a fold change of >2. In order to identify a robust fibrogenic gene signature, the web-based application, Venn Diagram Generator (http://www.bioinformatics.lu/venn.php), was used to select the differentially expressed genes (DEGs) commonly shared between the two datasets.

Functional analysis was performed on the selected genes using the Database for Annotation, Visualization and Integrated Discovery (DAVID) Bioinformatics Database (http://david.abcc.ncifcrf.gov/), as previously described (11). The biological significance of these genes was queried by performing enrichment analysis for Gene Ontology (GO) categories, including biological processes (BPs) and molecular functions (MFs), for domain architectures [Simple Modular Architecture Research Tool (SMART)] and for the identification of significant biological pathways. The Kyoto Encyclopedia of Genes and Genomes (KEGG), a collection of manually curated pathways and molecular networks, was used as reference. In order to construct a series of interconnected networks, the selected genes were subjected to k-means clustering using the STRING database (http://string-db.org/). STRING provides a global view of all the available interaction data by creating networks, which represent the current knowledge on the functional interconnections among genes. All associations are provided with a confidence score that represents an estimation on how likely a given association describes a functional linkage between two genes. Clustering was performed for five clusters and disconnected nodes were excluded. The web server, REVIGO (http://revigo.irb.hr), was used to summarize the lists of GO terms in order to find representative subsets of the terms that rely on semantic similarity.

We identified 1,800 genes which were upregulated in the activated vs. the resting HSCs in the GSE34949 dataset and 1,354 genes in the E-MTAB-2081 dataset. Among the significantly upregulated genes, 254 genes were common between the two datasets and were selected for further analysis (Fig. 1A). As expected, some genes had well-established associations with fibrosis. These included ECM components [collagen, type I, alpha 1 (COL1A1), collagen, type I, alpha 2 (COL1A2), collagen, type III, alpha 1 (COL3A1), collagen, type V, alpha 2 (COL5A2), collagen, type VIII, alpha 1 (COL8A1), collagen, type XII, alpha 1 (COL12A1), fibronectin 1 (FN1) and tenascin C (TNC)], metalloproteinases (MMP11, MMP12, MMP19 and MMP23) and growth factors [TGFB2, TGFB3, fibroblast growth factor (FGF)18 and FGF2].

We utilized DAVID to identify the major functions, networks and pathways relevant to the 254 significantly upregulated genes.

Pathway analysis revealed that the top three significant pathways enriched by the 254 genes were ʽfocal adhesion' (Figs. 1B and 3), ʽECM-receptor interaction' and ʽpathways in cancer'. These pathways shared several genes, i.e., the FGFs and ECM components already indicated. Other notable genes were represented by protein kinase C, alpha (PRKCA), insulin-like growth factor (IGF1), filamin C, gamma (FLNC), myosin, light chain 9, regulatory (MYL9), p21 protein (Cdc42/Rac)-activated kinase 3 (PAK3), platelet derived growth factor C (PDGFC) and c-fos-induced growth factor/vascular endothelial growth factor D (FIGF) in the ʽfocal adhesion' pathway, syndecan 2 (SDC2) in the ʽECM-receptor interaction' pathway and PRKCA, cyclin-dependent kinase inhibitor 2A (CDKN2A), cyclin-dependent kinase inhibitor 2B (CDKN2B), aryl-hydrocarbon receptor nuclear translocator 2 (ARNT2), IGF1, frizzled class receptor 2 (FZD2) and FIGF in the ʽpathways in cancer'. Among the other top enriched KEGG pathways, the presence of those genes related to hypertrofic and dilated cardiomyopathy was indicative of the myofibroblastic trans-differentiation of activated HSCs (Fig. 1B). In addition, the involvement of the p53, renal cell carcinoma, melanoma and mitogen-activated protein kinase (MAPK) signaling pathways seemed to correlate with liver fibrosis to liver cirrhosis and the evolution toward hepatocellular carcinoma (Fig. 1B), and as previously demonstrated (12).

The use of SMART led to the identification of the ʽcalcium-binding EGF-like domain', ʽEGF', ʽfibrillar collagens C-terminal domain', ʽzinc-depen dent metalloprotease' and ʽvon Willebrand factor, type C' as the principal protein domains overrepresented by the commonly upregulated genes during fibrogenesis (Fig. 2A and B).

A characteristic of HSCs is the ability to undergo cell cycle progression and cell mitosis (13,14). Indeed, biological processes strongly associated with the upregulated genes were found to be related to the regulation of cell proliferation, cartilage development, regulation of cell growth and division (Fig. 2B).

The most significant molecular functions (Fig. 2C) associated with the commonly upregulated genes were ʽcalcium ion binding', ʽgrowth factor binding', ʽextracellular matrix structural constituent' and ʽgrowth factor activity' (Fig. 2C).

The sole analysis of the gene expression profile provides an incomplete picture of the biological mechanisms involved in a process, such as liver fibrosis, since it does not include the interrelations and communication paths among genes. The construction of a gene network is the logical continuation of gene expression profiling. Describing a biological system using a network helps to identify characteristics that could potentially lead to the development of novel treatment strategies, by a better understanding of the process. In a gene network, genes form the network nodes and relationships serve as network edges.

Genes were clustered using the k-means algorithm for a number of five clusters (Fig. 4). The central hub was occupied by FN1 which interconnected all of the five clusters. Other important nodes were represented by caldesmon 1 (CALD1), COL1A1 and FGF2.

CALD1 is a protein expressed in smooth muscle cells, which combines with calmodulin, tropomyosin and actin thus regulates cellular contraction (15). Although it has been suggested to be a biomarker for smooth muscle tumors, its role in liver fibrosis remains unexplored.

Another important gene, at the interface of three different clusters, was represented by S100A4. S100A4, also known as fibroblast-specific protein 1 (FSP1), belongs to the S100 gene family, a multi-gene family of highly conserved Ca²⁺-binding proteins. S100A4 is known to be upregulated in pulmonary fibrosis (16), renal fibrosis (17) and in portal tract fibrosis following bile duct damage (18). Data on endometrial cancer have indicasted that S100A4 is induced by TGF-β1 and that it is necessary for the effects of TGF-β1 on cancer cell migration and invasion (19). S100A4 could thus play a central role in the induction and/or maintenance of liver fibrosis in chronic hepatic disease. The development of drugs specifically targeting S100A4 may thus be equally effective as anti-TGF-β treatment, but with less potential adverse effects, given the more downstream target.

Another gene that was identified as a relevant hub in the network analysis was Runt-related transcription factor 2 (Runx2). The Runx family of transcription factors is composed of three distinct genes, Runx1, Runx2 and Runx3, which play a key part in normal development and cancer (20). In particular, Runx2 is involved in osteoblast differentiation and bone formation and it has been reported that Runx2 is a downstream target of the TGF/Smad4 pathway. In addition, the MAPK pathway plays an essential role in the phosphorylation and activation of Runx2 (21). Hattori et al (22) demonstrated that the inhibition of p38 is able to ameliorate liver cirrhosis through the down-regulation of Runx2. Taken together, these data suggest the central involvement of this gene in the progression of hepatic fibrosis.

In order to define the complete mechanisms underlying the development of liver fibrosis, we identified commonly downregulated genes between the GSE34949 and E-MTAB-2081 datasets. We found 1,200 and 1,122 genes downregulated in these databases, respectively, with 132 genes being shared between the two datasets. GO analysis revealed five major biological processes enriched by these genes, represented by ʽresponse to wounding', ʽpositive regulation of developmental process', ʽcoagulation', ʽcarboxilic acid catabolism' and ʽpolysaccharide metabolism'. Fig. 5 shows a scatterplot of the cluster representatives obtained by applying multidimensional scaling to a matrix of the semantic similarities of the GO terms. Along the same lines, KEGG pathway analysis identified as the most enriched pathways, the ʽcomplement and coagulation cascades' [with kininogen 1 (KNG1), fibrinogen gamma chain (FGG), fibrinogen alpha chain (FGA), fibrinogen beta chain (FGB), coagulation factor II (thrombin) (F2), serpin peptidase inhibitor, clade C (antithrombin), member 1 (SERPINC1) and coagulation factor IX (F9) being involved], the ʽNOD-like receptor signaling pathway' [which included MAPK12, interleukin (IL)1B, nuclear factor of kappa light polypeptide gene enhancer in B-cells inhibitor, alpha (NFKBIA) and MAPK11], ʽcytokine-cytokine receptor interaction' [including pro-platelet basic protein (chemokine (C-X-C motif) ligand 7) (PPBP), IL1B, nerve growth factor receptor (NGFR), interleukin 13 receptor, alpha 1 (IL13RA1), IL10, IL1A and kinase insert domain receptor (KDR)] and the ʽRIG-I-like receptor signaling pathway' [with MAPK12, NFKBIA, MAPK11 and DEXH (Asp-Glu-X-His) box polypeptide 58 (DHX58)]. These data are concordant with those presented in the study by Ji et al (23) who found that downregulated genes were primarily related to immune response and lipid metabolism. Of note, GO analysis revealed that the major molecular functions involved were the endopeptidase inhibitor, the vitamin transporter, protein bridging and DNA-directed DNA polymerase activity.