Identification of Cuproptosis-Related Genes in Nonalcoholic Fatty Liver Disease

Nonalcoholic fatty liver disease (NAFLD) is the most prevalent hepatic pathology worldwide. However, the precise molecular mechanisms for NAFLD are still not sufficiently explained. Recently, a new mode of cell death (cuproptosis) is found. However, the relationship between NAFLD and cuproptosis remains unclear. We analyzed three public datasets (GSE89632, GSE130970, and GSE135251) to identify cuproptosis-related genes stably expressed in NAFLD. Then, we performed a series of bioinformatics analyses to explore the relationship between NAFLD and cuproptosis-related genes. Finally, 6 high-fat diet- (HFD-) induced NAFLD C57BL/6J mouse models were established to carry out transcriptome analysis. The results of gene set variation analysis (GSVA) revealed that the cuproptosis pathway was abnormally activated to a certain degree (p = 0.035 in GSE89632, p = 0.016 in GSE130970, p = 0.22 in GSE135251), and the principal component analysis (PCA) of the cuproptosis-related genes showed that the NAFLD group separated from the control group, with the first two principal components accounting for 58.63%-74.88% of the variation. Among three datasets, two cuproptosis-related genes (DLD and PDHB, p < 0.01 or 0.001) were stably upregulated in NAFLD. Additionally, both DLD (AUC = 0.786–0.856) and PDHB (AUC = 0.771–0.836) had favorable diagnostic properties, and the multivariate logistics regression model further improved the diagnostic properties (AUC = 0.839–0.889). NADH, flavin adenine dinucleotide, and glycine targeted DLD, and pyruvic acid and NADH targeted PDHB in the DrugBank database. The DLD and PDHB were also associated with clinical pathology, especially with steatosis (DLD, p = 0.0013–0.025; PDHB, p = 0.002–0.0026) and NAFLD activity score (DLD, p = 0.004–0.02; PDHB, p = 0.003–0.031). What is more, DLD and PDHB were correlated with stromal score (DLD, R = 0.38, p < 0.001; PDHB, R = 0.31, p < 0.001) and immune score (DLD, R = 0.26, p < 0.001; PDHB, R = 0.27, p < 0.001) in NAFLD. Furthermore, Dld and Pdhb were also significantly upregulated in the NAFLD mouse model. In conclusion, cuproptosis pathways, especially DLD and PDHB, could be potential candidate genes for NAFLD diagnostic and therapeutic options.


Introduction
Nonalcoholic fatty liver disease (NAFLD), as a metabolic disease, is the most prevalent hepatic pathology worldwide with a prevalence of approximately 25% [1]. It ranges from simple steatosis (SS) to advanced stage of nonalcoholic steatohepatitis (NASH), the latter rapid progression toward liver cirrhosis, and hepatocellular carcinoma (HCC) [2,3]. Cur-rently, sensitive biomarkers for the diagnosis of NAFLD are still lacking, and liver biopsies, as an invasive method, are still considered the gold standard for diagnosis and prognosis [4]. Although numerous studies have attempted to study the pathogenesis and progression of NAFLD, there are still no effective drugs for NAFLD other than lifestyle changes [5]. Moreover, the development of NAFLD is a complex process and is still not sufficiently explained [6].
Consequently, it is crucial to explore the mechanisms involved in the pathogenesis of NAFLD to identify new potential targets for diagnosis and therapy.
There has been recognition that copper metabolism disorder can lead to a variety of chronic liver diseases, such as Wilson disease (WD), NAFLD, and liver cirrhosis [9]. Zhang et al. indicated that serum copper levels were positively correlated with body mass index (BMI), leptin, and insulin resistance, which are all risk factors for NAFLD [2]. Dev et al. found a dramatic increase in hepatic copper levels, resulting in obesity and hepatic steatosis in the hepatocytespecific knockout of Atp7b WD mouse model [10].
In this research, we seek to comprehensively investigate the molecular alterations and clinical relevance of the cuproptosis-related genes in NAFLD. This study highlights the importance of cuproptosis-related genes in NAFLD and lays a foundation for future studies of cuproptosis in NAFLD.

