A novel immune-related genes signature after bariatric surgery is histologically associated with non-alcoholic fatty liver disease

ABSTRACT Increasing evidence shows that immune-related genes (IRGs) play an important role in bariatric surgery (BS). We identified differentially expressed immune-related genes (DEIRGs) of adipose tissue after BS by analysing the two expression profiles of GEO (GSE59034 and GSE29409). Subsequently, enrichment analysis, GSEA and PPI networks were examined to identify the hub IRGs and related pathways. The performance of the signature was evaluated by area under the curve (AUC) of the receiver operating characteristic (ROC). CIBERSORT algorithm was used to evaluate the relative abundance of infiltrated immune cells.42 DEIRGs were found between the GSE59034 and GSE29409 datasets. The AUC of the signature was 0.904 and 0.865 in the GSE58979 and GSE48452, respectively. Interestingly, the signature also showed good performance in diagnosing non-alcoholic fatty liver disease (NAFLD) (AUC was 0.834 and 0.800, respectively). The number of neutrophils, macrophages M2, macrophages M0 and dendritic cells activated decreased significantly. After BS, the infiltration of T cells regulatory, monocytes, mast cells resting and plasma cells in adipose tissue increased. The novel proposed IRGs signature reveals the underlying immune mechanism of BS and is a promising biomarker for distinguishing the severity of NAFLD. This will provide new insights into strategies for treating obesity and NAFLD.


Introduction
Globally, from 1975 to 2016, the number of obese patients (age >5 years) increased from 111 million to 797 million [1]. Obesity has clearly become a global burden with serious public health implications, and obesity-related metabolic diseases are spreading rapidly as obesity increases [2,3]. Obesity is closely related to non-alcoholic fatty liver disease (NAFLD) [4]. NAFLD is a clinicopathological stage from hepatic steatosis to non-alcoholic steatohepatitis (NASH) and cirrhosis [5]. With the rise of obesity, the prevalence of NAFLD has risen sharply [6]. Currently, the global prevalence of NAFLD can reach 24% [7]. Adipose tissue is an endocrine organ involved in the regulation of energy and inflammation [8]. Study has shown that adipose tissue inflammation is a prerequisite for the development of NAFLD in mice [9]. The immune cells infiltrating into fat can regulate the production and secretion of proinflammatory and anti-inflammatory factors in adipose tissue to induce NAFLD [9,10].
Currently, bariatric surgery (BS) is a more effective treatment for morbid obesity and NAFLD than nonsurgical treatment [11]. Some studies have also reported changes in the metabolic and immune status of adipose tissue following BS [10,12]. BS can promote the significant transformation of adipose tissue from pro-inflammatory state to anti-inflammatory state, and ameliorate adipose tissue function [13]. In addition, BS results in sustained weight loss and may reduce liver fat, inflammation, and fibrosis, reversing changes in hepatic pathology in patients with NAFLD and NASH [14]. However, the relationship between adipose tissue and NAFLD after BS remains unclear. It can be concluded from previous studies that immune-related genes (IRGs) play an important role in BS. Therefore, the analysis of IRGs in adipose tissue after BS is of great significance for exploring the potential value of new biomarkers for NAFLD.
With the rapid development of gene sequencing technology, IRG-based signatures have been constructed and validated in different types of diseases [15,16]. Some of these IRGs signatures have good sensitivity and specificity in predicting patient prognosis and can be used as potential tools to guide individualized therapy. However, IRGs signatures for clinical application have not been developed for BS.
In this study, we investigated the potential immune mechanism of BS in ameliorating NAFLD through a comprehensive bioinformatics analysis of adipose tissue. A novel IRGs signature was constructed and used to predict the histological degree of NAFLD. This study provides new insights into the immune regulation mechanism of BS and provides information for personalized diagnosis and treatment of NAFLD.

