Abstract Recent research has highlighted the therapeutic potential of citrus-derived dietary 5,6,7,4′-tetramethoxyflavone (TMF) against HeLa cancer. Our study aims to elucidate its mechanisms of action through proteomics analysis, network pharmacology, and molecular docking. The results suggested that TMF demonstrated efficacy by upregulating CD40, CD40L, Fas, Fas-L, HSP27, HSP60, IGFBP-1, IGFBP-2, IGF-1sR, Livin, p21, p27, sTNFR2, TRAILR2, TRAILAR3, TRAILR4, XIAP, p-Sre, p-Stat1, p-Stat2 p-c-Fos, p-SMAD1, p-SMAD2, p-SMAD4, p-SMAD5, p-IκBα, p-MSK1, p-NFκB, p-TAK1, p-TBK1, p-ZAP70, and p-MSK2, while downregulating p-EGFR, p-ATF2, p-cJUN, p-HSP27, p-JNK, and p-GSK3A. These targets are primarily involved in MAPK, apoptosis, and TNF signaling pathways. Notably, p21, p27, EGFR, SMAD4, JNK, ATF2, and c-JUN merged as pivotal targets contributing to TMF's anti-cancer efficacy against HeLa cells. This study is first to delineate the potential signaling pathways and core targets of TMF in treating of HeLa cancer, paving the way for further exploration of TMF's medical potential. Keywords: Polymethoxyflavones; 5,6,7,4′-tetramethoxyflavone; Tetrametlglscutellarein; Proteomics; Mechanisms; Network pharmacology; Molecular docking; Cervical cancer; HeLa 1. Introduction Cervical cancer (CCA) ranks as the fourth most common malignancy among women globally, with approximately 604,000 new cases and 342,000 deaths reported in 2020 [[29]1]. The vast majority of CCA cases are linked to persistent infections with high-risk human papillomavirus [[30]2]. Due to the absence of organized screening and vaccination programs, developing countries account for nearly 90 % of the global burden of new cases and deaths [[31]3]. Chemotherapy remains a cornerstone in the treatment of CCA, particularly in recurrent and advanced metastatic case, with cisplatin, gemcitabine, and paclitaxel commonly employed to improve patient prognosis and overall survival. However, the severe toxicities associated with these first-line chemotherapeutics, including neutropenia, vomiting, diarrhea, and even fatal outcomes, often undermine their clinical benefits [[32]4]. Therefore, there is a critical need for research into effective CCA chemotherapeutics agents, with reduced toxicity. Nature products have long been a significant source of anticancer agents [[33]5]. An analysis of small-molecule antitumor drugs approved globally between the 1940s and 2011 worldwide revealed that approximately 50 % were natural products or their derivatives [[34]6]. Among these, 5,6,7,4′-tetramethoxyflavone (TMF), a dietary polymethoxyflavone (PMF), is derived from various plant species, including Citrus species [[35]7], Nervilia concolor [[36]8], Marrubium vulgare [[37]9], Cissus assamica [[38]10] and Clerodendranthus spicatus [[39]11], as well as the traditional Chinese herbal formula Gehua Jiecheng decoction [[40]12]. Notably, TMF is abundantly found in the peel of citrus fruits, such as mandarin, orange, lemon and grapefruit. In different varieties of loose-skin mandarins, TMF content varies significantly, ranging from 0 to 13.7 mg/g, with the highest levels in 'Cupigoushigan' tangerine, while in sweet oranges, it generally ranges from 1 to 2 mg/g [[41]7] ([42]Fig. 1). TMF has been reported to possess broad-spectrum antitumor activity. For instance, Manthey et al. [[43]13]isolated TMF from orange peel and found that it significantly inhibited the proliferation of various human cancer cell lines, including DMS-114, HT-29, MCF-7, DU-145, and MDA-MB-435. Subsequently, Ohtani et al. [[44]14]demonstrated that TMF significantly enhanced the antitumor activity of vincristine in K562/ADM cancer cells by increasing vincristine uptake via P-glycoprotein inhibition. Additionally, TMF was shown to improve the inhibitory activity of 5-fluorouracil against HT29 colon cancer cells by upregulating CDH1 and downregulating ZEB1 and SNAI1 [[45]15]. However, the current research on TMF's antitumor activity is largely confined to a limited number of in vitro studies. One of the main challenges has been the lower concentration of TMF compared to other well-studied PMFs, such as nobiletin and tangeretin. Our previous research [[46]16] synthesized 12 PMFs with high purity at gram scale through the methylation of polyhydroxyflavones and assessed their pharmacokinetic parameters and antitumor activity against HeLa, A549, HepG2, and HCT116 cancer cells. Among these PMFs, only TMF significantly inhibited HeLa cells and was detectable in rat plasma, suggesting its potential as a promising therapeutic agent for CCA. Nevertheless, the underlying mechanisms of TMF's action against HeLa cancer remain unclear. Fig. 1. [47]Fig. 1 [48]Open in a new tab The chemical structure of 5,6,7,4′-tetramethoxyflavone (TMF) and its primary sources. Antibody microarray, a cutting-edge proteomic profiling technology, allows for the rapid and simultaneous detection of multiple target protein expression levels. Compared to traditional western blotting, this technique offers advantages such as miniaturization, high throughput, enhanced sensitivity, and reduced sample requirement, making it a valuable tool in assessing tumor-related protein expression [[49]17]. For instance, a previous study [[50]18] utilized an antibody microarray to investigate the mechanisms behind baicalin-induced apoptosis in gastric cancer SGC-7901 cells. The study found that baicalin's apoptotic effect were primarily mediated through the upregulation of BAD, caspase-3, caspase-8, FasL, HSP60, TRAIL R4, IGF-I, IGF-II, P21, P27, p53, and SMAC, along with a marked downregulation of XIAP. These findings were subsequently confirmed by western blotting, which yielded results consistent with those of the antibody microarray. This study aims to elucidate the molecular mechanisms underlying TMF's action against HeLa cancer by integrating proteomics analysis with network pharmacology and molecular docking approaches. 2. Materials and method 2.1. Chemicals and reagents The TMF (HPLC purity >98 %) used in this study was synthesized according to our previously published method [[51]16]. Dimethyl sulfoxide (DMSO) (biochemical grade, BCCD7238) was obtained from Sigma-Aldrich (Shanghai) Trading Co., LTD. Fetal bovine serum (20200408) and trypsin-EDTA (70070200) solution were sourced from Biosharp (Beijing, China), while the penicillin-streptomycin (J190033) solution was procured from HyClone Laboratories (UT, USA). 2.2. Cell cultures and treatments HeLa cells (ATCC Cat# CCL-2, RRID: CVCL_0300) were acquired from Chengdu University of Traditional Chinese Medicine. The cells were seeded in 5 mL of DMEM containing 10 % fetal bovine serum, penicillin (100 IU/mL), and streptomycin (100 μg/mL). Cultures were maintained in a humidified incubator at 37 °C with 5 % CO[2]. TMF was dissolved in DMSO and prepared at a concentration of 50 mM. The cells were then passaged into six 60 mm Petri dishes at a density of 0.5 × 10^6 cells/well, with three dishes allocated to the experimental group and three to the control group. After 24 h of incubation, when the cell density reached a density of 1 × 10^6 cells/well, the experimental group was treated with 5 μL of TMF, resulting in a final concentration of 50 μM [[52]16], while the control group received 5 μL of DMSO, yielding a final concentration of 0.1 % (v/v). Following an additional 24 h of incubation, the cells were washed twice with PBS and harvested into cryopreservation tubes. 2.3. Antibody array analysis To explore the signaling pathways triggered by TMF, apoptosis-related proteins and phosphorylation sites were analyzed using the Human Phosphorylation Pathway Profiling Array C55 (Raybiotech, Inc., USA, AAP-PPP-1), encompassing 55 proteins, including EGFR (Ser1070), JAK1 (Tr1-22), JAK2 (Tyr1007/1008), SHP1 (Ser591), SHP2 (Tyr542), Src (Tyr419), Stat1 (Ser727), Stat2 (Tyr689), Stat3 (Tyr705), Stat5 (Tyr694), Stat6 (Tyr641), TYK2 (Tyr1054), ATF2 (T69/71), C-Fos (T232), c-JUN (S73), SMAD1 (S463/465), SMAD2 (S245/250/255), SMAD4 (T277), SMAD5 (S463/465), ATM (S1981), eIF2a (S51),HDAC2 (S394), HDAC4 (S632), IκBα (S32), MSK1 (S376), NFκB (S536), TAK1 (S412), TBK1 (S172), ZAP70(Y292), CREB (S133), ERK1 (T202/Y204)/ERK2 (Y185/Y187), HSP27 (S82), JNK (T183), MEK (S217/221), MKK3 (S189), MKK6 (S207), MSK2 (S360), p38 (T180/Y182), p53 (S15), RSK1 (S380), RSK2 (S382), AKT (S473), AMPKa (T172), BAD (S112), 4E-BP1 (T36),GSK3a (S21), GSK3b (S9), mTOR (S2448), P27 (T198), P70S6K (T421/S424), PDK1 (S241), PRAS40 (T246), PTEN (S380), RAF-1 (S301), and PRS6 (S235/236). Additionally, the Human Apoptosis Array C1 (Raybiotech, Inc., USA, AAP-APO-1) was used, covering 43 proteins, including BAD, BAX, BCL-2, BCL-w, BID, BIM, Caspase 3, Caspase 8, CD40, CD40 Ligand, cIAP-2, Cytochrome-C, DR6, Fas Ligand, HSP27, HSP60, HSP70, HTRA, IGF-1, IGF-2, IGFBP-1, IGFBP-2, IGFBP-3, IGFBP-4, IGFBP-5, IGFBP-6, IGF-1sR, Livin, p21, p27, p53, SMAC, Survivin, sTNFRI, sTNFR2, TNF-α, TNF-β, TRAIL R1, TRAIL R2, TRAIL R3, TRAIL R4, and XIAP. The antibody experiment followed the kit's instructions meticulously. Total protein was extracted from HeLa cells using a cell lysis buffer, and protein concentrations were measured using a BCA Protein Assay Kit (Thermo Scientific, Pierce™ 23227). The antibody array membrane was placed in an incubator box, followed by addition of 2 mL of blocking buffer. After a 1-h incubation at room temperature, the blocking buffer was removed, and 350 μg of protein (2 mL) extracted from these HeLa cells was loaded onto the membrane. Following a 12-h incubation at 4 °C, the sample was removed, and the membrane was washed 3–4 times with 2 mL of Wash Buffer I-II for 5 min per wash. Subsequently, 1 mL of Biotin-Conjugated Anti-Cytokine was added and incubated at room temperature for 2 h, followed by another round of washing with Wash Buffer I-II. Lastly, 2 mL of HRP-Streptavidin Concentrate was applied. After completing these steps, Detection Buffer C/D (0.25/0.25 mL) was added, and protein visualization was performed using the ImageQuant LAS4000 Scanner (General Electric Company, USA). For significance analysis, fold change (FC) was used as the basic statistic. The formula Ln FC = |E − C| was used, where E represents the average expression level in the experimental group and C represents the control group. Differentially expressed proteins (DEPs) were identified with thresholds of FC ≥ 1.2 [[53]19]. 2.4. Gene ontology (GO) functional annotation GO analysis is a widely utilized system for gene function classification in biology, providing descriptions of gene products across three categories: biological processes (BP), cellular components (CC), and molecular functions (MF) [[54]20]. In this study, GO functional annotation was applied to the 37 DEPs identified via antibody array analysis using the online DAVID database ([55]https://david.ncifcrf.gov/). Terms were ranked and selected based on an adjusted p-value threshold of <0.05, with the p-value, DEP count, and enrichment factor serving as metrics to assess the significance and relevance of the enriched terms. Visualization of these terms were conducted using the free online platform Weishengxin ([56]https://www.bioinformatics.com.cn [last accessed on 30 July 2022]). 2.5. Kyoto Encyclopedia of genes and Genomes (KEGG) pathway enrichment For KEGG pathway enrichment analysis of the 37 DEPs, Fisher's exact test was employed through the data package clusterProfiler from R/Bioconductor, with Homo sapiens as the background reference. Pathway terms were selected if the number of DEPs enriched in a given pathway exceeded five and the p-value was less than 0.05, with subsequent visualization also performed via Weishengxin ([57]https://www.bioinformatics.com.cn [last accessed on 30 July 2022]). 2.6. Protein-protein interaction (PPI) analysis To identify key DEPs, protein-protein interaction (PPI) analysis was conducted using the STRING online platform, version 11.5 ([58]https://cn.string-db.org/), with the analysis restricted to Homo sapiens. The significance of DEPs within the network was determined by calculating the number of node connections using R software (version 3.6.0), reflecting their importance within the PPI network. 2.7. Molecular docking and molecular dynamics (MD) simulation Molecular docking analysis was conducted to identify the binding sites and calculate the binding energy (BE) between TMF and the DEPs using AutoDock Vina 1.1.2. The 2D and 3D chemical structures of TMF were created with ChemDraw Professional 19.0. Before performing the docking analysis, the 3D structure of TMF was optimized using MM2 minimization. The 3D structures of the DEPs were obtained from the Protein Data Bank (PDB) ([59]https://www.rcsb.org/) and the AlphaFold Protein Structure Database ([60]https://alphafold.ebi.ac.uk/). BE serves as an indicator of bond strength, with BE < −5 kcal/mol suggesting a stable interaction between TMF and the DEP, and lower BE values reflecting stronger binding affinities [[61]21]. Receptor-ligand interactions were visualized using PyMol 2.5.2 [[62]22]. To validate the docking results, MD simulations were performed using Gromacs 2023. CHARMM 36 force field parameters were applied for the protein [[63]23], while the ligand topology was constructed using GAFF2 force field parameters. The TIP3P water model was employed to solvate the system. Electrostatic interactions were handled using the particle mesh Ewald method and the Verlet algorithm. MD simulations were executed under conditions of constant temperature (300 K) and pressure (1 bar) over 2,500,000 steps, with a step size of 2 fs and a total simulation time of 50 ns. Root mean squared deviation (RMSD) were calculated to assess the stability of complex compounds. 3. Results 3.1. Effect of TMF on protein expression level in HeLa cells As shown in [64]Fig. 2 and [65]Table 1, treatment of HeLa cells with 50 μM of TMF for 24 h led to the identification of 37 DEPs. Specifically, 17 apoptosis-related proteins were upregulated, including CD40, CD40L, Fas, Fas-L, HSP27, HSP60, IGFBP-1, IGFBP-2, IGF-1sR, Livin, p21, p27, sTNFR2, TRAILR2, TRAILAR3, TRAILR4, and XIAP ([66]Fig. 2A). In the JAK/STAT signaling pathway, EGFR (Ser1070) was downregulated, while Sre (Tyr419), Stat1 (ser727), and Stat2 (Tyr689) were upregulated ([67]Fig. 2B). In the TGFb signaling pathway, ATF2 and c-JUN were downregulated, whereas c-Fos, SMAD1, SMAD2, SMAD4, and SMAD5 were upregulated ([68]Fig. 2C). Within the NFκB signaling pathway, IκBα (S32), MSK1 (S376), NFκB (S536), TAK1 (S412), TBK1 (S172), and ZAP70 (Y292) were upregulated ([69]Fig. 2D). In the MAPK signaling pathway, HSP27 and JNK were downregulated, and MSK2 was upregulated ([70]Fig. 2E), while in the AKT signaling pathway, GSK3A was downregulated ([71]Fig. 2F). Fig. 2. [72]Fig. 2 [73]Open in a new tab Effects of TMF on the proteome of HeLa cells using Human cytokine antibody microarray analysis. (A) Apoptosis-associated antibody array, (B) JAK/STAT signaling pathway array, (C) TGF-β signaling pathway array, (D) NFκB signaling pathway array, (E) MAPK signaling pathway array, (F) AKT signaling pathway array, and (G) Graphic representation of the 37 differentially expressed proteins (DEPs) in HeLa cells treated with TMF (50 μM) (FC > 1.2). Among the 37 DEPs, 12 were downregulated, and 25 were upregulated. Note: Red columns and frames indicate significant upregulation, while blue columns and frames indicate significant downregulation. Table 1. Thirty-seven differentially expressed proteins in HeLa cells treated with TMF (50 μM) for 24 h. Protein AveExp.E AveExp.C logFC FC (>1.2) regulation Entrez ID Uniport ID Fas 15.75511 15.43281 0.322298 1.25032 up 355 [74]P25445 CD40 15.56755 15.20546 0.362094 1.28529 up 958 [75]P25942 HSP60 15.20289 14.85137 0.351527 1.275911 up 3329 [76]P10809 CD40L 14.92035 14.57327 0.347079 1.271982 up 959 [77]P29965 p21 14.89228 14.57968 0.312606 1.241949 up 1026 [78]P38936 IGF-1sR 14.44195 13.92932 0.512632 1.426651 up 3479 [79]P05019 IGFBP-2 14.41411 14.0321 0.382004 1.303151 up 3485 [80]P18065 FasL 14.20255 13.89918 0.303369 1.234022 up 356 [81]P48023 p27 14.09551 13.80288 0.292625 1.224867 up 1027 [82]P46527 TRAILR-2 14.03571 13.6818 0.353904 1.278014 up 8795 [83]O14763 sTNF-R2 13.98099 13.65722 0.323763 1.251591 up 7133 [84]P20333 IGFBP-1 13.88088 13.39709 0.483791 1.398413 up 3484 [85]P08833 Livin 13.75238 13.2091 0.543285 1.457287 up 79444 [86]Q96CA5 TRAILR-4 13.67824 13.24379 0.434446 1.351391 up 8793 [87]Q9UBN6 XIAP 13.62999 13.35344 0.276551 1.211295 up 331 [88]P98170 TRAILR-3 13.62021 13.34216 0.278043 1.212549 up 8794 [89]O14798 HSP27 13.53669 13.22952 0.307176 1.237283 up 3315 [90]P04792 IκBα 8.35116 9.186733 0.835573 1.784566 down 4792 [91]P25963 MSK1 8.871351 9.655978 0.784627 1.722647 down 9252 [92]O75582 NFκB 12.08366 12.80398 0.720321 1.647548 down 5970 [93]Q04206 c-Jun 9.157372 9.781032 0.62366 1.540779 down 3725 [94]P05412 ZAP70 11.69695 12.20304 0.506088 1.420194 down 7535 [95]P43403 GSk3a 12.46839 12.94065 0.472255 1.387276 down 2931 [96]P49840 TBK1 10.05922 10.50457 0.445348 1.361643 down 29110 [97]Q9UHD2 ATF2 8.379682 8.81314 0.433458 1.350467 down 1386 [98]P15336 JNK 9.685204 10.07231 0.387102 1.307764 down 5599 [99]P45983 EGFR 8.381586 8.761551 0.379965 1.30131 down 1956 [100]P00535 HSP27 9.799492 10.09028 0.290792 1.223312 down 3315 [101]P04792 TAK1 10.17156 10.45558 0.284019 1.217582 down 6885 [102]O43318 c-Fos 9.039687 8.770829 0.268858 1.204853 up 2353 [103]P01100 SMAD2 7.826231 7.474112 0.352119 1.276434 up 4087 [104]P15796 SRC 10.07318 9.584023 0.489155 1.403622 up 6714 [105]P12931 STAT1 11.23691 10.74147 0.495444 1.409755 up 6772 [106]P42224 SMAD1 8.438251 7.887525 0.550726 1.464823 up 4086 [107]Q15797 SMAD5 10.96276 10.40652 0.556235 1.470427 up 4090 [108]Q99717 STAT2 11.63225 11.07146 0.560783 1.47507 up 6773 [109]P52630 SMAD4 11.34809 10.64773 0.700358 1.624908 up 4089 [110]Q13485 [111]Open in a new tab E, experiment; C, control. 3.2. GO functional annotation analysis A total of 924 BP terms were enriched, including those related to the transmembrane receptor protein serine/threonine kinase signaling pathway, transforming growth factor beta receptor signaling pathway, SMAD protein signal transduction, response to transforming growth factor beta, response to muscle stretch, and response to mechanical stimulus. Notably, the response to mechanical stimulus had the highest number of enriched genes, with 12 genes identified, including NFKBIA, RELA, JUN, MAPK8, EGFR, FOS, FAS, TNFRSF10B, CD40, IGFBP2, SRC, and STAT1. The necroptotic signaling pathway and regulation of cardiac muscle tissue regeneration both exhibited the highest enrichment factor (EF = 172.87), each associated with two genes: FASLG/FAS and CDKN1B/CDKN1A, respectively. Regarding CC, seven GO terms were enriched, including the SMAD protein complex, transcription regulator complex, membrane raft, membrane microdomain, membrane region, RNA polymerase II transcription regulator complex, and nuclear chromatin. Among these, the transcription regulator complex encompassed the highest number of genes, while the SMAD protein complex exhibited the highest enrichment degree (EF = 312.97). For MF, 57 GO terms were identified, such as ubiquitin-like protein ligase binding, ubiquitin protein ligase binding, TRAIL binding, R-SMAD binding, DNA-binding transcription activator activity, RNA polymerase II-specific, DNA-binding transcription activator activity, and DNA-binding transcription factor binding. Ubiquitin-like protein ligase binding was the most significant term (p = 1.69E-14) and included the largest number of genes (13 genes: NFKBIA, RELA, JUN, EGFR, CDKN1A, TNFRSF1B, HSPD1, SMAD2, CD40, SRC, STAT1, SMAD5, STAT2), while TRAIL binding had the highest enrichment degree (EF = 294.93). The top 20 most significant terms, based on adjusted p-value, were visualized and are presented in [112]Fig. 3 and [113]Table 2. Fig. 3. [114]Fig. 3 [115]Open in a new tab GO functional annotation for the 37 DEPs, encompassing biological processes, cellular components, and molecular functions. Table 2. Top 20 most significant GO terms enriched based on the 37 differentially expressed proteins. GO ID Description P. adjust Gene ID No Enrich factor Biological processes (BP) 0007178 transmembrane receptor protein serine/threonine kinase signaling pathway 2.67E-06 JUN/MAP3K7/FOS/XIAP/SMAD2/SRC/SMAD1/SMAD5/SMAD4 9 13.37 0007179 transforming growth factor beta receptor signaling pathway 1.05E-06 JUN/MAP3K7/FOS/SMAD2/SRC/SMAD1/SMAD5/SMAD4 8 20.85 0060395 SMAD protein signal transduction 1.05E-06 JUN/FOS/SMAD2/SMAD1/SMAD5/SMAD4 6 44.45 0071559 response to transforming growth factor beta 3.56E-06 JUN/MAP3K7/FOS/SMAD2/SRC/SMAD1/SMAD5/SMAD4 8 16.27 0035994 response to muscle stretch 4.61E-06 NFKBIA/RELA/JUN/FOS 4 115.25 0009612 response to mechanical stimulus 6.24E-12 NFKBIA/RELA/JUN/MAPK8/EGFR/FOS/FAS/TNFRSF10B/CD40/IGFBP2/SRC/STAT1 12 29.63 0050670 regulation of lymphocyte proliferation 1.19E-05 ZAP70/CDKN1A/TNFRSF1B/CD40LG/CD40/IGFBP2/IGF1 7 17.45 0043122 regulation of I-kappaB kinase/NF-kappaB signaling 2.56E-06 RELA/TBK1/HSPB1/MAP3K7/FASLG/TNFRSF10B/CD40/STAT1 8 17.51 0051090 regulation of DNA-binding transcription factor activity 1.38E-06 NFKBIA/RPS6KA5/RELA/JUN/ATF2/MAPK8/MAP3K7/FOS/CD40LG/CD40 10 12.00 2001233 regulation of apoptotic signaling pathway 6.54E-06 RELA/GSK3A/MAPK8/HSPB1/FASLG/FAS/TNFRSF10B/SRC/IGF1 9 11.50 0061614 pri-miRNA transcription by RNA polymerase II 4.07E-06 RELA/JUN/FOS/SMAD1/SMAD4 5 55.17 1901522 positive regulation of transcription from RNA polymerase II promoter involved in cellular response to chemical stimulus 9.24E-06 RELA/SMAD1/SMAD5/SMAD4 4 94.29 0001819 positive regulation of cytokine production 2.46E-06 RELA/TBK1/ATF2/HSPB1/MAP3K7/CD40LG/HSPD1/CD40/SRC/STAT1 10 11.18 2001237 negative regulation of extrinsic apoptotic signaling pathway 5.46E-06 RELA/FASLG/FAS/TNFRSF10B/SRC/IGF1 6 29.92 0032943 mononuclear cell proliferation 7.79E-07 ZAP70/TBK1/CDKN1A/TNFRSF1B/CD40LG/HSPD1/CD40/IGFBP2/IGF1 9 17.03 0046651 lymphocyte proliferation 7.79E-07 ZAP70/TBK1/CDKN1A/TNFRSF1B/CD40LG/HSPD1/CD40/IGFBP2/IGF1 9 17.16 0070661 leukocyte proliferation 1.05E-06 ZAP70/TBK1/CDKN1A/TNFRSF1B/CD40LG/HSPD1/CD40/IGFBP2/IGF1 9 15.66 0007249 I-kappaB kinase/NF-kappaB signaling 7.79E-07 NFKBIA/RELA/TBK1/HSPB1/MAP3K7/FASLG/TNFRSF10B/CD40/STAT1 9 17.35 0097191 extrinsic apoptotic signaling pathway 3.29E-07 RELA/GSK3A/TNFRSF10C/FASLG/FAS/TNFRSF1B/TNFRSF10B/SRC/IGF1 9 20.84 0071560 cellular response to transforming growth factor beta stimulus 3.19E-06 JUN/MAP3K7/FOS/SMAD2/SRC/SMAD1/SMAD5/SMAD4 8 16.66 Cellular components (CC) 0005667 transcription regulator complex 5.14E-04 RELA/JUN/FOS/SMAD2/SMAD1/SMAD5/SMAD4 7 9.24 0071141 SMAD protein complex 3.91E-08 SMAD2/SMAD1/SMAD5/SMAD4 4 312.97 0090575 RNA polymerase II transcription regulator complex 4.63E-03 JUN/FOS/SMAD2/SMAD4 4 13.20 0000790 nuclear chromatin 9.76E-03 RELA/JUN/SMAD2/STAT1/SMAD4 5 7.28 0098589 membrane region 6.21E-04 ZAP70/EGFR/FASLG/FAS/TNFRSF1B/SRC 6 10.02 0045121 membrane raft 6.21E-04 ZAP70/EGFR/FASLG/FAS/TNFRSF1B/SRC 6 10.43 0098857 membrane microdomain 6.21E-04 ZAP70/EGFR/FASLG/FAS/TNFRSF1B/SRC 6 10.40 Molecular functions (MF) 0044389 ubiquitin-like protein ligase binding 2.83E-12 NFKBIA/RELA/JUN/EGFR/CDKN1A/TNFRSF1B/HSPD1/SMAD2/CD40/SRC/STAT1/SMAD5/S TAT2 13 20.75 0031625 ubiquitin protein ligase binding 6.65E-10 NFKBIA/RELA/JUN/EGFR/CDKN1A/TNFRSF1B/HSPD1/SMAD2/CD40/SRC/SMAD5 11 18.65 0032813 tumor necrosis factor receptor superfamily binding 1.04E-03 FASLG/CD40LG/STAT1 3 32.06 0005164 tumor necrosis factor receptor binding 4.42E-04 FASLG/CD40LG/STAT1 3 47.57 0045569 TRAIL binding 4.32E-06 TNFRSF10C/TNFRSF10B/TNFRSF10D 3 294.93 0046332 SMAD binding 1.18E-05 JUN/FOS/SMAD2/SMAD1/SMAD4 5 30.72 0070412 R-SMAD binding 4.86E-06 JUN/FOS/SMAD2/SMAD4 4 85.49 0061629 RNA polymerase II-specific DNA-binding transcription factor binding 1.93E-05 NFKBIA/RELA/JUN/ATF2/FOS/SRC/STAT1 7 12.42 0001102 RNA polymerase II activating transcription factor binding 1.49E-03 JUN/ATF2/FOS 3 27.82 0030291 protein serine/threonine kinase inhibitor activity 4.42E-04 HSPB1/CDKN1B/CDKN1A 3 46.08 0019903 protein phosphatase binding 1.53E-03 TBK1/EGFR/CDKN1B/STAT1 4 14.04 0070878 primary miRNA binding 1.04E-03 SMAD2/SMAD1 2 122.89 0019902 phosphatase binding 4.42E-04 TBK1/EGFR/CDKN1B/SMAD2/STAT1 5 13.29 0070411 I-SMAD binding 1.93E-05 SMAD2/SMAD1/SMAD4 3 134.06 0031995 insulin-like growth factor II binding 1.04E-03 IGFBP2/IGFBP1 2 122.89 0140297 DNA-binding transcription factor binding 9.89E-06 NFKBIA/RELA/JUN/ATF2/FOS/SMAD2/SRC/STAT1 8 11.20 0001228 DNA-binding transcription activator activity, RNA polymerase II-specific 4.86E-06 RELA/JUN/ATF2/FOS/SMAD2/STAT1/SMAD1/STAT2/SMAD4 9 10.08 0001216 DNA-binding transcription activator activity 4.86E-06 RELA/JUN/ATF2/FOS/SMAD2/STAT1/SMAD1/STAT2/SMAD4 9 10.05 0017151 DEAD/H-box RNA helicase binding 1.04E-03 SMAD1/SMAD5 2 122.89 0033613 activating transcription factor binding 1.42E-05 RELA/JUN/ATF2/FOS/SMAD2 5 28.92 [116]Open in a new tab 3.3. KEGG pathway enrichment analysis Based on an adjusted p-value of <0.05, 115 signaling pathways were identified as enriched. The top 20 pathways revealed that the TMF's inhibitory effects on HeLa cells were predominantly associated with the MAPK signaling pathway, TNF signaling pathway, apoptosis, PD-L1 expression and PD-1 checkpoint pathway in cancer, pancreatic cancer, colorectal cancer, NF-kappa B signaling pathway, and FoxO signaling pathway. Of these, the MAPK signaling pathway compassed the highest number of genes, totaling 12 (P = 1.79E-08, EF = 9.93, including RPS6KA5, RELA, JUN, ATF2, MAPK8, EGFR, HSPB1, MAP3K7, FOS, FASLG, FAS, IGF1), while the apoptosis-multiple species pathway had the highest enrichment score (EF = 22.35). These pathways are primarily involved in environmental information processing, cellular processes, organismal systems, and human diseases ([117]Fig. 4B). Pathways with more than six enriched genes were visualized and are presented in [118]Fig. 4A and [119]Table 3. Fig. 4. [120]Fig. 4 [121]Open in a new tab (A) KEGG pathway enrichment analysis of the 37 DEPs, and (B) KEGG classification of these pathways. Table 3. KEGG signaling pathway enrichment for the 37 differentially expressed proteins. ID Signaling pathway p.adjust geneID Count Enrich factor hsa04010 MAPK 1.79E-08 RPS6KA5/RELA/JUN/ATF2/MAPK8/EGFR/HSPB1/MAP3K7/FOS/FASLG/FAS/IGF1 12 9.73 hsa04668 TNF 5.94E-10 NFKBIA/RPS6KA5/RELA/JUN/ATF2/MAPK8/MAP3K7/FOS/FAS/TNFRSF1B 10 21.28 hsa04620 Toll-like receptor 6.75E-09 NFKBIA/RELA/JUN/TBK1/MAPK8/MAP3K7/FOS/CD40/STAT1 9 20.63 hsa04659 Th17 cell differentiation 7.64E-09 NFKBIA/RELA/JUN/ZAP70/MAPK8/FOS/SMAD2/STAT1/SMAD4 9 20.05 hsa04210 Apoptosis 4.01E-08 NFKBIA/RELA/JUN/MAPK8/FOS/XIAP/FASLG/FAS/TNFRSF10B 9 15.78 hsa04621 NOD-like receptor 2.82E-07 NFKBIA/RELA/JUN/TBK1/MAPK8/MAP3K7/XIAP/STAT1/STAT2 9 11.85 hsa01522 Endocrine resistance 5.31E-08 JUN/MAPK8/EGFR/FOS/CDKN1B/CDKN1A/SRC/IGF1 8 19.46 hsa05142 Chagas disease 6.88E-08 NFKBIA/RELA/JUN/MAPK8/FOS/FASLG/FAS/SMAD2 8 18.70 hsa04660 T cell receptor 7.59E-08 NFKBIA/RELA/JUN/ZAP70/MAPK8/MAP3K7/FOS/CD40LG 8 18.34 hsa04380 Osteoclast differentiation 2.83E-07 NFKBIA/RELA/JUN/MAPK8/MAP3K7/FOS/STAT1/STAT2 8 14.90 hsa05164 Influenza A 2.00E-06 NFKBIA/RELA/TBK1/FASLG/FAS/TNFRSF10B/STAT1/STAT2 8 11.09 hsa04060 Cytokine-cytokine receptor interaction 7.43E-05 TNFRSF10C/FASLG/FAS/TNFRSF1B/CD40LG/TNFRSF10B/CD40/TNFRSF10D 8 6.46 hsa05212 Pancreatic cancer 1.76E-07 RELA/MAPK8/EGFR/CDKN1A/SMAD2/STAT1/SMAD4 7 21.96 hsa05210 Colorectal cancer 3.43E-07 JUN/MAPK8/EGFR/FOS/CDKN1A/SMAD2/SMAD4 7 19.40 hsa05235 PD-L1 expression and PD-1 checkpoint pathway in cancer 4.20E-07 NFKBIA/RELA/JUN/ZAP70/EGFR/FOS/STAT1 7 18.75 hsa04658 Th1 and Th2 cell differentiation 5.11E-07 NFKBIA/RELA/JUN/ZAP70/MAPK8/FOS/STAT1 7 18.14 hsa04657 IL-17 5.73E-07 NFKBIA/RELA/JUN/TBK1/MAPK8/MAP3K7/FOS 7 17.75 hsa04933 AGE-RAGE in diabetic complications 8.51E-07 RELA/JUN/MAPK8/CDKN1B/SMAD2/STAT1/SMAD4 7 16.69 hsa04064 NF-kappa B 1.05E-06 NFKBIA/RELA/ZAP70/MAP3K7/XIAP/CD40LG/CD40 7 16.04 hsa04625 C-type lectin receptor 1.05E-06 NFKBIA/RELA/JUN/MAPK8/SRC/STAT1/STAT2 7 16.04 hsa04068 FoxO 4.39E-06 MAPK8/EGFR/CDKN1B/FASLG/CDKN1A/IGF1/SMAD4 7 12.74 hsa04932 Non-alcoholic fatty liver disease 1.01E-05 RELA/JUN/GSK3A/MAPK8/FOS/FASLG/FAS 7 11.12 hsa04217 Necroptosis 1.42E-05 MAPK8/XIAP/FASLG/FAS/TNFRSF10B/STAT1/STAT2 7 10.49 hsa05203 Viral carcinogenesis 6.17E-05 NFKBIA/RELA/JUN/ATF2/CDKN1B/CDKN1A/SRC 7 8.18 hsa05205 Proteoglycans in cancer 6.25E-05 EGFR/FASLG/CDKN1A/FAS/SMAD2/SRC/IGF1 7 8.14 hsa04014 Ras 1.30E-04 RELA/ZAP70/TBK1/MAPK8/EGFR/FASLG/IGF1 7 7.19 hsa05132 Salmonella infection 2.00E-04 NFKBIA/RELA/JUN/MAPK8/MAP3K7/FOS/TNFRSF10B 7 6.70 hsa04151 PI3K-Akt 1.26E-03 RELA/ATF2/EGFR/CDKN1B/FASLG/CDKN1A/IGF1 7 4.71 hsa05120 Epithelial cell signaling in Helicobacter pylori infection 1.93E-06 NFKBIA/RELA/JUN/MAPK8/EGFR/SRC 6 20.43 hsa05140 Leishmaniasis 3.23E-06 NFKBIA/RELA/JUN/MAP3K7/FOS/STAT1 6 18.58 hsa04012 ErbB signaling pathway 5.51E-06 JUN/MAPK8/EGFR/CDKN1B/CDKN1A/SRC 6 16.83 hsa05222 Small cell lung cancer 8.58E-06 NFKBIA/RELA/XIAP/CDKN1B/CDKN1A/BIRC7 6 15.55 hsa05215 Prostate cancer 1.11E-05 NFKBIA/RELA/EGFR/CDKN1B/CDKN1A/IGF1 6 14.75 hsa04722 Neurotrophin signaling pathway 3.34E-05 NFKBIA/RPS6KA5/RELA/JUN/MAPK8/FASLG 6 12.02 hsa05418 Fluid shear stress and atherosclerosis 7.04E-05 RELA/JUN/MAPK8/MAP3K7/FOS/SRC 6 10.29 hsa04062 Chemokine signaling pathway 0.000367 NFKBIA/RELA/GSK3A/SRC/STAT1/STAT2 6 7.45 hsa05202 Transcriptional misregulation in cancer 0.000367 RELA/CDKN1B/CDKN1A/CD40/IGF1/SMAD1 6 7.45 hsa04510 Focal adhesion 0.00045 JUN/MAPK8/EGFR/XIAP/SRC/IGF1 6 7.12 hsa05022 Pathways of neurodegeneration - multiple diseases 0.021689 RELA/TBK1/MAPK8/FASLG/FAS/TNFRSF1B 6 3.01 [122]Open in a new tab 3.4. PPI analysis for these 37 DEPs The analysis identified 44 nodes within the network, comprising 37 known and seven predicted DEPs. Analysis of the interactions between nodes revealed that the network had 318 edges, with the number of edges per node in descending order as follows: JUN (36), TP53 (35), FOS (29), MAPK8 (27), XIAP (26), SPC (25), STAT1 (24), RELA (23), CCND1 (21), and CDKN1A (20) ([123]Fig. 5). Fig. 5. [124]Fig. 5 [125]Open in a new tab (A) Protein-protein interaction network of the 37 DEPs, and (B) graphic representation of the degree of connectivity for the top 19 DEPs. 3.5. Molecular docking and molecular dynamics simulation Following the PPI analysis, molecular docking was conducted to evaluate the binding strength between TMF and these 44 nodes. Using AGFR software (Version 1.2) for active site prediction, the optimal docking sites were determined based on the docking scores. The results indicated that 23 of the 44 proteins exhibited BE of less than 0 kcal/mol with TMF ([126]Table 4), including CDKN1B, SMAD4, NFKBIA, SMAD5, XIAP, ZAP70, IRF3, CD40, HSPB1, FAS, TBK1, BCAR1, IGFBP2, FASLG, CDKN1A, BIRC7, TNFRSF10C, CCND1, JUN, GSK3A, STAT1, MAPK9, and TNFRSF1B. Among these 23 proteins, only MAPK9 and TNFRSF1B had BEs greater than −5 kcal/mol ([127]Fig. 6G). The interactions of six proteins with the strongest binding affinities to TMF are illustrated in [128]Fig. 6A–F. Table 4. Docking parameters of TMF with target proteins. No Targets PDB ID Binding energy (kcal/mol) No Targets PDB ID Binding energy (kcal/mol) 1 JUN 1JNM −5.1 23 IGFBP1 1ZT3 0 2 BIRC7 1OXN −5.6 24 IRF3 1J2F −6.7 3 IGFBP2 2H7T −5.9 25 MAP3K7 4GS6 0 4 RELA 3RC0 0 26 MAPK8 3ELJ 0 5 ATF2 4H36 0 27 MAPK9 3E7O −4.5 6 BCAR1 1WYX −5.9 28 NFκBIA 1NFN −6.9 7 CCND1 2W96 −5.4 29 RPS6KA5 3KN6 0 8 CD40 1LB6 −6.6 30 SMAD1 1KHU 0 9 CD40LG 1ALY 0 31 SMAD4 1YGS −7.4 10 CDK2 1AQ1 0 32 SMAD5 6FZS −6.9 11 CDK3 7XQK 0 33 SRC 1A07 0 12 CDKN1A 2ZVV −5.7 34 STAT1 3WWT −5.0 13 CDKN1B 1H27 −8.0 35 STAT2 6UX2 0 14 FAS 3TJE −5.9 36 TAB1 2J4O 0 15 FASLG 5L19 −5.9 37 TBK1 4IM0 −5.9 16 FOS 1A02 0 38 TNFRSF1B 1CA9 −1.5 17 GSK3A 7SXF −5.1 39 TNFRSF10B 1D0G 0 18 HSPB1 3Q9B −6.1 40 TNFSF10 1D4V 0 19 HSPD1 4PJ1 0 41 TP53 1AIE 0 20 HSPE1 6HT7 0 42 XIAP 1NW9 −6.8 21 IGF1 2DSP 0 43 ZAP70 5O76 −6.8 22 TNFRSF10D AF-[129]Q9UBN6 0 44 TNFRSF10C AF-Q05D7 −5.5 [130]Open in a new tab Fig. 6. [131]Fig. 6 [132]Open in a new tab Molecular docking analysis of TMF with the 44 node proteins. Interactions between TMF and the top six proteins: (A) CDKN1B, (B) SMAD4, (C) NFKBIA, (D) SMAD5, (E) XIAP, and (F) ZAP70. (G) Graphic representation of the proteins with binding energies less than 0 kcal/Mol. To validate the molecular docking results, two proteins including CDKN1B and SMAD5 were selected for molecular dynamics simulation. The findings indicated that both the TMF-CDKN1B and TMF-SMAD5 complexes achieved equilibrium after 50 ns, with RMSD values of 2.07 ± 0.23 Å and 2.46 ± 0.34 Å, respectively, suggesting high stability for both two complexes, with TMF-CDKN1B exhibiting greater stability than TMF-SMAD5 ([133]Fig. 7A). The radius of gyration (Rg) analysis further confirmed the stability of both complexes throughout the simulation ([134]Fig. 7B). Additionally, the fluctuation range of solvent-accessible surface area (SASA) values indicated that TMF binding to CDKN1B and SMAD5 did not induce significant conformational changes ([135]Fig. 7C). Moreover, the root mean square fluctuation (RMSF) values for both complexes were relatively low, ranging from 1 to 4 Å except for the C-terminal region, with the TMF-SMAD5 complex displaying greater flexibility compared to TMF-CDKN1B ([136]Fig. 7D and E). In summary, these results suggest that TMF likely exerts its effects by directly interacting with these two proteins, supporting the reliability of the molecular docking result. Fig. 7. [137]Fig. 7 [138]Open in a new tab Molecular dynamics simulation of TMF-CDKN1B and TMF-SMAD5 complexes. (A) RMSD values of the two complexes over time, (B) Rg values of the two complexes, (C) SASA values of the two complexes, (D) RMSF values of amino acid backbone atoms in TMF-CDKN1B, and (E) RMSF values of amino acid backbone atoms in TMF-SMAD5. 4. Discussion TMF, a potent bioflavonoid belonging to the PMF family, is predominantly found in Citrus plants. Sytrinol, a U.S.-patented food supplement rich in citrus PMFs, is widely utilized globally for the prevention of cardiovascular and cerebrovascular diseases. PMFs have gained significant attention for their medicinal potential, becoming a major focus of research in recent years. However, compared to nobiletin and tangeretin—the two most extensively studied citrus PMFs—TMF has received relatively little attention. This discrepancy may be attributed to the lower concentrations of TMF in Citrus plants, which has historically hindered its availability for research. To address this issue, our previous study [[139]16] successfully overcame the challenge of TMF availability through a semi-synthetic method, demonstrating for the first time that TMF exhibited greater potency in inhibiting HeLa cell growth than nobiletin and tangeretin, and establishing its pharmacokinetic parameters in rats, thereby positioning TMF as a promising therapeutic candidate for CCA treatment. This makes it imperative to further elucidate the mechanisms by which TMF exerts its effects against HeLa cancer cells. In the present study, antibody arrays were employed to assess the proteomic profiles of HeLa cells following TMF treatment. The analysis revealed the upregulation of 17 apoptosis-related proteins, including p21, p27, Fas, and FasL. Both p21 and p27 are critical members of the CIP/KIP family of cyclin-dependent kinase inhibitors, with numerous studies indicating that their expression levels are inversely correlated with cancer invasion and recurrence [[140]24]. The upregulation of these proteins can inhibit CCA cell division and proliferation while promoting differentiation and apoptosis [[141][25], [142][26], [143][27]]. Additionally, whereas Fas and FasL, two pro-apoptotic proteins are known to significantly induce apoptosis in HeLa cells when activated [[144]28,[145]29]. Among the 55 identified phosphorylation sites, several key proteins were affected by TMF treatment. EGFR (Ser1070), ATF2, c-JUN, GSK3A, HSP27, and JNK were downregulated, whereas Sre (Tyr419), Stat1 (Ser727), Stat2 (Tyr689), c-Fos, SMAD1, SMAD2, SMAD4, SMAD5, IκBα (S32), MSK1(S376), NF-κB (S536), TAK1 (S412), TBK1 (S172), ZAP70 (Y292), and MSK2 were upregulated. EGFR, a member of the HER family, is closely associated with tumor development and progression, with clinical trials showing significantly higher EGFR expression levels in CCA tissues compared to adjacent tissues [[146]30,[147]31]. EGFR has been established as an effective target for cancer treatment, and inhibitors such as gefitinib, erlotinib, and osimertinib are frequently prescribed to patients with cancers [[148]32]. The findings from the phosphorylated antibody array suggest that TMF can suppress EGFR activation, mirroring the action mechanism of these EGFR inhibitors. GO functional annotation of the 37 DEPs revealed that TMF significantly impacted various functions in HeLa cancer cells, including the positive regulation of transcription from the RNA polymerase II promoter in response to chemical stimuli, SMAD protein complex formation, and TRAKL binding. Notably, the SMAD protein complex—a core GO term—was further substantiated by PPI analysis and molecular docking, with SMAD4 and SMAD5 exhibiting lower BE with TMF. As key signal transduction molecules downstream of the transforming growth factor β (TGF-β) receptor family, SMAD proteins, particularly SMAD2/3, are crucial in transmitting TGF-β signals from cell-surface receptors to the nucleus [[149]33]. TGF-β itself has been implicated in the formation, invasion, progression, metastasis, and treatment of CCA [[150]34]. Proteomics results showed that SMAD1, SMAD2, SMAD4, and SMAD5 were upregulated, with SMAD2 expression notably reduced in CCA [[151]35]. This study suggests that the TGF-β/SMAD2 signaling pathway may play a significant role in the anti-CCA efficacy of TMF. Further KEGG pathway enrichment analysis indicated that TMF primarily affects the MAPK, TNF, and apoptosis signaling pathways. Similar PMFs, such as tangeretin, sudachitin, and nobiletin, have also demonstrated inhibitory effects on tumor cells through the MAPK signaling pathway [[152][36], [153][37], [154][38]]. Molecular docking results revealed that TMF binds effectively to 21 targets, with BEs below −5 kcal/mol. Despite identifying 10 core nodes via PPI analysis, only half showed direct binding with TMF, suggesting that TMF may influence some targets indirectly by acting on upstream signaling molecules rather than directly interacting with them. Subsequent MD simulations confirmed that TMF directly on targets these proteins. These findings underscore that TMF's efficacy against HeLa cells is achieved through its multi-target, multi-pathway, and multi-functional mechanisms. In our previous study, four permethylated tetramethoxyflavones were evaluated for their anti-cancer bioactivity: 5,6,7,4′-tetramethoxyflavone (TMF), 5,7,3′,4′- tetramethoxyflavone, 3,7,3′,4′-tetramethoxyflavone and 3,5,7,4′-tetramethoxyflavone. Among these, only TMF exhibited strong inhibitory effects against HeLa cancer cells, while the other three compounds showed no significant inhibition of HeLa, A549, HepG2, and HCT116 cancer cells (IC[50] > 100 μM). Despite their lack of efficacy in cancer suppression, the other three PMFs demonstrated substantial biological activity in other areas. For instance, 5,7,3′,4′-tetramethoxyflavone significantly inhibited the expression of α-glucosidase [[155]39] and iNOS [[156]40] in vitro and has shown therapeutic potential in treating inflammatory diseases by inhibiting neuropeptide-stimulated proinflammatory mediator release via mTOR activation [[157]41], as well as by suppressing miR-29a/Wnt/β-catenin signaling through the upregulation of Foxo3a activity [[158]42]. Similarly, 3,7,3′,4′-tetramethoxyflavone was found to suppress inflammation by reducing NO production [[159]43], while 3,5,7,4′-tetramethoxyflavone exhibited potential in preventing cataracts by inhibiting MMP9 and MAPKs [[160]44]. Although these four compounds share the same number of methoxy group and have highly similar chemical structures, their mechanisms of action and biological targets differ significantly. 5. Conclusion This study suggests that the inhibitory effects of TMF against HeLa cancer were primarily mediated through the MAPK, TNF, and apoptosis signaling pathways. These findings provide a foundation for further exploration of TMF's medicinal potential against HeLa cancer. However, this study assessed the effect of TMF on HeLa cell proteomics only at the cellular level in vitro, with the current findings primarily derived from in silico bioinformatics and molecular docking analyses. To strengthen the reliability of these results, additional research is needed using animal models, along with transcriptomics, proteomics western blotting, and qPCR experiments. Funding This work was supported by the Natural Science Foundation of Hainan Province (No. 823MS148), China; the Foundation of Hainan Health Committee (No. 22A200257), China; and the Foundation of Hainan Educational Committee (No. Hnky2023-25), China. Ethics declarations Not applicable. Consent for publication Not applicable. Data availability statement The authors confirm that all data supporting the findings of this study are included in the article. CRediT authorship contribution statement Qiang You: Writing – original draft, Funding acquisition, Conceptualization. Lan Li: Writing – review & editing, Writing – original draft. Haiyan Ding: Visualization, Formal analysis. Youping Liu: Writing – review & editing, Supervision. Declaration of competing 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. References