Identification and validation of immune related core transcription factors GTF2I in NAFLD

Patients

Liver tissues used in the study were obtained from patients with hepatic hemangioma who underwent hepatic hemangioma resection in Tongji Hospital (Shanghai, China). The use of specimens was approved by the Research Ethics Committee of Tongji Hospital, and the Institutional Review Board approval number is K-KYSB-2020-139. All patients signed an informed consent to participate in this research.

Cell culture and hepatic steatosis model constructed

At the condition of 37 °C, 5%CO2, HepG2 (HB-8065, ATCC) was cultured in 6-well plates (2*10⁵ cells/well) of DMEM Medium (KGM12800-500, Keygene BioTECH, Nanjing, China) supplemented with 10% FBS (Gibco, South America) and 1% penicillin-streptomycin solutions. Cells were starved for 6 h in serum-free DMEM, 0.5% BSA (Sigma-Aldrich, St. Louis, MO, USA). Bovine albumin (BSA, 1mM) and free fatty acids (FFA,1mM) were involved in the medium for 24 h to induce control/steatosis cells (Kanuri & Bergheim, 2013; Chavez-Tapia, Rosso & Tiribelli, 2011).

Data acquisition

The inclusion and exclusion criteria of participants were as described in previous studies. Specifically, microarray data and clinical information of 26 liver samples from healthy participants and 109 liver samples from NAFLD patients were achieved from the Gene Expression Omnibus (GEO) database (accession number: GSE126848, GSE130970) (Brenner et al., 2013; Sun & Gao, 2004). The participants included in this analysis were all diagnosed with histologically normal and NAFLD liver by biopsy, and interference from excessive alcohol consumption, diabetes mellitus, drug injury and other factors were excluded. In addition, the annotation files were acquired from GPL18573 Illumina NextSeq 500 (Homo sapiens) and GPL16791 Illumina HiSeq 2500 (Homo sapiens) platforms (Hoang et al., 2019; Suppli et al., 2019). Hallmark signaling pathways were collected from the Molecular Signatures Database. 29 immune-associated gene sets which represented diverse immune cell types, functions, and pathways were used to identify immune signatures (He et al., 2018).

Data prerecession and differential expression analysis

Participants with incomplete samples and clinical information were excluded. Affy package was utilized for reading original microarray data. Robust multi-array average (RMA) background was to correct, standardize, probe-specific background correction, and summarize probe set values in expression measure. Meanwhile, normalization was used to correct the integrated microarray data in preparation for differential expression analysis. Differential expression analysis of genes between NAFLD and normal liver were conducted using the Linear Models for Microarray Data (limma) package. False discovery rate (FDR) was utilized for multiple testing. FDR < 0.05 as well as absolute log₂[Fold Change (FC)] ≥ 1.0 were cut off criteria.

