Genes related to the very early stage of ConA-induced fulminant hepatitis: a gene-chip-based study in a mouse model

Background Due to the high morbidity and mortality of fulminant hepatitis, early diagnosis followed by early effective treatment is the key for prognosis improvement. So far, little is known about the gene expression changes in the early stage of this serious illness. Identification of the genes related to the very early stage of fulminant hepatitis development may provide precise clues for early diagnosis. Results Balb/C mice were used for ConA injection to induce fulminant hepatitis that was confirmed by pathological and biochemical examination. After a gene chip-based screening, the data of gene expression in the liver, was further dissected by ANOVA analysis, gene expression profiles, gene network construction and real-time RT-PCR. At the very early stage of ConA-triggered fulminant hepatitis, totally 1,473 genes with different expression variations were identified. Among these, 26 genes were finally selected for further investigation. The data from gene network analysis demonstrate that two genes, MPDZ and Acsl1, localized in the core of the network. Conclusions At the early stages of fulminant hepatitis, expression of twenty-six genes involved in protein transport, transcription regulation and cell metabolism altered significantly. These genes form a network and have shown strong correlation with fulminant hepatitis development. Our study provides several potential targets for the early diagnosis of fulminant hepatitis.


Background
Fulminant hepatitis mainly develops from chronic or acute hepatitis B virus (HBV) infection. It has been reported that world widely about 400 million people are living with chronic HBV infection. While 500,000 people die from HBV-related liver cirrhosis and hepatocellular carcinoma, acute fulminant hepatitis B causes additional 40,000 deaths each year [1].
Despite the significant improvement in the antiviral treatment targeting HBV infection, the morbidity and mortality of fulminant hepatitis keep at a high level. This is at least partly due to our poor understanding of the pathogenesis of the disease and the lack of early diagnosing/warning system which lead to a delayed medical intervention. As a serious illness, fulminant hepatitis is not only a threat to our health, but a big financial burden to the patient's family. Therefore, an accurate diagnosis at an early stage is linked directly to an effective treatment, better recovery and lowered mortality of the disease.
Because it resembles very much the human fulminant hepatitis, concanavalin A (ConA)-induced hepatitis has long been used as an experimental model of immunemediated liver disease [2]. It is characterized by massive hepatocellular degeneration and lymphoid infiltration of the liver [3]. The development of ConA-induced hepatitis depends mainly upon T cell activation [4] which raises plasma levels of various cytokines including TNF-α, interferon (IFN)-γ, and interleukin (IL)-6 [5]. To investigate the pathogenic mechanism of human fulminant hepatitis, in this study, we utilized the ConA-induced hepatitis model and selected microarray as our primary detection tool.
Recently, high-throughput gene chip has increasingly been used to explore the transcription dynamics on a genome-wide scale. This raises a possibility that similar approach could be used to identify gene expression at a specific condition such as the very early stage of fulminant hepatitis. Then the alternations in gene expression could possibly be used as diagnostic measures. Given this aim, it is necessary to use statistical methods to help identify genes whose transcript profiles respond to liver damage. In addition, such approach involves modeling the association of a generic response with a specific experimental variable, such as timing, tissues impacted, temperature, or drug dosage. Therefore, whole genome expression analysis provides a system-level insight for understanding crucial factors that control the liver damage and in turn may help us elucidate the molecular mechanisms of liver injury.
In this study, by using the ConA-induced mouse fulminant hepatitis model, we analyzed the hepatic gene expression with whole genome mouse GeneChip (Affymetrix). Through an integrative analysis that combines changes in gene expression with gene function within a genetic network, we identified 26 genes that may be related to ConA-induced liver injury.

ConA-injection induces fulminant hepatitis in Balb/C mice
Pathological examination at indicated time periods after ConA challenge showed that fulminant hepatitis mouse model was successfully made. As shown in Figure 1, a single injection of ConA led to a clear liver damage as presented by large area of necrosis ( Figure 1). In order to study the changes of very early stage of fulminant hepatitis, liver specimen were taken from the mice at 1 h, 3 h, and 6 h after ConA administration for the further investigation.

Biochemical evidence for ConA-induced hepatitis
To observe the liver damage in our ConA-induced hepatitis model, we measured by using the classical Reitman-Frankel method the serum levels of alanine aminotransferase (ALT), aspartate aminotransferase(AST), serum total bilirubin (STB) and serum lactate dehydrogenase (LDH). See details in Table 1.