Microarray data
The gene expression profile was obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo). Gene expression profiles GSE59034 [17], GSE29409 [18], GSE83452 [19], GSE58979 [20] and GSE48452 [21] were selected from the GEO database (Table 1). The platform of GSE59034 was GPL11532 [HuGene-1_1-st] Affymetrix Human Gene 1.1 ST Array (Affymetrix, Inc., Santa Clara, CA). A total of 32 adipose tissue samples were obtained from patients before BS (n = 16) and patients after BS (n = 16) for subsequent analysis. The platform of GSE29409 was GPL7020, NuGO array (human) NuGO_Hs1a520180. Adipose tissues were obtained from five before BS patients and five after BS patients, with a total of 10 samples for subsequent analysis.

Identification of differentially expressed genes (DEGs)
GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/) is the tool used to identify DEGs in the GEO database. GEO2R was used to identify DEGs before and after BS. |log FC| >0.5 and p < 0.05 were used as the threshold for DEGs screening.
The overlapping DEIRGs were selected from DEGs and IRGs for further analysis. The online tool Venny (https://bioinfogp.cnb. csic.es/tools/venny/index.html) was used to select the common DEIRGs of the two datasets.

Functional and pathway enrichment analysis
DAVID database (DAVID, http://david.abcc.ncifcrf. gov/) is considered to be the most commonly used function annotation tool [22]. David was used for biofunctional analysis, including Gene Ontology (GO) analysis and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway enrichment analysis. DEIRGs were uploaded to investigate their potential functions. The cut-off criteria for p-value and false discovery rate (FDR) were 0.05 for both KEGG and GO analysis.

Gene Set Enrichment analysis (GSEA)
Further GSEA were performed on all IRGs to detect significant differences in biological functions resulting from BS. The R package 'fgsea' was used for GSEA, and 0.05 was the cut-off P value of the GSEA [23]. Ten thousand times per analysis was performed for gene set.

Protein-Protein Interaction (PPI) network analysis and identification of hub IRGs
Hub genes are generally considered to be functionally critical among other genes. We uploaded the DEIRGs to the STRING database (http://www.string-db.org/) and then selected the confidence >0.4 as the cut-off criterion for PPI network analysis. Cytoscape software was used to plot PPI network. Cytohubba plug-ins in Cytoscape was used to study key nodes in the PPI network [24]. The degree method was used to identify the hub genes in PPI network. The top 10 genes were identified as hub IRGs.

IRGs signature construction and validation
Hub IRGs between preoperative and postoperative BS were identified from GSE59034 and GSE29409. Logistic regression was applied to establish IRGs signature for predicting the histological severity of NAFLD, and all hub IRGs were covariates. GSE58979 (n = 53) and GSE48452 (n = 32) were used for signature validation.
According to the histological severity of NAFLD, 53 patients from GSE58979 were divided into mild steatosis group (group I and group II) and severe steatosis group (group III and group IV). The performance of the signature was evaluated by area under the curve (AUC) of the receiver operating characteristic (ROC) and consistency index (C index). The Hosmer-Lemeshow test was used to evaluate goodness-of-fit of the signature. Decision curve analysis (DCA) was used to quantify the net benefit of patients under different threshold probabilities to determine the clinical value of signature.

Construction of a diagnostic signature for NAFLD
To investigate whether IRGs signature could be used to predict the occurrence of NAFLD, GSE48452 (n = 73) and GSE83452 (n = 76) were used to validate the signature. The predictability of the signature was then evaluated by AUC of ROC and DCA. Goodness-of-fit was examined by using the Hosmer-Lemeshow test.

Comparison of immune cells infiltration in preoperative and postoperative BS
CIBERSORT algorithm was used to evaluate the relative abundance of infiltrated immune cells before and after BS. CIBERSORT can estimate the composition of 22 immune cell types, mainly including plasma cells, macrophages, eosinophils, natural killer cells, B cells, dendritic cells, T cells, and neutrophils [25]. The R package 'CIBERSORT' was used. Wilcox test was used to compare the differences of immune cell types between the pre-BS and the post-BS.

DEIRGs after bariatric surgery
1162 DEGs were identified from the GSE59034 dataset, including 272 upregulated DEGs and 890 downregulated DEGs with |log FC| >0.5 and p < 0.05. Six hundred and fifty-eight DEGs were identified from the GSE29409 dataset, including 184 upregulated DEGs and 474 downregulated DEGs with |log FC| >0.5 and p < 0.05. IRGs were extracted from the ImmPort database. Next, DEIRGs were determined between the extracted IRGs and DEGs. A total of 164 DEIRGs were identified from the GSE59034 dataset (Figure 1a), which 17 were upregulated and 147 were downregulated. A total of 105 DEIRGs were identified from the GSE29409 dataset (Figure 1b), which 10 were upregulated and 95 were downregulated. We used the online tool Venny to select 42 co-expressed DEIRGs ( Figure 1c and Table 2).

Biofunctional analysis of the DEIRGs
The biological functions of 42 DEIRGs were studied by GO and KEGG analysis ( Figure 2). In GO analysis, for molecular function (MF), these genes were significantly enriched in receptor activity, protein tyrosine kinase activity, and RAGE receptor binding. For biological process (BP), these genes were significantly enriched in signal transduction, G-protein coupled receptor signalling pathway, chemotaxis, positive regulation of cell proliferation and neutrophil chemotaxis. And they played a key role in host inflammatory response, immune response, and innate immune response. For cellular component (CC), these genes were enriched in plasma membrane, integral component of membrane, external side of plasma membrane, extracellular space, extracellular region, mast cell granule, membrane raft, and Toll-like receptor 1-Toll-like receptor 2 protein complex. KEGG analysis (Figure 2d) found that DEIRGs were enriched in immune pathways, such as cytokinecytokine receptor interaction, chemokine signalling pathway, JAK-STAT signalling pathway, natural killer cell mediated cytotoxicity, and Fc epsilon RI signalling pathway.

Potential signalling pathways were identified by GSEA
GSEA was performed to determine the underlying mechanism by which BS may relieve inflammation. According to its normalized enrichment score, the most significant enriched signalling pathway was identified. Correlative pathways, including graft-versus-host disease, Fc gamma R-mediated phagocytosis, toll-like receptor signalling pathway, type 1 diabetes mellitus, leukocyte transendothelial migration, chemokine signalling pathway, viral myocarditis, haematopoietic cell lineage, antigen processing and presentation, cytokine-cytokine receptor interaction, and renin angiotensin system were highly enriched after BS ( Figure 3 and Table 3).

Performance of the IRGs signature
The AUC of the IRGs signature was 0.904 (GSE58979) and 0.865 (GSE48452), respectively (Figure 5a, b). osmer--Lemeshow test was not significant (0.554 and 0.710, respectively), indicating that the signature had a good fit. The IRGs signature indicated favourable prediction with a C index of 0.904 and 0.865 (training and validation cohort, respectively). DCA had shown that IRGs signature could obtain a high net benefit in predicting histological severity of NAFLD within the most reasonable threshold probability ranges (Figure 5c, d).

Immune cell infiltration in preoperative and postoperative BS
We measured the relative abundance of 22 different immune cells in each patient before and after BS using CIBERSORT. Comparison of immune cell composition in adipose tissue before and after BS was shown in Figure 6a. The results showed that neutrophils, macrophages M2, macrophages M0 and dendritic cells activated decreased significantly after BS. However, T cells regulatory (Tregs), monocytes, mast cells resting and plasma cells infiltrated more in the adipose tissue after BS (Figure 6b).

Discussion
Obesity can cause chronic low-grade inflammation of adipose tissue, which is manifested in various systems [26]. Obesity causes changes in adipokine production and other tissue inflammation through the proliferation and hypertrophy of adipocytes, as well as an increase in free fatty acids and extensive tissue remodelling [27]. Notably, inflammation of adipose tissue may be a prerequisite for NASH development [20]. BS is an effective treatment for obesity, reducing inflammation and NAFLD, and improving long-term survival [13,28]. However, the molecular mechanism and clinical significance of BS in improving NAFLD remain to be explored to a certain extent. Therefore, it is of profound significance to search for specific diagnostic markers and analyse the pattern of immune cell infiltration after BS to ameliorate the prognosis of NAFLD patients.
In this study, we identified 42 DEIRGs from two GEO datasets (GSE59034 and GSE29409), revealing the potential immune pathways that BS reduces adipose tissue inflammation. We developed and validated a new signature based on 10 IRGPs in multiple independent cohorts. The signature had a reliable predictive value and accuracy for predicting the histological severity of NAFLD. The results showed that BS can cause the  changes of IRGs and a variety of signalling pathways in adipose tissue, and such changes may be involved in the regulation of the development of NAFLD.
Functional enrichment analysis to clarify the biological function of DEIRGs. DEIRGs were primarily enriched in chemokine signalling pathway, G-protein coupled receptor signalling pathway, JAK-STAT signalling pathway, natural killer cell mediated cytotoxicity, and Fc epsilon RI signalling pathway, etc. Confirming with previous evidence, obesity caused a low degree of chronic inflammation characterized by infiltration of B-cells, neutrophils, natural killer T cells, Th1 CD4 + T cells, macrophages, and CD8 + T cells into adipose tissue [29]. In addition to the above pathways, GSEA indicated that BS-associated IRGs were mainly enriched in Toll-like receptor signalling pathway, leukocyte transendothelial migration, antigen processing and presentation, and renin angiotensin system. These pathways have been shown to be important mechanisms of adipose tissue inflammation and to be involved in the pathogenesis of NAFLD [30][31][32][33]. More and more evidence showed that chemokine signalling pathway played an important role in the pathogenesis of NAFLD and was an important factor linking obesity and NAFLD [33]. JAK-STAT signalling pathway was involved in the regulation of inflammation associated with metabolic abnormalities in adipose tissue [34]. Mice with deficiency of Jak3 and Stat6 can display obesity and liver steatosis [35,36]. Other pathways such as Toll-like receptor, G-protein coupled receptor and natural killer cell mediated cytotoxicity were important trigger factors to promote chronic inflammation of adipose tissue and related metabolic disorders, and inhibition of these pathways played a protective role in NAFLD [30][31][32].
By constructing PPI network and Cytohubba analysis in Cytoscape, 10 hub IRGs were identified as IL6, TLR8, TLR2, IL10, TLR1, CXCR4, HCK, FCER1G, LYN and SYK. These genes were all downregulated after BS. The results were similar to previous studies [13,37,38]. IL-6 from adipose tissue, which travels through the portal vein to the liver, promotes NAFLD development [39]. IL10 is an anti-inflammatory cytokine released by monocytes/macrophages in adipose tissue. With the decrease of IL6 after BS, IL10 as the antagonistic anti-inflammatory cytokine of IL6 also decreased [40]. Abnormal activation of Toll-like receptors (TLRs) had been found to cause obesity-related systemic inflammation and associated comorbidities in obese rats [41]. Inhibition of TLRs could reduce the activation of MAPKs and NF-κB and ameliorate obesity-induced NAFLD [42]. CXCR4 was a G proteincoupled chemokine receptor and plays a key role in improving NAFLD after BS [43]. By increasing the affinity of CXCR4, CD4 + T-cells deposition were increased in the livers of NAFLD patients [44]. Previous study had shown that FCER1G played a key role in the occurrence and progression of NAFLD [38]. Inhibition of SYK could ameliorate the inflammation and steatosis of NAFLD by reducing the activation of macrophages and the production of CCl2, IL-6 and TNF-α [37]. Notably, FCER1G, LYN and SYK downregulation were able to reduce obesity induced by highfat diet and the growth of lipid droplets in adipocytes  [45,46]. This may be the underlying mechanism by which BS can ameliorate NAFLD.
NAFLD is one of the leading causes of liver transplantation and promotes the development of hepatocarcinoma, and BS appears to be the most durable and effective treatment for NAFLD [7,47]. BS can manipulate adipocyte-derived adipocytokines (adiponectin, leptin, IL-6, etc.) to regulate NAFLD [43]. For example, BS could inhibit the activation of Kupffer cells and hepatic stellate cells (HSC) and reduces the expression of inflammatory genes by increasing adiponectin levels [48]. In addition, Johannie DP et al. found that adipose tissue may directly promote the progression of NAFLD through inflammatory cytokines and immune cells, and they built a model using five genes to accurately predict the histological severity of NAFLD [20]. Fu C et al. recently analysed the liver tissue after BS by bioinformatic method, and they found the potential key genes and pathways for Roux-en-Y gastric bypass to ameliorate NASH [43]. However, the liver tissue samples collected include hepatocytes and mesenchymal cells, which will affect the cellular composition of the liver tissue. It has been concluded from previous studies that IRGs play an important role in obese adipose tissue and NAFLD, and inflammation in adipose tissue may be a prerequisite for NASH development. Therefore, we selected adipose tissue after BS as the training cohort to construct the signature, and liver tissue from another gene profile was selected as the validation cohort for verification. This signature had good performance in diagnosing NAFLD and its histological severity. DCA has shown that when determining NAFLD and its histological severity within the most reasonable threshold probability range, IRGs signature can provide a relatively high overall net benefit. At the same time, it also suggests that BS may ameliorate NAFLD by reducing the degree of inflammation in adipose tissue.
Increasing evidence supports the fact that immune cell infiltration plays a significant role in the initiation and progression of adipose tissue inflammation in obesity [10,49]. To further investigate the infiltration of immune cells in adipose tissue after BS, CIBERSORT algorithm was used to calculate the changes in the composition of immune cells. Neutrophils can recruit macrophages and dendritic cells to adipose tissue and play a major role in the initiation of inflammation [10]. Hypertrophy of adipocyte can secrete MCP-1/CCL2, promote the macrophages collect and produce proinflammatory cytokines TNFα, IL-6 and IL-10, resulting in NAFLD [43]. Dendritic cells in adipose tissue induce Th17 to secrete cytokines IL-6, IL-10 and IL-23, leading to more inflammation [50]. Tregs, a type of anti-inflammatory immune cells, were found that the reduction of Tregs in adipose tissue may increase the severity of NAFLD [51,52]. The above immune cell infiltration, IRGs and related pathways indicate that immune process is an important pathway for BS to treat obesity and NAFLD.
However, there were several notable limitations to this study. First, because we used data accessed from public databases, we were unable to obtain sufficient demographic and clinical information. Gender, body weight, BMI, surgical methods and other factors may influence the immune regulation mechanism of BS. Second, the follow-up time of the two datasets in this study was different. GSE59034 and GSE29409 were followed up for 2 and 1 years, respectively. In this study, some biological information may have been overlooked due to a lack of detailed consideration of follow-up time. Finally, although BS can ameliorate NAFLD, not all patients can be completely cured. The IRGs signature of this study was not used to predict the prognosis of NAFLD after BS. This limited the clinical application of the signature in this study.
In conclusion, this study analysed the potential immune mechanism of adipose tissue after BS and revealed the change of immune cells infiltration in adipose tissue after BS. In addition, this study constructed a novel IRGs signature that could accurately distinguish the severity of NAFLD. This study suggests that adipose tissue may play an important role in ameliorating NAFLD after BS.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This research did not receive any specific grant from funding agencies in the public,commercial, or not-for-profit sectors.

Data availability statement
The data that support the findings of this study are available in GEO database (http://www.ncbi.nlm.nih.gov/geo), reference number [GSE59034, GSE29409, GSE83452, GSE58979 and GSE48452].