The GSE profiles (GSE84044) based on GPL570 (Affymetrix Human Genome U133 Plus 2.0 Array) was selected and downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geoprofiles/). GSE84044 dataset was established and evaluated by Wang M et al. (11). Unlike their focus on fibrosis-related genes, the present study would focus on the immune mechanisms in the development of HF. GSE84044 dataset contained a gene expression matrix and clinical information of 124 patients with HF. Clinical information for each HF patient included gender, age, grading of inflammation (Scheuer score ‘G’) and Histological stage of fibrosis (Scheuer score ‘S’). Scheuer score was used to evaluate liver inflammation and fibrosis. G0-G1, G2 and G3-4 were defined as no or mild, moderate and severe liver inflammation respectively; S0-S1, S2-S3 and S4 indicated no significant hepatic fibrosis, moderate hepatic fibrosis, severe hepatic fibrosis or cirrhosis, respectively. After standardization of GSE84044 gene expression matrix, the gene expression matrix and clinical information of GSE84044 were used for the construction of a gene co-expression network.

The Weighted Gene Co-Expression Network Analysis (WGCNA) of GSE84044 gene expression matrix was performed using R software (version: 3.6.3). Briefly, the top 50% of genes with the smallest MAD (Median Absolute Deviation) were eliminated and the outlier genes and samples were also removed. A weighted adjacency matrix was constructed. Then, the appropriate soft-thresholding value β was chosen and a topological overlap matrix (TOM) was transformed from the adjacency matrix. The average linkage hierarchical clustering and dynamic tree cut algorithm were used to identify gene co-expression modules. To obtain modules related to clinical features, module eigengenes (MEs), as the first principal component of the module, were used to compute the correlation with clinical traits and the modules with pvalue <  0.001 were identified to be clinical traits associated modules. In the study, Common significant modules related to fibrosis and inflammation were screened. Finally, gene significance (GS) and module membership (MM) were calculated, the module genes with GS > 0.4 and MM > 0.5 and pvalue < 0.05 were identified as crucial module genes and the common crucial module genes related to fibrosis and inflammation were identified for subsequent analysis.

