Identification and characterization of dynamically regulated hepatitis-related genes in a concanavalin A-induced liver injury model

Identifying dynamic gene expression changes during concanavalin A-induced acute liver injury using Trendy

We used the Trendy software to analyze global gene expression changes in liver tissues from ConA-induced hepatitis model mice and identify genes with breakpoints/segments (upregulation or downregulation) at 3 h and 24 h after ConA treatment. We identified 115 top dynamic genes with an adjusted R² >0.98 (Supplementary Table 1). We observed breaks or changes in gene expression at 3 h and 24 h after ConA treatment in all the top dynamic genes compared to their corresponding gene expression at the 0 h time point (Figure 1A). This showed at least two time points or segments at which the top dynamic genes showed differential expression in response to ConA treatment. The differential expression patterns (upregulation or downregulation) of all the top dynamic genes by combining both the time points together are shown in Figure 1B.

Figure 1. Dynamics of global gene expression after ConA treatment. (A) The breakpoint distribution of 115 top dynamic genes identified by Trendy with R² > 0.98 is shown. (B) Heatmap shows the expression of top dynamic genes in the liver samples harvested at 0h, 3h, and 24h from 10 mg/Kg ConA treated mice (GSE45413). (C) The most significantly enriched GO terms and KEGG pathways for the top dynamic genes are shown. (D) The GO Chord plot shows the enriched biological function terms for the top dynamic genes. The genes are listed on the left side and their fold change values are shown according to the color scale.

Next, we performed functional enrichment analysis of the top dynamic genes to determine the plausible pathogenetic mechanisms underlying ConA-induced acute liver injury. The most significantly enriched gene ontology (GO) terms related to the top dynamic genes were metabolic process, oxidation-reduction process, response to toxic substances, and positive regulation of apoptotic process (Figure 1C, 1D). The most significantly enriched KEGG pathways were metabolic pathways, chemical carcinogenesis and drug metabolism (Figure 1C). Peng et al. showed that exposure of hepatic cells to toxic substances induced chronic inflammation-related dynamic changes in metabolism [13]. Therefore, we postulate that the top dynamic genes dysregulate metabolic pathways in the ConA-induced hepatitis model mice.

Verification of the top dynamic genes in the ConA-induced hepatitis model mice

We established the ConA hepatitis model mice as shown in Supplementary Figure 1. ConA-treated mice showed significantly enlarged liver, spleen and kidneys and elevated serum ALT and AST levels upon ConA treatment. H&E stained liver sections of ConA-treated mice showed significant infiltration of inflammatory cells, massive hepatocyte necrosis, and disordered hepatic sinusoid structures (Supplementary Figure 1). These findings demonstrated AIH-like characteristics in the ConA hepatitis model mice.

Next, we examined the gene expression of two candidate genes having the extreme breakpoints identified by Trendy, namely, matrix metallopeptidase 3 (Mmp3) and arylacetamide deacetylase (Aadac) in the liver tissues of ConA-induced hepatitis model mice. Trendy analysis showed that the expression of Mmp3 significantly increased at 3 h and peaked at 24 h in the livers of ConA-treated mice (Figure 2A). The mRNA expression of Aadac significantly decreased in the ConA-treated mice livers at 24 h post-ConA treatment (Figure 2A). Both Mmp3 and Aadac play significant roles in the liver functions [14]. Aadac is involved in lipolysis of cellular triacylglycerol stores and the assembly of very low-density lipoprotein (VLDL) [15]. The top 5 dynamic genes (Cd63, Saa3, Slc10a1, Nrxn1, Ugt2a3) differentially expressed in the ConA-treated mice livers were verified by qRT-PCR analysis (Figure 2B). These results confirm that significant changes in the expression of the top dynamic genes correlate with acute hepatitis in the ConA-induced hepatitis model mice.

Figure 2. Validation of gene expression of top dynamic genes in the ConA hepatitis model mice. (A) The expression of two genes with the earliest or latest breakpoint time (Mmp3 at 3 h, and Aadac at 24 h) as determined by Trendy (upper panel) were verified by qRT-PCR analysis (bottom panel) of the liver tissues from the in vivo ConA-treated hepatitis model mice. (B) qRT-PCR results show the mRNA levels of the top 5 dynamic genes (Cd63, Saa3, Slc10a1, Nrxn1, Ugt2a3) in the liver tissues of the ConA-liver injury model mice at 0h, 3h, and 24 h respectively. All data are shown as means ± SEM (n = 5 per group). * p < 0.05; ** p < 0.01; *** p < 0.001.

WCGNA identifies two top modules with hub genes that overlap with top dynamic genes identified by Trendy