Materials and Methods
2.1. Data Retrieving and Processing from GEO. Data from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih .gov/geo) fulfilled the inclusion criteria below: ① publication date from 2010 to 2022; ② containing NAFLD and normal tissue samples; ③ providing detailed clinicopathological information. The exclusion criteria were as follows: ① duplicated research; ② patients who underwent bariatric surgery or with severe diseases; ③ data that could not be analyzed; ④ animal or cell experiments. Subsequently, the gene expression profiles of GSE89632 [11], GSE130970 [12], and GSE135251 [13] were downloaded from GEO. Eventually, 24 control patients, 19 SS patients, and 19 NASH patients in GSE89632, 6 control patients and 72 NAFLD patients in GSE130970, 10 control patients, and 216 NAFLD patients were included in this study (Table 1). Since the datasets were from a public database, patient consent and ethics committee approval were not required. Then, the gene array data (GSE89632) was converted into the quantile normalized values, and the read count data (GSE130970 and GSE135251) was standardized by the transcripts per million (TPM). The overall research design is shown in Figure 1.

Gene Set Variation Analysis (GSVA) and Single-Sample
Gene Set Enrichment Analysis (ssGSEA). GSVA, a pathway enrichment method that estimated variation of pathway activity, was performed to evaluate the role of the cuproptosis pathway in NAFLD using R package "GSVA" [14]. In addition, the stromal score, immune score, and immune cells' marker enrichment were ssGSEA and they were calculated by "estimate" (http://bioinformatics.mdanderson.org/ estimate/), and "GSVA" R packages. The immune gene sets were downloaded from Charoentong et al. [15]. Then, the spearman correlation analysis between DLD/PDHB expression and immune cells was performed. The results were visualized using the "ggplot2" R package.

