Introduction
Hepatocellular carcinoma (HCC) is the second most deadly cancer with a low 5-year survival rate and poor prognosis.1 Previous studies have shown that liver cancer has an incidence of approximately 850,000 cases annually.2 It is well-known that the primary risk factors for liver cancer include chronic hepatitis B virus (HBV) infection, hepatitis C virus (HCV) infection, water pollution, and chemical carcinogens. Patients with liver cirrhosis can progress to HCC. Hence, the prevention and prediction of liver cirrhosis are crucial to prevent liver cancer development. A recent transcriptome meta-analysis indicated that targeting the lysophosphatidic acid pathway was the primary chemical preventive measure in >500 patients with liver cirrhosis.3 Of note, the survival rate of liver cancer has decreased significantly recently, although the 5-year survival rate of advanced liver cancer remains approximately 5%.4 The prognosis for patients with HCC remains controversial. Therefore, the discovery of biomarkers for the prognosis of HCC is urgently required.5,6
Ferroptosis is a novel form of a regulated cell death pathway, and its discovery has helped to explore the role of iron metabolism in the development and progression of cancers. It is characterized by the excessive accumulation of lipid peroxidation products that cause oxidative damage, which is dependent on iron ions and the mitochondrial chain.7 Several factors contribute to ferroptosis. First, glutathione peroxidase (GPX4) can convert lipid peroxides into non-toxic lipid alcohols and its deficiency results in the accumulation of reactive oxygen species (ROC) and a decline in the antioxidant capacity in cells, leading to ferroptosis and oxidative cell death.8 In addition, ferroptosis is regulated by defective oxidoreductase FSP1 (AIFM2), GCH1, and its product tetrahydrobiopterin (BH4). Functionally, ferroptosis contributes to the pathogenesis of diseases, including various cancers. Significantly, the nuclear receptor coactivator 4 (NCOA4) is one of the selective cargo receptors for ferritinophagy. NCOA4 overexpression increases the degradation of ferritin and promotes ferroptosis.9,10 Previous research found that many types of tumor cells are prone to drug-induced ferroptosis.11 Therefore, it is essential to explore and determine the ferroptosis-related biomarkers for the prognosis of HCC patients.
Recently, target therapy has become attractive and benefits HCC patients.12 Sorafenib can inhibit the growth of HCC by targeting the Raf/MEK/ERK signaling,13,14 and attenuates tumor angiogenesis by blocking the vascular endothelial and platelet-derived growth factor receptors. Immunotherapy by targeting immune checkpoints is a promising method to kill tumor cells by enhancing T cell activity. Immune checkpoint inhibitors enhance T cell immunity against tumors and benefit patients with hematologic malignancies and those with some types of solid tumors.15 Therefore, it is vital to discover the immune biomarkers for the response to immunotherapy and the prognosis of HCC patients.
In the present study, the ferroptosis-related genes (FRGs) in HCC patients were screened and identified in the Cancer Genome Atlas (TCGA) database and a prognostic signature was constructed using 10 genes for the prognosis of HCC patients, which was validated in the International Cancer Genome Consortium (ICGC) database. Furthermore, the signature was used to establish a model and determine its value in the prognosis of HCC patients. Then, the immunotumor environment of HCC was analyzed. Finally, a nomogram was constructed using the signature genes to predict the 1-, 3-, and 5-year overall survival (OS) of HCC. Our data indicated that the model was valuable for the prognosis of HCC. Our findings might help to predict the therapeutic response to immunotherapies in HCC patients and provide new insights into the pathogenic role of ferroptosis in the pathogenesis of HCC.
Materials and methods
Data collection
The demographic, clinical, and transcriptomic data of HCC patients were acquired from the TCGA database, which included 374 HCC and 50 adjacent non-tumor samples. The data was validated from the ICGC portal and these samples were from patients with chronic HBV or HCV infection in Japan. The data were analyzed using the algorithm in the limma R package.16
Identification of differentially expressed genes
The differentially expressed genes (DEGs) relative to ferroptosis were defined when one gene expression between the control and HCC tissues reached significant (|logFC| > 1 and adjusted p < 0.05) with a false discovery rate (FDR) <0.05, determined using the R package limma in R software (version 4.0.5).17
Construction and validation of a prognostic ferroptosis-related gene signature
The significance of individual DEGs was determined by univariate Cox analysis for the OS of HCC patients, and only DEGs with a p-value < 0.001 were chosen for the following analysis. Finally, the R package glmnet was applied for least absolute shrinkage and selection operator (LASSO) regression analysis, and the prognostic candidates for ferroptosis-related DEGs were obtained by tenfold cross-validation lambda (λ) value.18 The expression levels of different genes were the independent variables for LASSO regression and the survival rate and status were the dependent variables for patients with HCC. A prediction model was established for the analysis of the risk of FRGs using the Cox regression coefficient weighted estimation method.19 The prognostic FRGs signatures were manifested as:20
Risk score=∑i=1ncoefficient∗Expressionnn of FRG.
To better predict the prognosis of HCC, an ROC curve was used to predict the OS of HCC patients.21 The survival of patients with HCC in TCGA and ICGC was estimated by the Kaplan-Meier (K-M) method and analyzed by the log-rank test. Then, an independent prognostic analysis was performed to evaluate the clinical features of OS-independent prognostic factors.22 The data were further analyzed by principal component analysis (PCA) to reduce their dimensionality. In addition, the distribution of HCC patients in each group was analyzed by t-distributed stochastic neighbor embedding (t-SNE).23
Function enrichment analysis of DEGs in TCGA and ICGC cohorts
To explain the biological enrichment pathway related to the risk score, gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were conducted based on the R package clusterProfiler with high and low-risk groups of DEGs (|log2FC| ≥ 1, FDR < 0.05).16 The DEGs involved in the immune function-related pathways were further analyzed.24,25 Finally, the activity of 13 immune-related pathways and 16 types of immune infiltrates were analyzed by the single-sample gene set enrichment analysis (ssGSEA algorithm) using the gsva R package.26
Construction of predictive nomogram
A nomogram was established for the prognosis of HCC patients.27 The R package rms was used to assess the value of the nomogram in predicting the OS. Finally, the nomogram included risk score and clinical stage and was used to predict 1-, 3-, and 5-year OS of HCC patients. The calibration curves for 1-, 3- and 5-year OS were performed for verification of the power and accuracy of prognosis prediction.28 The C index was used to evaluate the performance of each model. The different models were compared using the rms R package.
Statistical analysis
The HCC patients in the TCGA training cohort were stratified into high and low-risk groups, based on the median risk score. Univariate and multivariate Cox regression analyses determined the independent prognostic parameters. The survival of HCC patients was estimated by the K-M method and analyzed by a log-rank test using the survival R package. The sensitivity and specificity of the nomogram to predict the OS of HCC patients were analyzed by ROC curves using the R package pROC. The difference between groups was analyzed by a Chi-squared test and Student’s t-test. All statistical analyses were performed using R 4.0.5, and a p-value < 0.05 was considered statistically significant.
Results
Research process and patients’ characteristics
For simplicity, a flow chart (Fig. 1) was constructed that started with a cohort of 422 HCC patients in the TCGA-LIHC and a validation cohort of 260 HCC patients in the ICGC. Their demographic and clinical data are shown in Table 1. The patients’ characteristics included age, tumor Gleason grade, gender, T, M and N stage, and survival status. In addition, 291 FRGs were acquired from FerrDb database (http://www.zhounan.org/ferrdb ) and the human gene database (Gene Cards) with the keywords “Ferroptosis” (https://www.genecards.org/ ).
Table 1The demographic and clinical characteristics of HCC patients
Case | TCGA cohort
| LIRI-TP cohort
|
---|
377 | 260 |
---|
Age(%) | | |
≤60 | 180(47.7%) | 55(21.2%) |
>60 | 197(52.3%) | 205(78.8%) |
Gender(%) | | |
Female | 122(32.4%) | 68(26.2%) |
Male | 255(67.6%) | 192(73.8%) |
Grade(%) | | |
Grade 1 | 55(14.6%) | NA |
Grade 2 | 180(47.8%) | NA |
Grade 3 | 124(32.9%) | NA |
Grade4 | 13(3.4%) | NA |
Unknown | 5(1.3%) | NA |
Stage(%) | | |
I | 175(46.4%) | 40(15.4%) |
II | 87(23.1%) | 117(45.0%) |
III | 86(22.8%) | 80(30.8%) |
IV | 5(1.3%) | 23(8.8%) |
Unknown | 24(6.4%) | 0(0.0%) |
T(Tumor)(%) | | |
T1 | 185(49.1%) | NA |
T2 | 95(25.2%) | NA |
T3 | 81(21.5%) | NA |
T4 | 13(3.4%) | NA |
N(Lymph Node)(%) | | |
N0 | 257(68.2%) | NA |
N1 | 4(1.1%) | NA |
M(Metastasis)(%) | | |
M0 | 272(72.1%) | NA |
M1 | 4(1.1%) | NA |
Survival status(%) | | |
Alive | 245(65.0%) | 214(82.3%) |
Dead | 132(35.0%) | 46(17.7%) |
Identification of prognostic ferroptosis-related DEGs in the TCGA cohort
To determine the association of ferroptosis-related DEGs with the OS of HCC patients using univariate Cox regression analysis, 32 ferroptosis-related DEGs were significantly associated with the OS of HCC patients (p < 0.001) (Table 2). Based on the criteria of FDR < 0.05 and |log2FoldChange| > 1, these FRGs were differentially expressed between 374 HCC and 50 non-tumor liver samples in the TCGA database. Further univariate Cox regression analyses indicated that these 19 ferroptosis-related DEGs were associated with the OS of HCC patients (Fig. 2e). The forest plot of hazard ratios (HRs) displayed that the majority of these DEGs were risk factors for a worse prognosis in HCC patients (Figs. 2b, 2d). These DEGs connected to form a network (Fig. 2c) and among these genes, SLC2A1, SLC38A1, and SLC7A11 were the hub genes (Fig. 2a).
Table 2Univariate Cox regression analysis of 32 ferroptosis-related DEGs for OS of HCC patients
Id | HR | HR.95L | HR.95H | P value |
---|
ABCC1 | 1.373 | 1.156 | 1.630 | <0.001 |
ACACA | 1.613 | 1.229 | 2.117 | <0.001 |
ATG3 | 2.333 | 1.529 | 3.558 | <0.001 |
ATG7 | 2.904 | 1.677 | 5.027 | <0.001 |
AURKA | 1.321 | 1.120 | 1.558 | <0.001 |
CARS1 | 1.685 | 1.243 | 2.284 | <0.001 |
EIF2S1 | 2.462 | 1.607 | 3.773 | <0.001 |
FANCD2 | 1.712 | 1.293 | 2.266 | <0.001 |
G6PD | 1.416 | 1.265 | 1.584 | <0.001 |
GCLM | 1.445 | 1.193 | 1.750 | <0.001 |
HELLS | 1.554 | 1.196 | 2.019 | <0.001 |
HILPDA | 1.435 | 1.240 | 1.661 | <0.001 |
MAFG | 1.760 | 1.385 | 2,239 | <0.001 |
MT3 | 1.426 | 1.164 | 1.746 | <0.001 |
MYB | 3.101 | 1.618 | 5.941 | <0.001 |
NCF2 | 1.314 | 1.128 | 1.531 | <0.001 |
NRAS | 1.799 | 1.377 | 2.350 | <0.001 |
PCBP2 | 2.179 | 1.538 | 3.086 | <0.001 |
PGD | 1.502 | 1.247 | 1.810 | <0.001 |
PRDX1 | 1.659 | 1.303 | 2.112 | <0.001 |
RRM2 | 1.417 | 1.206 | 1.666 | <0.001 |
SLC1A4 | 1.456 | 1.170 | 1.810 | <0.001 |
SLC1A5 | 1.364 | 1.215 | 1.531 | <0.001 |
SLC2A1 | 1.533 | 1.313 | 1.789 | <0.001 |
SLC38A1 | 1.333 | 1.161 | 1.530 | <0.001 |
SLC7A11 | 1.466 | 1.233 | 1.744 | <0.001 |
SQSTM1 | 1.381 | 1.175 | 1.623 | <0.001 |
SRXN1 | 1.697 | 1.343 | 2.142 | <0.001 |
STMN1 | 1.426 | 1.220 | 1.666 | <0.001 |
TXNRD1 | 1.361 | 1.170 | 1.583 | <0.001 |
ZFP69B | 4.076 | 2.486 | 6.684 | <0.001 |
Construction of a prognostic model with good performance in the TCGA cohort
A risk prognosis model was constructed based on 10 FRGs with a risk for worse OS of HCC patients by LASSO regression analysis. To avoid collinearity, the prognostic DEGs were entered into a LASSO regression. Furthermore, a 10-gene signature was generated that was significantly associated with the prognosis with an optimal value of λ (Fig. 3d, e). The coefficient of these 10 genes is exhibited in Table 3. A prognostic model was generated to evaluate the prognosis of each HCC patient as follow:
Risk score = ACACA * 0.004266 + G6PD * 0.071354 + MYB * 0.117039 + SLC2A1 * 0.109265 + SLC38A1 * 0.031264 + SLC7A11 * 0.121545 + SQSTM1 * 0.020266 + SRXN1 * 0.227654 + STMN1 * 0.126750 + ZFP69B * 0.496697.
Table 3LASSO-penalized Cox analysis identifies 10 genes to build a prognostic model
Gene | Coef |
---|
ACACA | 0.004266 |
G6PD | 0.071354 |
MYB | 0.117039 |
SLC2A1 | 0.109265 |
SLC38A1 | 0.031264 |
SLC7A11 | 0.121545 |
SQSTM1 | 0.020266 |
SRXN1 | 0.227654 |
STMN1 | 0.126750 |
ZFP69B | 0.496697 |
The K-M curve exhibited worse survival in the high-risk group of HCC patients (p < 0.001, Fig. 4a). The condition of risk score and survival status in the TCGA cohort is shown in Figure 4c–d. The PCA and t-SNE of the TCGA cohort separated both groups (Fig. 4e–f). Then, the predictive ability of the risk score was estimated for the time-dependent ROC curve in the TCGA cohort, and the areas under the curve (AUCs) reached 0.790, 0.705, and 0.678 for the 1-, 2-, and 3-year OS in this population, respectively (Fig. 4b). Therefore, this risk model had good reliability and specificity.
Validation of the prognostic model in the ICGC cohort
To validate the value of the prediction model from the TCGA database, 260 HCC samples from ICGC were used to validate the ferroptosis-related prognostic model. These HCC patients were stratified into two high and low-risk groups (Fig. 5c). The patients in the high-risk group had significantly worse OS than those in the low-risk group (Fig. 5a). In addition, the ROC curve demonstrated that the AUCs of HCC patients for 1-, 2-, and 3-year OS rates were 0.704, 0.695, and 0.684 (Fig. 5b). Compared with the TCGA cohort, the PCA and t-SNE analysis indicated that patients in both subgroups were different (Fig. 5e–f). In addition, the OS was significantly lower, and the mortality rate was higher in the high-risk group than those in the low-risk group of HCC patients (Fig. 5d).
Functional enrichment analysis of the ferroptosis-related signature
To further investigate the potential functions and pathways related to risk prognosis, GO and KEGG analyses revealed that the DEGs associated with HCC prognosis were enriched in extracellular matrix receptor (ECM-receptor) interaction, rheumatoid arthritis, phagosome, human T cell leukemia virus 1 infection, and others (Fig. 6b–d). Of note, the DEGs in the TCGA cohorts were significantly enriched in many immune-related pathways, such as humoral immune response, immunoglobulin complex, and humoral immune response (adjusted p < 0.05, Fig. 6a). In addition, GO pathway analysis indicated that the collagen-containing and integrin-binding were enriched in both cohorts (adjusted p < 0.05, Fig. 6a, c). The GO circle graph from the TCGA cohort suggested that the biological functions were related to the FRGs (Fig. 3b, c).
Correlation estimation of the clinical characteristics
Then, the interrelationship between the 10 genes and clinical characteristics (e.g., survival status, survival time, gender, age, tumor Gleason, grade, T, M and N stage) were assessed. We showed the expression difference of high-risk and low-risk groups in different stages (T, M and N) through heat map. In addition, the risk scores for males and HCC patients at >65 years old were elevated (Fig. 3a).
Independent prognostic analysis of the ten genes signature
To verify whether the prognostic model was independent of other clinical parameters in HCC prognosis, whether the risk score was an independent prognostic parameter for OS was assessed. Univariate and multivariate Cox regression analyses revealed that high-risk scores were significantly associated with worse OS with a higher HR (HR = 3.606) in the TCGA (Fig. 7a) and high-risk score and HCC stage were independent risk factors for worse OS (risk score: HR = 3.294, 95% confidence interval (CI) = 2.294–4.730) and stage (HR = 2.123, 95% CI = 1.453–3.104) (Fig. 7b). Similar findings were achieved in the ICGC cohort (HR = 3.632, 95% CI = 1.627–8.106; Fig. 7c, d).
Analysis of immunological function
To further investigate the relationship between the HCC and prognostic and immune status, the enrichment scores for different immune cell subpopulations were adopted, which related to immune cells and immunity-related functions with ssGSEA. The frequency of antigen-presenting autologous dendritic cells (aDCs), macrophages (mø), natural killer cells (NK cells), helper T cells 2 (Th2 cells), and regulatory T cells (Tregs) differed significantly between the low and high-risk groups of HCC patients in the TGCA and ICGC cohorts (adjusted p < 0.05, Fig. 8a, c). Furthermore, there was a significant difference in the levels of type II interferon (IFN) responses, major histocompatibility complex (MHC I), human leukocyte antigen (HLA), checkpoint, chemokine receptor (CCR) expression and antigen-presenting cell (APC) costimulation between the high and low-risk groups of HCC patients in the ICGC database (adjusted p < 0.05, Fig. 8b, d).
Construction and evaluation of nomogram
To strengthen the model’s predictive power, a comprehensive prognostic nomogram was generated based on the patients’ risk score and clinical characteristics by integrating the metabolic risk signature, age, stage, grade, and gender (Fig. 9a). The nomogram had an excellent agreement with the ideal curve at 1-, 3-, and 5-years survival rates in the TCGA cohort (Fig. 9b–d).
Discussion
Currently, the incidence and mortality rates of HCC remain high, although there has been improvement in the treatment of HCC.29 Recently, immunotherapy has become a strategy for HCC patients, in addition to surgical resection and target therapies. However, the rate of HCC patients’ response to immune checkpoint inhibitors remains low. Therefore, the discovery of reliable biomarkers for evaluating the prognosis of HCC patients is urgently needed.30,31 Therefore, it is crucial to understand the pathological mechanisms of HCC development and progression.19,32
In this study, 291 DEGs were identified related to ferroptosis and HCC patient’s OS and generated a prognostic risk model with 10 DEGs after evaluating their association with the prognosis of HCC. These 10 DEGs were associated with an increased risk for worse OS in HCC patients, suggesting that they might contribute to the development and progression of HCC. Of these, SLC7A11 was of importance for the pathogenesis of HCC. The prognostic risk model was valuable to predict OS in TCGA and ICGC cohorts, because the high-risk score and HCC stage were independent risk factors for worse OS in HCC patients.33,34 Then, a nomogram was developed by combining the risk scores and clinical features. In addition, the nomogram effectively predicted the 1-, 3-, 5-year OS of HCC patients.27 Therefore, this nomogram might be valuable for the prognosis of HCC patients in the clinic.
Of interest, the 10 ferroptosis-related DEGs were associated functionally with tumorigenesis and development because they were significantly enriched in many immune-related pathways, such as lymphocyte-mediated immunity, humoral immune response, and immunoglobulin complex.24 This suggested that ferroptosis might be closely related to antitumor immunity, which agreed with a previous observation that antigen presentation-related genes were differentially expressed between the high and low-risk groups.35 It is well-known that the glutamate/cystine antiporter system Xc- is the most important target for erastin during erastin-induced ferroptosis.36 Xc- can transport cysteine (the major form of cysteine in the atmosphere) into cells.37 The inhibition of the Xc- system reduces cellular cysteine, which means that it is unavailable for the synthesis of GSH (glutathione, r-glutamyl cysteingl + glycine), which is an important antioxidant. GSH deficiency in proteins or membranes can lead to an accumulation of reactive oxygen species (ROS) and subsequent ferroptotic cell death. The combination of erastin into SLC7A5 interferes with cystine uptake by the SLC3A2/SLC7A11 complex in the trans form. The ferroptosis-related DEGs in HCC were highly enriched in the phagocytosis pathway, which suggested a close relationship between tumor induction and phagocytosis.
The prognostic model was composed of 10 FRGs (ACACA, SLC7A11, SQSTM1, SRXN1, SLC2A1, SLC38A1, MYB, ZFP69B, G6PD, STMN1). Significantly, these 10 genes were associated with a worse OS in HCC patients. ACACA is the rate-limiting enzyme of fatty acid synthetase (FAS). Inhibition of ACACA expression inhibits proliferation and induces apoptosis of prostate cancer LNCaP cells.38 Previous studies have shown that silencing ACACA decreases prostate cancer cell proliferation by affecting the NAD+/NADH balance. SLC7A11 is a vital subunit in the Xc-system and is crucial for ferroptosis and oxidation resistance.39 Furthermore, IFN can downregulate SLC7A11 expression in tumor cells and a combination of IFN and radiotherapy can decrease cystine transport and antioxidant storage in tumor cells, which leads to ferroptosis.40 A previous study identified different p62-positive structures in patients with HCC, steatohepatitis, and alcoholic hepatitis, such as Mallory-Denk bodies and intracellular hyaline bodies.41SRXN1 plays a vital regulatory role in the antioxidant response in eukaryotic cells. High expression of SRXN1 might lead to a poor prognosis in HCC, further reducing the overall survival rate of HCC patients.42 However, the role of SRXN1 in HCC remains unclear. SLC2A1 is a member of the solute carrier family 2 and upregulated SLC2A1 expression is significantly associated with invasiveness and late clinical stage of gastric cancer.43 High SLC38A1 expression is associated with worse survival and was an independent prognostic factor for a poor prognosis in gastric cancer.44MYB and ZFP69B are highly expressed in most tumor tissues. However, the function of MYB in HCC requires further research.45 Glucose-6-phosphate dehydrogenase (G6PD) is an important metabolic enzyme in glycolysis, and it can regulate the proliferation and apoptosis of HCC.46 Increased G6PD activity promotes tumor cell growth and proliferation.47 The activity of G6PD is closely related to inflammation, infection, and the progression of the tumor and G6PD might be a therapeutic target for cancer therapies. STMN1 enhances the crosstalk between HCC and hepatic stellate cells by promoting the mesenchymal-o-epithelial transition (MET) process of HCC cells.48 Therefore, these DEGs are crucial for the development and progression of HCC and are potential therapeutic targets for interventions in HCC.
This study has several limitations. First, the prognostic model was established and validated on public databases, without validation in new samples and prospective experimental trials. In addition, the relationship between immune and risk scores was not experimentally confirmed for the prognosis of HCC.
Future directions
Ferroptosis is a new type of cell death with unique properties and recognition functions. Cell ferroptotic death could be significant for the treatment of cancer in the future. This study screened biomarkers related to the prognosis of liver cancer based on ferroptosis-related regulatory genes, to provide help for the clinical treatment of liver cancer in the future. In addition, the prognosis of liver cancer based on ferroptosis-related non-coding RNA will be studied in the future, to explore the regulatory relationship between ferroptosis-related non-coding RNA in the tumor immune microenvironment of liver cancer, and antitumor drugs will be screened to provide help for immunotherapy in liver cancer.
Conclusions
In total, 10 FRGs were identified that were associated with a worse OS in HCC patients. This generated model and nomogram effectively predicted the OS of HCC patients in both databases and might be valuable for the prognosis of HCC in the clinic. In the future, the potential mechanisms that underlie ferroptosis in antitumor immunity against HCC will be investigated.
Abbreviations
- AUC:
area under the curve
- CI:
confidence interval
- DEGs:
differentially expressed genes
- FDR:
false discovery rate
- FRGs:
ferroptosis-related genes
- GO:
gene ontology
- HCC:
hepatocellular carcinoma
- HR:
hazard ratio
- ICGC:
International Cancer Genome Consortium
- KEGG:
Kyoto Encyclopedia of Genes and Genomes
- K-M:
Kaplan-Meier
- LASSO:
least absolute shrinkage and selection operator
- LIHC:
liver hepatocellular carcinoma
- OS:
overall survival
- PCA:
principal component analysis
- ROC:
receiver operating characteristic
- ssGSEA:
single-sample gene set enrichment analysis
- TCGA:
The Cancer Genome Atlas
- t-SNE:
t-distributed stochastic neighbor embedding
Declarations
Data sharing statement
The data used to support the findings of this study are included in the article.
Funding
This project was supported by the grants from Administration of Traditional Chinese Medicine of Guangdong Province (20201180; 20211223); Science and Technology Special Project of Zhanjiang (2019A01009); Basic and Applied Basic Research Program of Guangdong Province (2019A1515110201); Program of Department of Natural Resources of Guangdong Province (No. GDNRC [2020]038 and [2021]53); Discipline Construction Project of Guangdong Medical University (4SG21004G).
Conflict of interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Authors’ contributions
Data curation (XSH and ZZ), Funding acquisition (LXL), Investigation (XSH and XLL), Project administration (LXL), Writing – original draft (LXL), Writing – review & editing (ZZ, HL and LXL).