Weighted gene co-expression network analysis (WGCNA) was widely utilized to identify the hub genes of diseases. To determine if the top dynamic genes identified by Trendy could also be reproduced by the hub gene analysis, we first constructed gene coexpression networks of ConA-induced liver injury using global transcriptome. We identified the brown module with 910 genes and the turquoise module with 2992 genes as the top 2 modules (Figure 3A).

Figure 3. Identification of top gene modules and hub genes related to liver injury in ConA-treated mice by WGCNA. (A) The cluster dendrogram of 6936 genes in the ConA-treated murine liver samples is shown. The colored bars at the bottom show the color that are designated for specific gene clusters (3 modules). (B) Heatmap shows the correlation between module eigengenes (ME) and the trait (time of liver injury). Each row corresponds to a single module eigengene and the corresponding column represents a trait. Each cell contains the corresponding correlation and p value. (C) Module significance (MS) of each module based on the average absolute gene significance values of all genes in a module are shown for the 3 h and 24 h time points. (D) Venn diagram shows the overlap between the top genes identified by Trendy and the module genes detected by WGCNA. The lower panel shows the significant overlap between the top dynamic genes and the module genes in the turquoise module. (E) Top 10 hub genes in the turquoise module and their gene significance values. (F) Gene co-expression network shows the relationship between Cd63 and its co-expressed genes. The node colors are denoted from green to red (low to high) based on the fold change of gene expression between 3-24 h after ConA administration relative to their expression at 0h.

At 3 h after ConA treatment, the brown module showed the highest module significance (MS) value; the turquoise module showed the highest MS value at 24 h after ConA treatment (Figure 3B–3C). Functional enrichment analyses showed that the brown module was significantly enriched in inflammatory response and cytokine-cytokine receptor interaction pathways, whereas, the turquoise module was significantly enriched in the oxidation-reduction process and chemical carcinogenesis pathway (Supplementary Table 2).

We identified the hub genes in the turquoise module using the network feature selection of WGCNA (Figure 3D). The genes in the turquoise module were ranked according to their gene significance (GS) values. Most of the hub genes in the turquoise module overlapped with the top dynamic genes identified by Trendy (Figure 3E). The highest ranked hub gene in the turquoise module was Cluster determinant 63 (Cd63, GS = 0.996). Cd63 is an exosomal marker in the drug-resistant HCC-derived exosomes [16]. Cd63 and Cd63-related genes were enriched in pathways related to chemical carcinogenesis and metabolic pathways (Figure 3F).

Cd63 silencing ameliorates ConA-induced hepatic injury in mice

We then investigated the effects of silencing Cd63 in the ConA-induced hepatitis model mice. The mice injected with lentiviruses carrying Cd63-specific shRNAs showed 50% reduction in the Cd63 mRNA levels in the liver tissues compared to the corresponding controls (Figure 4A). The survival rates of the Cd63-silenced ConA-induced mice were significantly higher than the corresponding control group mice (Figure 4B). We also analyzed serum ALT and AST levels in the control and Cd63-silenced groups of mice at 8 h after ConA injection. The serum ALT and AST levels, which were estimated at 8 h after ConA injection, were significantly lower in the Cd63-silenced group mice compared to the control group mice (Figure 4C).

Figure 4. Cd63 silencing protects against in vivo ConA-induced liver injury. (A) qRT-PCR analysis shows Cd63 mRNA expression in the liver tissues of control (sh-Ctrl) and Cd63 knockdown (sh-Cd63) group mice. (B) Survival curves show overall survival rates of sh-Ctrl and sh-Cd63 group mice after 15mg/Kg ConA treatment. (C) Comparison of serum ALT and AST levels in the sh-Ctrl and sh-Cd63 group mice treated with ConA. (D) qRT-PCR analysis shows the relative mRNA levels of pro-inflammatory genes, IL-1β, IFN-γ, IL-2, IL-6, and TNF-α in the liver tissues of sh-Ctrl and sh-Cd63 group mice. All data are shown as means ± SEM (n = 5 per group). * p < 0.05; ** p < 0.01; *** p < 0.001. (E) Representative images show H&E stained liver sections of sh-Ctrl and sh-Cd63 group mice. (F) Representative images show TUNEL stained liver sections of sh-Ctrl and sh-Cd63 mice.

H&E stained liver sections estimated at 24h showed widespread tissue necrosis in ConA-treated control group mice, but necrosis was significantly reduced in the Cd63-silenced group mice treated with ConA (Figure 4E). TUNEL assay results showed that hepatocyte death was significantly reduced in the Cd63-silenced group mice treated with ConA compared to the ConA-treated control group mice (Figure 4F). qRT-PCR analysis showed that the mRNA levels of pro-inflammatory mediators such as IL-1β, IFN-γ, IL-2, IL-6, TNF-α were significantly lower in the liver tissues from the Cd63-silenced group mice compared to the liver tissues from the ConA-treated control group mice (Figure 4D). Overall, our results suggest that Cd63 deletion significantly reduces liver injury in the ConA-treated mice.