2483 immune-related genes (IRGs) were obtained from the ImmPort database (https://www.immport.org/shared/home). The overlapping genes between IRGs and crucial module genes related to fibrosis and inflammation were screened out and identified as hepatic fibrosis immune genes (HFIGs).

DAVID Database (https://david.ncifcrf.gov/) is a reliable program involved in functions like functional annotation, gene functional classification, gene ID conversion. In this study, DAVID Database was used to performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of HFIGs. The enrichment terms with pvalue < 0.05 were considered as significant enrichment results.

The 98 HFIGs were uploaded to the String database (https://string-db.org/) for protein-protein interaction (PPI) analysis. PPI networks were constructed using Cytoscape software. Then, 10 hub genes were identified from the PPI network using the CytoHubba plugin based on the Maximal Clique Centrality (MCC) algorithm. To further explore the biological functions of hub genes, 10 hub genes were uploaded to the Metascape database (http://metascape.org/gp/index.html#/main/step1) for functional enrichment analysis.

The role of the immune microenvironment in the formation of HF was analyzed using the CIBERSORTx database (https://cibersortx.stanford.edu/). Gene expression data from each patient in the GSE84044 dataset were uploaded to the CIBERSORTx database to calculate the level of infiltration of 22 immune cells, and then correlations between immune cells and clinical trait score (fibrosis stage and inflammatory grade) and between hub genes and immune cells were analyzed.

Male Sprague-Dawley rats (180–220 g) were provided by the Laboratory Animal Center of Ningxia Medical University. The animal experiment was approved by the Ethical Committee of Ningxia Medical University (IACUC-NYLAC-2021-009). Rats were maintained on a 12 h light/dark cycle under a favorable environment (22–24 °C, 45–50% humidity) and had free access to water and rodent chow.

After 1 week of acclimatization, the animals were randomly divided into the normal control and CCl4-induced (Bodi Chemical, Tianjin, China) model group of 6 animals each. The rats in CCl4-induced model group received olive oil with CCl4 subcutaneously at a dose of 23 mL/kg twice a week for 6 weeks to induce hepatic fibrosis, and olive oil was given to normal control rats, and the rats in the normal control simultaneously were received the same volumes of olive oil. 24 hours after the last subcutaneous injection, all rats were anesthetized with 3% pentobarbital sodium (Bodi Chemical, Tianjin, China). The liver tissues were taken for further examination including histological analysis, ELISA assay, transcriptome sequencing and q-PCR assay.

The study protocol involving the use of experimental rats was reviewed and approved by the Research Ethics Committee of Ningxia Medical University and was in accordance with the Animal Management Rules of the Ministry of Health of the People’s Republic of China.

The liver tissues were fixed in 10% formalin for 24 h, dehydrated, embedded in paraffin and then cut into 5 mm tissue sections. Liver sections were stained with hematoxylin and eosin (HE) (Solarbio, Beijing, China) staining for histological examination and stained with Masson’s Trichrome (Sangon, Shanghai, China) staining to assess the degree of hepatic fibrosis.

The liver collagenic parameters markers Hyaluronic acid (HA) and Laminin (LN) were measured by ELISA kit (Bioswamp, Wuhan, China) according to the manufacturer’s protocol. The levels of interleukin-1β (IL-1β), interleukin-6 (IL-6) and tumor necrosis factor α (TNF-α) in liver tissues were estimated by using ELISA kit (Lc-Bio Technologies, Hangzhou, China) according to the manufacturer’s instructions.

Here, immunofluorescence staining was performed to detect the infiltration of each immune cell. The expression of Cd10 (Abcam, Cambridge, UK), Cd11b (Abcam, Cambridge, UK), Cd86 (Abcam, Cambridge, UK) and Cd206 (Abcam, Cambridge, UK) were detected to assess the infiltration status of neutrophils, natural killer (NK) cells, macrophages M1 and macrophages M2, respectively. Briefly, After dewaxing to water and antigen repair, the paraffin-embedded liver tissue sections were incubated in 3% hydrogen peroxide for 10 minutes to quench endogenous peroxidase activity. The sections were then permeabilized with 0.2% Triton X-100 (Abcam, Cambridge, UK) for 15 min and blocked with 5% BSA (Abcam, Cambridge, UK) for 30 min. Tissue slices were incubated primary antibodies at 4°C overnight and then were treated with secondary antibodies at room temperature for 30 min. Nuclei were counterstained with 4,6-diamidino-2-phenylindole (DAPI) (Thermo Fisher Scientific, Wilmington, DE) and the slices were mounted in Pro-Long Diamond Antifade Mountant (Thermo Fisher Scientific, Wilmington, DE). Sections were observed and captured with a fluorescence microscope (Nikon Corporation, Tokyoratio.

Total RNA was extracted from liver tissue of 6 rats (3 each for the Model and control groups) using TRlzol Reagent (Life technologies, California, USA) according to the manufacturer’s instructions. The integrity of total RNA was determined using Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA) and the concentration and purity of total RNA were detected by NanoDrop 2000 (Thermo Fisher Scientific, Wilmington, DE). The cDNA library was generated using NEBNextUltra™ RNA library prep kit for Illumina (NEB, Ipswich, MA, USA) following manufacturer’s recommendations and the index codes were added to attribute sequences to each sample. The library quality was developed on Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). After cluster generation which was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v4-cBot-HS (Illumia) according to the manufacturer’s instructions, the library were sequenced on an Illumina platform.

Clean reads were obtained by removing reads containing adaptor sequences, adaptor sequences and reads containing ploy-N from raw data. Quantification of gene expression levels was then calculated using the reads per kilo bases per million reads (RPKM) method. | Fold change |≥1.5 and adjusted pvalue<0.01 were defined as the threshold for screening differentially expressed genes (DEGs) using the DESeq2. The expression of 10 hub genes was compared in normal control and model groups, and the expression correlations between hub genes were calculated.

To verify the reliability of the functional enrichment analysis results of HFIGs, the DEGs of normal rats and hepatic fibrosis rats in RNA-seq data were used for KEGG function enrichment analysis by DAVID database and GSEA analysis respectively. At first, DEGs were uploaded to the DAVID database for KEGG pathway enrichment analysis. pvalues were calculated using the Benjamini-corrected modified Fisher’s exact test and pvalue<0.05 was defined as a threshold of significance. Then, Gene Set Enrichment Analysis (GSEA) was performed on the transcriptome expression matrix using GSEA v4.0.3 software, the KEGG pathways with |NES|>1 (Normalized enrichment score) and pvalue < 0.05 were considered to be significantly enriched.

Total RNA of each sample was extracted using Trizol. The total RNA was reverse‐transcribed into cDNA using the PrimeScript RT Master Mix kit (Takara Bio, Kusatsu, Japan). The mRNA expression of hub genes was quantitated by q-PCR using SYBR Green PCR Mix (Monad Biotech, Wuhan, China) on an CFX Connect machine (Applied BIO-RAD, Hercules, CA, USA). β-actin was used as internal control. Primers used in q‐PCR are listed in Table 1 .

Connectivity Map (CMAP) database (https://clue.io) is a reliable platform for predicting potential small molecules that may reverse or induce the expression of genes encoded in specific biological states. The 10 hub genes upregulated in liver fibrosis were introduced into the CMAP database to screen for small molecules with potential therapeutic effects in HF. The results were ranked by connectivity score. Small molecules with negative connectivity scores indicating the ability to reverse the upregulation of hub gene expression in HF, were considered to have anti-hepatic fibrosis potential.

The “ggplot2” software package is used for image generation in bioinformatics analysis. GraphPad Prism 8.3.0 software (San Diego, CA, USA) was used for statistical analysis and image construction. Pearson correlation was used to evaluate the linear relationship between data. The experimental data were presented as mean ± standard deviation (SD). Comparisons between groups was evaluated using the Student’s t-test. *pvalue< 0.05, **pvalue< 0.01, ***pvalue< 0.001 were considered significant.

After standardization, The WGCNA analysis of GSE84044 data was performed. To ensure the network was scale-free and more biological significance, the optimal soft threshold β was set as 5 ( Figures 1A, B ). Then, based on the DynamicTreeCut algorithm, setting the minimum number of module gene as 30, the sensitivity dynamic tree cut (deep split) as 3, and the maximum module distance as 0.25, the gene modules were generated and the modules with high similarity were further merged. As showed in Figure 1C , 21 gene modules were finally generated after merging. Eigengene dendrogram and eigengene adjacency were ploted to analyze the connectivity of module eigengenes (MEs), and the results showed that the distance between modules was greater than 0.25 ( Figure 1D ).

The module-trait correlation coefficients were calculated and the results were showed in Figure 1E . The results showed that 5 modules including brown4, black, floralwhite, darkturquoise and red were highly relevant with fibrosis and inflammation (pvalue<0.001). Hence, the 5 modules significantly associated with fibrosis and inflammation were selected for further analysis. According to the criteria of |GS|>0.4 and |MM|>0.5 and pvalue<0.05, a total of 684 crucial genes which highly associated with both fibrosis and inflammation were identified from 5 significant modules ( Figures 2 , 3 ).

At last, comparing the 2483 immune genes obtained from the ImmPort database, a total of 98 overlapping genes between IRGs and crucial module genes related to fibrosis and inflammation were screened out and identified as hepatic fibrosis immune genes (HFIGs) ( Figure 4 ).

In order to determine the biological functions of HFIGs, 98 HFIGs were uploaded to DAVID Database to perform functional enrichment analysis including GO annotation and KEGG pathway enrichment analysis. GO annotation result showed that the 98 HFIGs were significantly enriched in 149 biological processes (BP), 32 cellular components (CC), and 39 molecular functions (MF), (pvalue<0.05). The BP involved in immune response, chemotaxis, neutrophil chemotaxis, inflammatory response, signal transduction, chemokine-mediated signaling pathway, adaptive immune response, antigen processing and presentationn, killing of cells of other organism, etc. The CC terms included external side of plasma membrane, extracellular space, cell surface, MHC class II protein complex, extracellular region, plasma membrane, ER to Golgi transport vesicle membrane, integral component of plasma membrane and so on. In terms of MF, 98 HFIGs were mostly enriched in chemokine activity, peptide antigen binding, CCR chemokine receptor binding, MHC class II receptor activity, cytokine activity, CXCR chemokine receptor binding, MHC class II protein complex binding, T cell receptor binding and so on. The top 10 most signifcantly enriched GO terms in BP, CC and MF were displayed in Figure 5A . In addition, a total of 24 KEGG pathways were obtained after removing the disease related signal pathways ( Figure 5B ). The obtained KEGG pathways involved in cytokine-cytokine receptor interaction, viral protein interaction with cytokine and cytokine receptor, Th17 cell differentiation, antigen processing and presentation, chemokine signaling pathway, TNF signaling pathway, NF-kappa B signaling pathway, Toll-like receptor signaling pathway and so on. Most of the enriched to pathways were immune-related pathways ( Figure 5C ). The connectivity relationships of the HFIGs and their enriched KEGG pathways were displayed in circle chart ( Figure 5D ).

Using String database and Cytoscape software, the PPI network was constructed ( Figure 6A ). According to MCC algorithm, the 10 top hub gene were screened out from the PPI network and sequentially ordered as follows: CXCL8, IL18, CXCL10, CD8A, IL7, PTPRC, CCL5, IL7R, CXCL9 and CCL2 ( Figure 6B ). The expression of these 10 bub genes was positively correlated with fibrosis stage (S0-S4) and inflammatory grade (G0-G4) ( Figures 6C, D ). All of the hub gene expressions were significantly correlated with each other (pvalue<0.001). High correlations between CXCL9 and CXCL10 (Correlation coefficient=0.96) and between CCL5 and CD8A (Correlation coefficient=0.95) were observed ( Figure 6E ).

Function enrichment analysis showed that the hub genes were enriched in a variety of KEGG pathways. Compared with KEGG pathway enrichment results of 98 HFIGs, 8 common signaling pathways were identified, including cytokine-cytokine receptor interaction, viral protein interaction with cytokine and cytokine receptor, chemokine signaling pathway, Toll-like receptor signaling pathway, NOD-like receptor signaling pathway, IL-17 signaling pathway, hematopoietic cell lineage and TNF signaling pathway ( Figure 6F ).

Firstly, we analyzed the infiltration of 22 immune cell types in the liver tissue of 129 patients with hepatic fibrosis. The bar graph and heat map showed the type and number of immune cell infiltrates in each sample ( Figures 7A, B ). 21 infiltrating immune cell types were identified in liver tissue, except for eosinophils, which were not expressed in each sample. Subsequently, we assessed the correlation between these immune cell populations ( Figure 7C ). A positive correlation between macrophages M1 and T cells gamma delta (r=0.56) and a negative correlation between macrophages M2 and macrophages M1 (r=-0.72) and T cells gamma delta (r=-0.58) was found.

In addition, the correlation between the infiltration level of 21 immune cells in the liver tissue of each patient and fibrosis stage (S0-S4) and inflammatory grade (G0-G4) was assessed ( Figures 8A, B ). Fibrosis stage (S0-S4) showed a positive correlation with the infiltration status of macrophages M1 (pvalue<0.001), T cells gamma delta (pvalue<0.001) and B cells naive (pvalue<0.001), and a negative link to the infiltration status of B cells memory (pvalue<0.001), macrophages M2 (pvalue<0.001), NK cells resting (pvalue<0.001), neutrophils (pvalue<0.01), dendritic cells activated (pvalue<0.05)and T cells regulatory (Tregs) (pvalue<0.05). Inflammatory grade (G0-G4) was positively correlated with T cells gamma delta (pvalue<0.001), macrophages M1 (pvalue<0.001) and B cells naive (pvalue<0.05), and negatively correlated with B cells memory (pvalue<0.001), macrophages M2 (pvalue<0.001), NK cells resting (pvalue<0.001), mast cells resting (pvalue<0.001), T cells regulatory (Tregs) (pvalue<0.01) and neutrophils (pvalue<0.01).

Furthermore, We further investigated the correlation of 10 hub genes with 21 immune cell infiltration ( Figures 9A–J ). Each hub gene expression was positively correlated with the infiltration status of macrophages M1, T cells gamma delta and B cells naive, and negatively correlated with the infiltration status of macrophages M2, B cells memory, T cells regulatory (Tregs) and NK cells resting.

To verify the reliability of the results of WGCNA analysis, we constructed a rat model of HF. The results of HE and Masson staining were shown in Figures 10A, B . The liver tissues in the control group showed a normal structure. However, The liver tissue of model group rats exhibited a large amount of hepatocyte steatosis, pseudolobule formation, hepatocyte swelling and necrosis, infiltration of immune cells and extensive distribution of collagen fibers.

Fibrosis and inflammation parameters of liver tissue in the normal control and model groups were measured using ELISA kits. The results were shown in Figures 11A–E , compared with the normal control group, The liver collagenic parameters Hyaluronic acid (HA) (pvalue<0.001) and Laminin (LN) (pvalue<0.001) and the liver inflammatory cytokines IL1β (pvalue<0.001), IL6 (pvalue<0.001) and TNF-α (pvalue<0.01) were significantly increased in the CCl4 model group.

Then, we performed immunofluorescence analysis on the liver tissue of the normal group and the model group. The results were shown in Figure 12 . The expression level of Cd86 (a marker of macrophages M1) were significantly higher in fibrotic rat liver than that in normal rat liver tissues, but Cd10 (a marker of neutrophils), Cd11b (a marker of NK cells), and Cd206 (a marker of macrophages M2) expression level was significantly decreased in the CCl₄ model group compared with the normal control group. The results of immunofluorescence staining verified the accuracy of the results of immune infiltration analysis to some extent.

The expression of Cd10, Cd11b, Cd86 and Cd206 were detected to assess the infiltration status of neutrophils, NK cells, macrophages M1 and macrophages M2, respectively.

Liver tissues from normal rats and CCl4-induced HF rats were examined using high-throughput sequencing technology, and then screened for differential genes. Setting |Fold change|>1.5 and pvalue<0.05 as the threshold, 1039 annotated differential genes were screened out, including 744 up-regulated genes and 295 down-regulated genes ( Figure 13A ). Using DAVID database, the 1039 DEGs were enriched in 14 KEGG pathways that were identified as the HFIGs related KEGG signaling pathways in the previous results, including osteoclast differentiation, chemokine signaling pathway, natural killer cell mediated cytotoxicity, complement and coagulation cascades, Toll-like receptor signaling pathway, PI3K-Akt signaling pathway, NOD-like receptor signaling pathway, Apoptosis, NF-kappa B signaling pathway and so on ( Figure 13B ).

Furthermore, we investigated the signaling pathways that were significantly differentially expressed in the model and normal control groups by using GSEA. Compared with 24 pathways that were HFIGs-enriched in the previous results, 19 common KEGG pathways were identified. As shown in Table 2 , all the common pathways were significantly overexpressed in the model group compared to the normal control group.

9 hub genes were screened from RNA-seq data, including Il18, Cxcl10, Cd8a, IL7, Ptprc, Ccl5, Il7r, Cxcl9 and Ccl2. Cxcl8 did not have any annotation information in rat species and may not be expressed in rat species. The 9 hub genes were significantly upregulated in liver tissue of the CCl4-induced hepatic fibrosis rats compared to the normal liver tissue ( Figure 13C ). Person correlation analysis showed good correlation between all the 9 hub genes ( Figure 13D ). Overall, the results of our transcriptomic analysis further confirm the key roles of these immune-related signaling pathways and hub genes in the progression of HF.

Since Cxcl8 is not expressed in rat species, the mRNA expression levels of nine hub genes including Il18, Cxcl10, Cd8a, IL7, Ptprc, Ccl5, Il7r, Cxcl9 and Ccl2 were detected by q-PCR experiments to verify the expression of identified hub genes in liver tissues of the normal group and CCl4 model group. The result was showed in Figure 14 . The mRNA expression of these 9 hub genes in the CCl4 model group were significantly up-regulated compared to the normal control group, consistent with previous WGCNA analysis and RNA sequencing results (pvalue<0.05, pvalue<0.01, pvalue<0.001).

10 hub genes were uploaded to CMAP to determine potential drugs for reversing HF. The top 15 small molecules with negative connectivity scores were screened out, including telotristat, lomitapide, malotilate, tranylcypromine, TG-101348, cordycepin, tafamidis-meglumine, fostamatinib, delcorine, phenylbutazone, endo-IWR-1, scopolamine, anidulafungin, simvastatin and terconazole ( Table 3 ). These small molecule compounds have the potential to exert anti-hepatic fibrosis effects by inhibiting the upregulation of hub genes.