Genes screened by Gene-chip study
Genome-wide transcriptional profiling of the liver has demonstrated that extensive gene expression occurs during very early-stage fulminant hepatitis. To investigate the possible gene expression change in our ConAinduced hepatitis model at a very early stage, we performed a gene chip study by using the Affymetrix probe dataset which includes totally 45101 probes. We found when compared with probes in the control group, 1559 probes increased and 1794 probes decreased in the 1 h group, 3040 probes increased and 4160 probes decreased in the 3 h group, and 5466 probes increased and 5445 probes decreased in the 6 h group. Our data indicate that by time after ConA injection more and more gene expression is detectable. The in detail gene expression data are also available at the GEO website under accession number GSE17184: http://www.ncbi.nlm.nih.gov/geo/

Genes further screened by ANOVA and model profile analysis
ANOVA corrected with the randomized variance model was performed to determine genes that were expressed separately and differentially. As a result, total 1,473 genes that had a p-value and the false discovery rate (FDR) less than 0.05 were declared to be significantly expressed. The gene expression value per group was the geometric mean of the Robust Multichip Average (RMA) normalized gene signals of 3 samples per time point.
To further narrow the target genes which harbor great significance among the declared 1,473 genes, we chose to use the twenty-six model profiles to summarize the expression pattern of the genes. As shown in Figure 2, among the 26 patterns, we identified 10 patterns of genes that shown very significant p-values (colored boxes).
Among these patterns, the two most significant patterns were profiles No. 9 and No.16 according to ascending p-values. While the profile No. 16 contained 118 genes whose expression increased constantly, the profile No. 9 contained 209 genes whose expression reduced constantly, after ConA injection ( Figure 3).  Genes in profiles No.9 and No.16 were then analyzed and identified by gene co-expression network with k-core algorithm to determine which gene or genes may play pivotal role in the early stage of fulminant hepatitis. Gene networks are constructed from functional gene associations. In the network, cycle nodes represent genes, and edges between two nodes represent interactions between genes, which were quantified by degree. Degrees within the network which describe the number of single gene that regulates other genes represent the size of the cycle node. The higher the degree, the more central the gene occurs within the network. The clustering coefficient can be used to estimate the complexity of interactions among genes that neighbor the core gene with the exception of core gene participation. The lower the clustering coefficient, the more independent of the core gene are the  interactions among genes in the neighborhood of the core gene [6]. A k-core of a gene co-expression network usually contains cohesive groups of genes. On the other hand, phylogenetic analysis classifies genes into groups based on the similarity of base sequences [7].
The k-core subnetwork with higher k-core level in our result is considered to have a core status within a largescale gene network made up of differential 26 genes (see Figure 4 and Table 2). Most of the genes are attributed to transport, transcription, regulation of transcription, metabolic processes, and carbohydrate metabolic processes, consistent with the gene ontology hierarchical category, lipoprotein catabolic processes, peptide biosynthetic processes, positive regulation of calcium-dependent cell-cell adhesion, fructosamine metabolic processes, male sex differentiation, and regulation of cholesterol biosynthesis. The core genes MPDZ and ACSL1 appear at the center of both the large-scale network and the 11 k-core subnetwork. They directly regulate 34 neighboring genes that interact according to their degrees. These interactions depend in large part on MPDZ and ACSL1 because the clustering coefficients of these genes are 0.38 and 0.42, which are lower than for other genes.
In this network, transmembrane proteins, such as MPDZ and Acsl1, were shown to play a role in the important processes of cell adhesion, carbohydrate metabolism and molecular transport. However, MPDZ expression constantly decreased over the four time points. It may contribute to damage to cell junctions in liver cells. The result of MPDZ down expression may be related to disorders of liver tissue and may lead to obstacles in the carbo-hydrate metabolic process. Acs11 was an essential gene that was associated with long-chain-fatty-acid-CoA synthase activity. The under expression of the gene gives rise to a series of obstructions in lipid and cholesterol biosynthetic and metabolic processes.

Real-time RT-PCR verification of the genes
Real-time RT-PCR was used to confirm the expression of the 26 genes. Compared with normal controls, the expression of most of the 26 genes showed a statistical significant change (p < 0.05). Expressions of PARD3, ALDH8A1, CMTM8 and FBXO3 were also higher than those of normal controls, but there was no statistical significance (p > 0.05) (See details in Table 3 and the primers see in Additional file 1).