Functional enrichment analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes & Genomes (KEGG) enrichment analyses were carried on to investigate intrinsic biological processes and pathways which were significantly related to DEGs with FDR p-value < 0.05 as cut off value. Furthermore, we obtained all transcription factors (TFs) from the Cistrome Cancer database (http://cistrome.org/CistromeCancer/), which were then merged with DEGs to acquire the key differentially expressed transcription factors (DETFs).

Identification of potential downstream hallmark pathways

Potential downstream hallmark pathways of DETGs were identified through GSVA and GSEA (Ferreira et al., 2021; Mei et al., 2017). 50 hallmark pathways (GSEA | MSigDB (gsea-msigdb.org)) were absolutely quantified and then evaluated to extract the hallmark differential expression pathways related to NAFLD through ClusterProfiler package and GSVA package (Subramanian et al., 2007; Yu et al., 2012). Significant enrichment of upregulated and downregulated hallmark pathways in NAFLD and normal livers was indication through GSEA.

Regulation network of transcription factors

Key DETFs ultimately associated with NAFLD were obtained through combing the Pearson correlation analysis of DETFs (correlation more than 0.85; significance less than 0.05), merging with the results of GSVA and GSEA analysis, and then plotted the regulatory network combining TFs, DEGs, and hallmark pathways using the igraph R package (https://igraph.org/r/).

Oil red O staining

After HepG2 was treated with FFA for 24 h, the medium was removed. Wash twice with phosphate buffered saline (PBS) and fix with 4% paraformaldehyde for 30 min. Then, add oil red O (0.6% oil red O in isopropanol: H₂O = 3:2) for 1 h and washed with PBS three times. The formation of lipid droplets was observed under optical microscope.

Hematoxylin-eosin staining (HE) and immunohistochemistry (IHC)

Hematoxylin-eosin Staining (HE) was performed on liver tissue (5 mm) of three healthy subjects and 10 patients with NAFLD. Following routine rehydration, antigen retrieval, and blocking procedures, the sections were incubated overnight with GTF2I antibody (1:50 dilution, ab248269, Abcam, Cambridge, UK) at 4 °C. All slides were labeled polymer HRP for half an hour and hematoxylin as a counterstain for 5 min at room temperature.

Western blot

30 mg of each sample was taken and lysed with RIPA at 4 °C to obtain the total protein. After protein quantification using a microplate reader, an appropriate amount of loading buffer was added and subsequently changed by heating in 95 °C for 15 min. After polyacrylamide gel electrophoresis (30 g protein usage) and subsequent routine western blot experimental steps, GAPDH (diluted 1:25,000, 60004-1, Proteintech, USA), CHD2 (diluted 1:1,000, 12311-1-AP, Proteintech, USA) were respectively incubated overnight with the corresponding PVDF membrane region with the GTF2I (1:1,000, ab248269, abcam, USA) antibody at 4 °C overnight. Then, the membranes were incubated with the horseradish peroxidase-conjugated goat anti-rabbit IgG secondary antibody (cat. no. A0208; 1:10,000; Beyotime Institute of Biotechnology, Jiangsu, China) or goat anti-mouse IgG (cat. no. A0216; 1:10,000; Beyotime Institute of Biotechnology) at room temperature for 1 h. BeyoECL Plus (Beyotime Institute of Biotechnology, Jiangsu, China) and QuantityOne v4.6.6 software (Bio-Rad Laboratories, Inc., Hercules, CA, USA) were used to observe the protein bands.

Quantitative reverse transcriptase-PCR (qRT-PCR)

Total RNA of liver specimens (30 mg) was extracted via Trizol reagent as described by the manufacturers, and RT-PCR were performed with total RNA (1 µg) according to the instructions of the RT-PCR reaction kit. The product was amplified in a reaction volume of 20 µl containing 1 µl RT product, 10 µl SYBR H (5X), and 0.5 µl of each primer (10X). PCRs were reacted for 40 cycles at 95 °C for 10 s, 60 °C for 30 s, 72 °C for 30 s. Samples’ threshold cycle values were standardized to GAPDH mRNA expression, while the fold-change for each mRNA was calculated using the 2-delta delta Ct method. Primers of qRT-PCR were listed in Table 1.

Immune profile assessment

According to the median of GTF2I in all samples, the samples were divided into high expression group and low expression group. CIBERSORT algorithm was used to analyze the effect of GTF2I on different immune cell composition. GSEA was used to analyze the differential pathways between high expression group and low expression group. Immune cell surface markers were used to further compare the difference of high expression group and low expression group.

Statistics analysis

Statistical analyses were performed with R version 3.6.1 (R Core Team, 2019). The continuous variables were presented as the mean ± standard deviation. For non-microarray data, data that satisfied normal distribution were analyzed by Student’s t-test. A Mann–Whitney U test was used to evaluate non-parametric data. P value < 0.05 was supposed to statistically significant.

Identification of DEGs and functional enrichment analysis

The analysis procedure of this study is shown in Fig. 1. The participants included in this analysis were all diagnosed with histologically normal and NAFLD liver by biopsy, and interference from excessive alcohol consumption, diabetes mellitus, drug injury and other factors were excluded. On this basis, we integrated the expression matrix of the microarray data for differential gene analysis. As shown in the volcano plot (Fig. 2A) and heat map (Fig. 2B), 617 genes were identified as DEGs from 18435 genes in 26 controls and 109 NAFLD patients.

Next, GO and KEGG enrichment analyses were used to explore the potential function of identified DEGs. In GO analysis, the most significant terms of biological process (BP), cellular components (CC), and molecular function (MF) were ‘establishment of protein localization to membrane’ (Fig. 2C), ‘cytosolic ribosome’ (Fig. 2D), and ‘structural constituent of ribosome’ (Fig. 2E), respectively, which were critical in NAFLD pathogenesis. In KEGG analysis, DEGs between NAFLD and normal participants were involved in cell signal transduction and immune related diseases, such as ‘Glycine, serine and threonine metabolism’, ‘Cell adhesion molecules’, ‘autoimmune thyroid disease’ and ‘Allograft rejection’ (Fig. 2F). Through the functional analysis of differential genes between NAFLD and normal participants, these results suggested that the expression pattern of differential genes in NAFLD may be closely related to the immune profile, which we will investigate in the next step.

Identification of key differentially expressed transcription factors (DETFs)

The results of the analysis above showed that differential genes were enriched in immunity and signal transduction. We further screened transcription factors on the basis of DEGs to obtain key regulatory genes.

All TFs were obtained through the Cistrome Cancer database, and then crossed with the DEGs to get the relevant transcription factors that are differentially expressed. As the volcano plot (Fig. 3A) and heat map (Fig. 3B) illustrated, six TFs (ESR1, CHD2, GTF2I, EGR1, HCFC1, SP2) were considered as the DETFs in patients with NAFLD. Among these transcription factors, CHD2, ESR1, and GTF2I were significantly increased in NAFLD. EGR1, HCFC1 and SP2 were significantly decreased in NAFLD. These transcription factors would prepare for the later construction of co-expression networks among immune signatures, pathway hallmarks, DEGs and DETFs.

Immune signatures and pathway hallmarks between NAFLD and normal participants

Based on the enrichment analysis results, we further deeply evaluated the distinction of immune profile in NAFLD and normal participants. 29 immune-associated gene sets which represented diverse immune cell types, functions, and pathways were used to identify immune signatures between NAFLD and normal participants. Firstly, we used the ssGSEA algorithm to rank individual sample gene expression by absolute value, compute cumulative empirical probability distributions for specific gene sets, and finally obtain enrichment score (ES) representing the abundance of immune infiltrating cells. The results were shown in the heatmap (Fig. 4A). and the most significant immune related gene set based on difference analysis were ‘APC co-inhibition’, ‘DCs’, ‘Cytolytic activity’, ‘T helper cells’ and ‘CD8+ T cell’ (Table 2).

In addition, at the biological function rather than a single gene level, we used the GSVA algorithm to analyze the biological status between NAFLD and normal participants to identify key biological processes. Hallmark gene sets, which summarize and represent specific well-defined biological states or processes, were used as reference datasets. The difference analysis results of GSVA in the pathway between NAFLD and normal participants were visualized by volcano map and heat map (Figs. 4B–4C). According to the differential ranking of each gene set between the two groups, the most significant gene sets were ‘E2F_targets’ and ‘ANGIOGENESIS’, which play an important role in cell proliferation and metabolism (Fig. 4D). The results of ssGSEA and GSVA analyses helped us to identify immune signatures and pathway hallmarks between NAFLD and normal participants, which would prepare us to construct co-expression networks among immune cells, signaling pathways, differential genes and transcription factors.

Construction of regulatory network and screening of hub genes based on DEGs, DETFs, immune signatures and pathway hallmarks

In the previous study herein, we obtained results for DEGs, DETFs, immune signatures, and pathway hallmarks. On this basis, we used DEGs to perform correlation analysis (correlation more than 0.85; significance less than 0.05) with transcription factor/immune signatures/pathway hallmarks, respectively (Fig. 5A). The common genes of DEG, DETF, immune signatures and pathway hallmarks were screened and visualized by heat map (Fig. 5B). Then, we constructed the correlation coefficient matrix of common genes, immune signatures and pathway hallmarks, which were shown to be highly correlated with each other (Fig. 5C).

Next, to explore the hub genes in common genes, we used correlation coefficient as links and common genes as nodes to construct the co-expression network of DEGs, DETFs, immune signatures and pathway hallmarks. As shown in the network, GTF2I and CHD2 were identified as hub genes between NAFLD and normal participants (Fig. 5D). We would validate the hub genes to clarify its function in the next step.

The expression of GTF2I is higher in patients with NAFLD than healthy people

Free fatty acids (FFA, 1mM) were used to stimulate HepG2 cells for 24 h to construct an in vitro model of NAFLD. Oil red O staining confirmed that the lipid droplets increased significantly under the stimulation of FFA (Fig. S1). Meanwhile, the mRNA and protein expression of GTF2I were significantly increased in the FFA group compared with the control group. Conversely, the expression level of CHD2 was not significantly changed by FFA stimulation (Figs. 6A–6B).

Next, a total of 13 histologically validated human samples (Normal =3; NAFLD =10) were used to examine the expression levels of GTF2I and CHD2 (Fig. 6C). The results showed that the mRNA and protein expression levels of GTF2I in NAFLD patients were significantly higher than those in healthy people (Figs. 6D–6E). In contrast, the mRNA expression level of CHD2 in NAFLD patients was significantly higher, while its protein expression was not detected (Figs. 6D–6E). Furthermore, IHC was also used to verify the relationship of GTF2I and NAFLD in human liver tissue samples. Compared with healthy individuals, GTF2I was significantly upregulated in NAFLD patients and appeared positive in the nucleus (Fig. 6F). These results suggested that GTF2I rather than CHD2 may be playing an important role in NAFLD.

GTF2I may be involved in the regulation of immune profile by regulating chemokine receptor family

Based on bioinformatics analysis and experimental validation, we screened out that GTF2I might be an important molecule in NAFLD. We further explored the relationship of GTF2I and the immune profile.

First, based on the NAFLD patient global expression matrix, we re-divided the high and low expression groups according to the expression level of GTF2I to clarify the effect of GTF2I on different immune cells. Interestingly, high expression of GTF2I group (NAFLD) was all NAFLD patients, while the GTF2I low expression group (Normal) was all non-NAFLD patients. CIBERSORT analysis revealed that GTF2I has a significant effect on the proportion of immune cells, including NK cells, CD4 + T cells, macrophages, and dendritic cells (Fig. 7A). On this basis, we further analyze the relationship of GTF2I and immune cell surface markers to clarify the specific mechanism of GTF2I in influencing the composition of immune cells. Surprisingly, GTF2I all affected the expression of surface molecules on NK cells, CD4 + T cells, macrophages, and dendritic cells (Fig. S2). This suggested that GTF2I might play core role in the regulation of immune cells. GSEA analysis showed that GTF2I is negatively correlated with ‘Cytokine–cytokine receiver interaction’ (Fig. 7B). Consistent with this, the correlation coefficient matrix and co-regulatory network showed that GTF2I is negatively correlated with chemokine receptors (CCR) (Fig. 5C). These results suggested that GTF2I may affect the immune profile in the liver by regulating the chemokine receptor family. Along this line, we compared the expression of chemokine receptor family between the two groups. As expected, GTF2I had a significant impact on the chemokine receptor family (227/261, count of significance) (Fig. 8), which may be the key pathways of GTF2I and need to be further confirmed in future work.