Transcriptome data analysis

The time-course transcriptome profile of the ConA hepatitis murine model was retrieved from the GSE45413 dataset in the GEO database [26]. This dataset was generated with liver tissues collected at 0 h, 3 h, and 24 h after 10mg/kg ConA injection into 8-12 week old C57BL/6 male mice. The microarray data was normalized by selecting probes for genes with a mean expression in the top 75% and removing the lowly-expressing genes. Then, we selected the genes with above average expression variance. We then collapsed different probes that targeted the same gene, resulting in 6936 genes [27]. Finally, we identified differentially expressed genes using the limma R package.

Identification of top dynamic genes using Trendy

We used the Trendy R package to characterize dynamic gene-specific expression patterns over a time-course during acute liver injury [28]. In brief, Trendy fits a set of segmented regression models with varying numbers of breakpoints for each gene. Each breakpoint represents a dynamic change of gene expression over time. Genes with high R² values were categorized as top dynamic genes. Then, the parameter estimates of the optimal model including the sign and p-value of the slope estimate were used to determine the direction (up, down, or no-change) of the changes in the expression of the top dynamic genes over time.

Weighted gene co-expression network analysis

The top dynamic genes were evaluated by weighted gene co-expression network analysis (WGCNA) to identify gene modules and hub genes [29]. In brief, the eigengene module was identified based on the weighted average of the gene expression profiles by evaluating the matrix of pair-wise Pearson's correlation coefficients. Then, the gene significance (GS) was computed for each gene within the eigengene module at all time points after ConA injection (3 h and 24 h). The geometric mean was then calculated for the absolute values of all the GS values within each module to determine the module significance (MS) of each module. Modules with higher MS values significantly correlated with the trait (time of liver injury). The network of the module genes was visualized using Cytoscape.

Functional enrichment analysis

We used the DAVID database (https://david-d.ncifcrf.gov/) [30], and GOplot [31] to determine significant gene ontology (GO) terms and Kyoto Encyclopedia of Gene and Genomes (KEGG) pathways related to the module genes.

Establishment of Cd63 knockdown mice

We obtained 8-12 week old C57BL/6 male mice weighing 20-25g from the Shanghai Slac Laboratory Animal Co. Ltd (Shanghai, China). The mice were housed in a specific-pathogen-free facility with a consistent room temperature and humidity. We generated lentiviruses carrying shRNA-Cd63 (CCAGGTGAAGTCAGAGTTTAA) or control scrambled shRNA (shRNA-Ctrl) vector as previously described [32]. Four weeks before the ConA injection, lentiviruses carrying the sh-Cd63 or sh-Control (sh-Ctrl) was injected into the tail vein of mice (n=5/group).

ConA-induced hepatitis model mice

We injected mice with 10 mg/kg ConA (prepared in saline) through the lateral tail vein. For the survival assay, 15 mg/kg ConA was used. At the indicated time points (0 h, 3 h, and 24 h), blood samples were obtained through retro-orbital bleeds and serum samples were prepared and stored at -80° C until measured by automatic biochemical analyzer (Hitachi Auto Analyzer 7170, Japan) for aspartate aminotransferase (AST) and alanine aminotransferase (ALT). Then, the mice were sacrificed and liver samples were harvested and stored for further experiments.

Quantitative real time PCR

Total RNA was extracted from the murine liver tissues using Trizol (Invitrogen). Equal amounts of RNA were reverse-transcribed into cDNA using the PrimeScript RT reagent kit (Takara Bio). Then, gene expression was analyzed by qPCR with the SYBR Premix Ex Taq kit (Takara Bio) in the ABI 7900 Real-Time PCR System. The expression of various genes relative to GAPDH (housekeeping gene) was determined using the 2^-ΔΔCt method.

Histopathology assays

The murine liver tissue samples were fixed in 4% paraformaldehyde for 48 h and then paraffin embedded. Then, 4-5 μm thick paraffin embedded liver sections were cut, deparaffinized with xylene, rehydrated with decreasing concentrations of ethanol, and stained with hematoxylin and eosin (H&E). Cryosections of liver tissues were fixed with 4% paraformaldehyde in PBS for 15 min. Then, relative levels of apoptosis were quantified in all samples using the TUNEL assay with the In Situ Cell Death Detection Kit (Roche) according to the manufacturer's instructions.

Statistical analysis

The data are expressed as means ± SEM. The statistical differences between samples were compared using 2-tailed Student’s t tests. p<0.05 was considered statistically significant.

Ethical standards

The experiments were carried out according to the protocols approved by the Ethical Committee of the Affiliated Hospital of Hangzhou Normal University (Approval No. 2019(02)-HS-51).