Correlation Analysis of Cuproptosis-Related Genes,
Protein-Protein Interaction (PPI) Network Construction, and the Prediction of Potential Drugs. The "ggcorrplot" R package was used to recognize the correlation between cuproptosis-related genes by the Spearman correlation analysis. The STRING database (https://string-db.org/) was utilized to construct a PPI network with an interaction score > 0:4. The prediction of potential drugs was performed in the Drug-Bank database (https://go.drugbank.com).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Analyses.
To explore the potential pathways of DLD and PDHB, the GO and KEGG analyses were performed using "clusterProfiler" [16] and "org.Hs.eg.db" R package. The median of the risk model scores of DLD and PDHB was selected as the cut-off value.

Animal Model and Experiment Design.
A total of 6 male C57BL/6J mice (weight, 23:47 ± 1:18 g; age, 6 weeks) were purchased from the medical laboratory animal center of Guangdong (Guangzhou, China). All the mice were acclimatized under a temperature of 24 ± 2°C, a relative humidity of 55 ± 10%, and a 12 h light/dark cycle for 10 days before the commencement of the animal experiment. All animal experiments were approved by the experimental animal ethics committee of Jinan University and were performed in accordance with the Guidelines for Care and Use of Laboratory Animals of Jinan University (IACUC-20190916-09).
After acclimation, the 6 mice were randomly divided into the normal control group (NC, n = 3) and the NAFLD group (n = 3). The NC group mice were fed with standard mouse diet, and the NAFLD group mice were induced by a high-fat diet (HFD) (containing 34% fat, 2% cholesterol, 26% carbohydrate, 26% protein, and 12% basic feed (w/w)) for 8 weeks before sacrificing.
2.6. Histological Analysis. At the end of the study, the mice were sacrificed to measure liver mass and liver mass index by using the following formula: live mass index = liver mass /body mass × 100%. Afterward, the liver samples were immersed in 10% formalin neutral buffer solution for 48 h, then processed routinely, embedded in paraffin, sectioned to 5 μm thickness and stained with hematoxylin and eosin (H&E). Lipid accumulation in liver was analyzed by Oil red O (ORO) staining (Sigma). Slides were observed with a light microscope (Leica DMi 8).

2
Oxidative Medicine and Cellular Longevity . Sequences were aligned to the mm10 mouse reference genome using STAR Aligner. Then, reads aligning to genes were counted using htseq-count, and analysis of differentially expressed genes (DEGs) was performed using the "limma" R package [17] and later visualized by volcano plot using "ggplot2" R package (https://ggplot2.tidyverse.org/). The threshold for the DEGs was set as p value < 0.05 and jlog 2 fold change ð FCÞj ≥ 0:5. The raw data was stored in the ArrayExpress database (accession: E-MTAB-11980, https://www.ebi.ac .uk/arrayexpress/).  Figure 1: The overall research design. The data were downloaded from the GEO database. Then, GSVA and PCA were performed to identify the activation of the cuproptosis pathway in NAFLD. The differential expression analysis (DEA) and Venn diagram showed that DLD and PDHB were stably upregulated in three datasets. Subsequently, correlation analysis and PPI were conducted to explore their relationship. Clinical pathology analysis and risk model construction were carried out. Next, the immune microenvironment and potential drugs for DLD and PDHB were explored. Finally, the NAFLD mouse model was used to further verify the expression of DLD and PDHB.  5 Oxidative Medicine and Cellular Longevity component analysis (PCA) was applied to display the overall differences of cuproptosis-related genes and visualized by the "ggplot2" R package. Venn diagram was performed using the "ggvenn" R package. The multivariate logistics regression risk model was constructed in DLD and PDHB. A receiver-operating characteristic (ROC) curve was conducted to assess the diagnostic value of the genes and the model by "pROC" R package and visualized by the "ggplot2" R package. A difference with p < 0:05 was considered significant.

Cuproptosis Pathway Plays a Role in NALFD.
To determine whether the cuproptosis pathway played a role in NAFLD, GSVA and PCA were performed. The results showed that the enrichment score of the cuproptosis pathways in GSE89632 and GSE130970 was significantly increased (p < 0:05) in the NAFLD group compared with the control group, but not in GSE135251 (p = 0:22, Figures 2(a)-2(c)). Next, the principal component analysis (PCA) of the cuproptosis-related genes showed that the NAFLD group separated from the control group, with the first two principal components accounting for 58.63%-74.88% of the variation (Figures 2(d)-2(f)). These results suggest that the cuproptosis pathway is activated to some degree in NAFLD.

Differential Expression of Cuproptosis-Related Genes.
Subsequently, the expression levels of the cuproptosis-related genes were explored in the three datasets, respectively. The cuproptosis-related genes were changed to varying degrees in the three datasets, in which DLD (p < 0:001 in GSE89632, p = 0:004 in GSE130970, p < 0:001 in GSE135251) and PDHB (p = 0:001 in GSE89632, p = 0:007 in GSE130970, p = 0:001 in GSE135251) were stably and significantly upregulated in the three datasets (Figures 3(a)-3(d)). The outcomes illustrate that DLD and PDHB play a part in the pathogenesis and progression of NAFLD.

DLD and PDHB Are Associated with Metabolic
Pathways and Copper Toxicity. Subsequently, GO and KEGG analyses were implemented for pathway analysis by comparing the low-and high-risk scores of DLD and PDHB in GSE130970. The GO results showed that the high-risk score of DLD and PDHB was associated with the acetyl-CoA metabolic process and fatty acid metabolic process, and detoxification of copper ion pathway was enriched in the low-risk score group (Figure 5(a)). Besides, the KEGG results illustrated that the high-risk score of DLD and PDHB was associated with fatty acid metabolism, glycolysis/gluconeogenesis, and pyruvate metabolism ( Figure 5(b)). These results indicate that high-expressed DLD and PDHB appear to increase copper toxicity and are related to the acetyl-CoA metabolic process and pyruvate metabolism, both important substrate sources for the TCA cycle.

Further Validation of DLD and PDHB in the NAFLD
Mouse Model. Mice fed HFD develop hepatic steatosis, mimicking the NAFLD of humans [18]. We conducted an animal experiment, and after HFD was fed to NAFLD mouse group for 8 weeks, a significant increase in the body mass, liver mass, and liver mass index was observed in the NAFLD   14 Oxidative Medicine and Cellular Longevity group compared with the NC group (Table 2). Besides, the results of liver pathology showed that liver sections from the NALFD group had hepatocyte swelling, ballooning degeneration, and different sizes of lipid droplets (Figures 7(a) and 7(b)). Oil Red O staining revealed the accumulation of neutral lipids in NAFLD (Figures 7(c) and 7(d)). These results demonstrate that the NAFLD mouse model was successfully established. Subsequently, transcriptome analysis of mice liver tissues indicated that Dld and Pdhb were also significantly upregulated in the NAFLD group compared with the NC group (Figure 7(e)).

Discussion
In this study, we first explored the relationship between 10 cuproptosis-related genes and NAFLD. Then, DLD and PDHB were stably upregulated in NAFLD among the three datasets. Besides, the DLD and PDHB were also associated with the clinical characteristics (especially steatosis and NAS) and the immune microenvironment in NAFLD. NADH, FAD, and glycine targeted DLD, and pyruvic acid and NADH targeted PDHB. In addition, the high-expressed DLD and PDHB appeared to increase copper toxicity and had impacts on the TCA cycle by affecting the acetyl-CoA and pyruvate, further suggesting that the cuproptosis might affect the development of NAFLD. Furthermore, the NAFLD mouse model was established to determine the expression signature of Dld and Pdhb. These outcomes suggest that DLD and PDHB promoted hepatic steatosis and trigger liver inflammation through the cuproptosis.
Copper, a cofactor for various enzymes, is an essential micronutrient required for normal cell function, and the cuproptosis is related to various cancer progressions, such as HCC [19,20]. Recent studies have determined that serum and liver copper are related to NAFLD, but the specific mechanism is still unclear [21]. Mitochondrial respiration is required for cuproptosis, but upregulation of mitochondrial activity often occurs in an early NAFLD stage, and the expression levels of cuproptosis-related genes (DLD and PDHB) are increased in our study, which suggests that cuproptosis may contribute to the progression of NAFL to NASH [7,22]. However, the mitochondrial function gradually decreases during the progression of NAFLD, suggesting that the effect of cuproptosis is progressively attenuated and  may lead to HCC eventually [22]. In addition, it is reported that suppression of DLD expression inhibits melanoma growth and tumor proliferation by increasing intracellular reactive oxygen species (ROS) production and thereby inducing autophagy cell death [23]. PDHB catalyzes the conversion of pyruvate to acetyl-CoA, acting as a central node that links glucose metabolism, lipid metabolism, and the TCA [24]. The dysfunction of PDHB can lead to metabolism alteration which is one of the hallmarks of cancer cells [25,26]. Recent studies show that glycine and hepatic NAD + decrease in NAFLD individuals, and their supplementation can ameliorate NAFLD, but their mechanisms have not been completely elucidated, and our study may provide new insight for them [27][28][29]. Besides, the altered pyruvate metabolism, particularly enhanced lactate production, plays an essential role in the NAFLD progression, but further studies are needed to determine whether pyruvic acid supplementation improves NAFLD [30]. FAD, a redox-active coenzyme, is involved in various metabolic pathways, including the beta-oxidation of fatty acid and TCA, and its role in NAFLD is still unknown [31]. Moreover, DLD and PDHB are positively correlated with gamma delta T cell and CD4+ T cell which often accumulate in the liver and subsequent stimulate inflammatory processes in NAFLD [32,33]. Hence, the cuproptosis, especially DLD and PDHB, may also play a vital role in NAFLD. The present study had several advantages. Firstly, this was the first study to explore the relationship between NAFLD and cuproptosis-related genes. Secondly, we comprehensively analyzed multiple datasets and screened out stably upregulated cuproptosis-related genes (DLD and PDHB), which were later verified in the NAFLD mouse model, whereas our study also had some limitations. First of all, we did not carry out an in-depth study of cuproptosis on NAFLD. For another, the NAFLD mouse model rather than human tissue was used in this study, which might have influenced the results. But HFD-induced NAFLD mouse model can mimic the NAFLD of humans most accurately [18].
In conclusion, our study presented a systematic analysis of molecular alterations and interactive genes of cuproptosis in NAFLD. Finally, we screened out two cuproptosis-related genes (DLD and PDHB) that were correlated with NAFLD prognosis. Although further research is still needed, we provide useful and novel information to explore the potential candidate genes for NAFLD diagnostic and therapeutic options.

Ethical Approval
All animal experiments were approved by the experimental animal ethics committee of Jinan University and was per-formed in accordance with the Guidelines for Care and Use of Laboratory Animals of Jinan University (IACUC-20190916-09).