Discussion
In this study, use of a fulminant hepatitis mouse model, we have focused on hepatic gene expression during early stages of fulminant hepatitis. By a series of biological and bioinformatics analysis, we identified 26 genes whose expression are significantly altered at the very early stage and have shown big correlations with fulminant hepatitis.
Using gene chips has the advantage of generating huge amounts of information. The high cost of microarray experiments and limited sample availability resulted in a small sample size for each treatment group. In addition, because we had to filter differential genes, our two-group comparison accuracy was reduced because variance estimates have few degrees of freedom. Consequently, we obtained results with more power and less variation by  using the random variance model corrected standard ttest or ANOVA methodology [8].
We predefined certain sets of profiles in the context of patterns that can be expected in the course of growth. On account of the noise and the small number of points for each gene, any one pattern will be shared by many genes and can be expected to appear by chance. First, we selected a set of potential expression profiles. These sets of profiles cover all of the possible profiles that can be generated by gene expression in the course of growth, and each represents a single temporal expression pattern. Enrichment of genes in each profile was then used to measure the profile significance level. Significant profiles may indicate that there were mostly common functions attributed to genes that are co-expressed. Such functions mainly influence the biological character [9].
In biological processes, macromolecular networks can be constructed by the experience results from such as Y-2H, and coimmunoprecipitation, or from the algorithmic predict base on the gene function correlation and expression profiles. Because of the flexibility of the network model base on the algorithmic predict from high throughput gene expression tests, we can look at snapshots of protein-protein interaction, gene expression regulatory networks and metabolism networks among different groups. The intrinsic gene networks of a phenotype can represent gene function propriety of sample. kcores and degrees of genes are key attributes in the network. A k-core of a network is a subgraph in which all genes are connected to at least k other genes in the subgraph [6]. In a protein-protein interaction network, a node represents a protein and an edge represents an interaction between proteins. In this work, we use k-cores Figure 4 Gene co-expression network. Genes from profiles No.9 and No.16 were analyzed and identified by gene co-expression network with kcore algorithm. Cycle nodes represent genes, the size of nodes represents the power of the interrelation among the nodes, and edges between two nodes represent interactions between genes. The more edges of a gene, the more genes connecting to it, the more central role it has within the network. of protein-protein interaction networks to define main gene function in main subgraph. The rank of k-core value describes the complexity of gene association relationship.
The complexity of gene relationship increased with kcore value rank. We wanted to find main Gene Ontology (GO) assigned by the maximum numbers of genes in separately k-core and then define the key gene functions at each complexity level of network [7]. For this analysis result, we conclude the core functions at the core status of network which have a top k-core level. Regulation of gene expression as well as protein-protein interactions can change in a dynamic process over time. Such processes may exhibit a regulatory character. These macromolecular networks respond to the events such as growth, metabolism, and apoptosis. Because certain gene interactions are only observable under very specific conditions, we combined experimental data to model coexpression networks from different types of experiments. One of the important reasons for doing this was the fact that differential macromolecular network structures have different sources. In this study, twenty-six genes were selected using the network, and validated by Real-time RT-PCR. Compared with normal controls, the expression of most of the 26 genes appeared a statistical significant change (p < 0.05).
A very important issue about the results from a genechip-based analysis is their correlation to the actual biological processes. Many of the 26 genes identified in this study including MPDZ and ACSL1 are involved in the metabolism of liver cells. Although the expression of these genes may not be specifically initiated or suppressed by fulminant hepatitis, the emergence or inhibition of the genes at the very early stage of fulminant hepatitis could be followed by severe liver damage, the pathological hallmark of fulminant hepatitis presented as large area of necrosis, thereby supplying potential diagnostic information for clinical fulminant hepatitis. In this study, we employed a ConA-induced animal fulminant hepatitis model instead of the use of human tissue which is extremely tough to obtain. As mentioned earlier, ConA-induced hepatitis has broadly been used as an experimental model of immune-mediated liver disease. Our pathological and biochemical data further confirmed that the model resembles very much the human fulminant hepatitis. Therefore, our results may provide important referential merit for clinical investigation. Nevertheless, in all conscience, the genes identified here are required to be further dissected and confirmed in human by other clinic-related studies.

Conclusions
At the early stages of fulminant hepatitis induced by ConA, expression of twenty-six genes involved in protein transport, transcription regulation and cell metabolism altered significantly. These genes form a network and have shown strong correlation with fulminant hepatitis development and a promising future for early diagnosis. Our study provides several potential targets for the early diagnosis of fulminant hepatitis in clinical scenario.

Animals
Eight-to 10-week-old male Balb/C mice were purchased from the Animal Center, Zhejiang Academy of Medical Sciences. The average weight was 25 g. All mice were kept at least 1 wk at 22°C and 55% relative humidity in a 12-h day/night lighting environment with free access to food and water. The animals were raised and cared for in accordance with guidelines established by the National Science Council of the Republic of China.

ConA administration
Mice were given a tail vein injection of ConA (Sigma Chemical Co., St. Louis, MO) in 100 μl PBS at a concentration of 20 mg/kg, and same dose of PBS was used as negative control. Degrees describe the number of single gene that regulates other genes represent the size of the cycle node. The higher the degree, the more central the gene occurs within the network. The clustering coefficient can be used to estimate the complexity of interactions among genes that neighbor the core gene with the exception of core gene participation. The lower the clustering coefficients, the more independent of the core gene are the interactions among genes in the neighborhood of the core gene. A k-core of a gene co-expression network usually contains cohesive groups of genes. The higher the k-core, the more central the genes occurs within the network. In this study, the top k-value is 11, which is considered to have a core status within a large-scale gene network made up of differential genes from profiles No.9 and No. 16.

Histology
The livers were removed, and part of each liver was frozen in liquid nitrogen, while another part was fixed with 4% phosphate-buffered paraformaldehyde and embedded in paraffin. Tissue sections (4 μm) were prepared and stained with hematoxylin/eosin; sections were then examined under light microscopy.

Biochemical detection
The extent of liver injury was assessed by determining serum ALT, AST, STB and LDH activity levels using the standard Reitman-Frankel method.

Isolation of RNA
Liver samples from very early stage of fulminant hepatitis mouse model, which include control group, 1 hour 3 hours and 6 hours after the injection of ConA, were ground into a fine powder in a mortar cooled by liquid nitrogen, and 100 mg was added to 1 ml prechilled Trizol reagent (Invitrogen, Carlsbad, CA). Total RNA extractions were performed according to the manufacturer's directions, and the RNA was further purified by passage through RNeasy mini-columns (QIAGEN, Valencia, CA) according to the manufacturer's protocols for RNA cleanup. Final RNA preparations were resuspended in RNasefree water and stored at -80°C. The RNA samples were quantified spectrophotometrically, and purity and integrity were assessed by agarose gel electrophoresis. All samples exhibited 260/280 absorbance ratios of approximately 2.0, and all showed intact ribosomal 28S and 18S RNA bands in an approximate ratio of 2:1 as visualized by ethidium bromide staining.

First-Strand cDNA Synthesis
Total RNA (1 μg to 15 μg) and 2 μL of 50 μM T7-Oligo (dT) Primer were mixed together in a 0.2 mL PCR tube, and RNase-free water was added to create a final volume of 12 μL. The tubes were gently flicked a few times to mix the solution, and then the tubes were briefly centrifuged to collect the reaction products and reagents at the bottom of the tube. The reaction was incubated for 10 minutes at 70°C. The sample was cooled at 4°C for at least 2 minutes. Subsequently, 4 μL of 5× 1st Strand Reaction Mix, 2 μL of 0.1 M DTT, 1 μL of 10 mM dNTP, were combined and mixed thoroughly by flicking the tube a few times. The tubes were centrifuged briefly to collect the reaction products and reagents at the bottom of the tube, and the tubes were immediately placed in a 42°C environment. SuperScript II (1 μL) was added to each RNA sample for a final volume of 20 μL. The RNA was incubated for 1 hour at 42°C; finally, the sample was cooled for at least 2 minutes at 4°C.

Second-Strand cDNA Synthesis
RNase-free water 91 μL, 30 μL of 5× 2nd Strand Reaction Mix, 3 μL of 10 mM dNTP, 1 μL of E. coli DNA ligase, 4 μL of E. coli DNA Polymerase I, and 1 μL of RNase H were mixed by gently flicking the tube a few times. The tubes were centrifuged briefly to collect solids at the bottom of the tube. We added 20 μL of the First-Strand cDNA synthesized as described above. The tube was flicked gently a few times to mix the solution, and then the samples were centrifuged briefly to collect any reaction products and reagents at the bottom of the tube. The cDNA was incubated for 2 hours at 16°C. 2 μL of T4 DNA Polymerase was added to each sample, and the tubes were : Compared with 6 h, the difference was significant (p < 0.05) Compared with normal controls, the expression of most of the 26 genes showed a statistical significant change (p < 0.05). Expressions of PARD3, ALDH8A1, CMTM8 and FBXO3 were also higher than those of normal controls, but there was no statistical significance (p > 0.05)

Table 3: The expression of 26 genes in the very early stage of fulminant hepatitits mice model by real-time RT-PCR (Continued)
incubated for 5 minutes at 16°C. Finally, we added 10 μL of 0.5 M EDTA.

Microarrays
The biotinylated cDNAs were hybridized to 12 individual Affymetrix GeneChips Mouse 430 2.0 Arrays (Affymetrix, Inc., Santa Clara, CA). This product analyzes the expression level for over 39,000 transcripts and variants, including over 34,000 well-substantiated mouse genes.

Data Analysis Significant Differential Gene Analysis
The per-group gene expression value was the geometric mean of RMA normalized gene signals from 3 samples per time point. Due to the fact that the number of samples per class (i.e., three samples per group) was far lower than the number of genes, there were few degrees of freedom for the gene expression signal variance. The low number of degrees of freedom makes an accurate estimation of variability difficult. Significant differential genes among the three groups were filtered by ANOVA and corrected by the random variance model (RVM) [8]. In this experiment, the differential genes from different time points were filtered by RVM ANOVA. To assess the significance of a particular sample having occurred by random chance, the false discovery rate (FDR) was estimated to allow us to determine whether certain genes were actually significant by repeating the comparison test and the permutation test 1000 times. The threshold of truly significant genes was taken to be p-value < 0.05 and FDR < 0.05.

Time Sequence Profile Analysis of Gene Expression
We first select a set of distinct and representative temporal expression profiles independent of the data. These

Construction and Topological Analysis of the Gene Co-Expression Network
In biological processes, macromolecular networks can be constructed by using results from experiments involving Y-2H, coimmunoprecipitation [10], or algorithmic predictions that may be based on gene function correlation and expression profiles [11]. Because of the flexibility of network models based on algorithmic predictions that derive from high-throughput gene expression tests, we can look at snapshots of protein-protein interactions, gene expression regulatory networks and metabolic networks among different groups [12]. Gene networks are constructed from functional gene associations. In the network, cycle nodes represent genes, and edges between two nodes represent interactions between genes. As network elements represent the ways in which genes may regulate other genes, large scale gene networks can be divided into certain subgraphs, named k-core networks, in which all genes are connected to at least k other genes in the subgraph [7]. As a result, the rank of k-core value describes the complexity of the gene association relationship. Cycles with identical colors are part of the same subgraph. In light of the definition of k-core networks, core status within large scale gene networks consists of subgraphs that are associated with higher k values. The complexity of gene relationships increases with k-core value rank. The intrinsic gene network of a given phenotype can represent gene function propriety within a sample. In this work, we use k-cores of protein-protein interaction networks to define major gene functions in the main subgraph. We wanted to identify the main GO on the basis of the maximum number of genes in separate k-cores and then we wished to define the key gene functions at each complexity level of the network [13]. We transformed the normalized expression value of the Pearson correlation into measures of pairwise connection strengths [12]. We selected genes from the two most significant profiles, namely profiles 9 and 16, to construct a co-expression network. The network edges were specified to feature correlation coefficients of above 0.9, to ensure strong gene co-expression relationships. The centrality of the network is represented by the central degree [6]. It is possible to identify characteristic distance variables among genes. We attempted to describe the contri- bution that gene i makes to the status (i.e., the determining central status) of gene j in the network. The maximum core order is termed the main core or the highest k-core of the graph [7]. A k-core subgraph can be generated by recursively deleting the vertices from a graph whose degree is less than k [7]. Larger values clearly correspond to vertices with larger degrees and more central positions in the network structure. In this work, we apply the notion of k-core subgraphing to predict gene function similarity [13]. The nodes labeled by the same color may represent similar gene ontology terms. The highlighted ontologies of same-colored nodes were assessed by counting genes that featured the same gene ontology and the same color.

Real-time RT-PCR confirmation of differentially expressed genes
Real-time RT-PCR was used to verify the differential expression of twenty-six genes that were detected by the Affymetrix GeneChip. The primers used are listed in Additional file 1. Each real-time RT-PCR reaction (in 25 μL) contained 2×SYBR Green Realtime PCR Master Mix, 0.4 μM primers, and 0.5 μL of template cDNA. The cycling conditions consisted of an initial, single cycle of 5 min at 95°C, followed by 40 cycles of 30 sec at 95°C, 30 sec at 54°C, 15 sec at 72°C, and fluorescence acquisition at 83°C for 1 sec. The cDNA was synthesized using reverse transcriptase (Promega), oligo(dT) and random primers with 5 μg RNA from the same samples as those used in the microarray. The PCR amplifications were performed in duplicate for each sample. The gene expression levels were quantified relative to the expression of β-actin by employing an optimized comparative Ct (ΔΔCt) value method. The differences in gene expression levels between groups were compared using the Student's t-test. A p value < 0.05 was considered significant.

Additional material
Authors' contributions FC and ZC conceived of the study and participated in its design and coordination. HHZ and LFZ performed the data analysis, LYZ and JW performed the experiment of animal model, FC and WL drafted the manuscript, and JL and SSW performed RT-PCR. All authors contributed to and approved the final manuscript.