Identification of Differentially Expressed Genes and Pathways in Human Atrial Fibrillation by Bioinformatics Analysis

Abstract

Introduction

Atrial fibrillation (AF) is the most prevalent sustained cardiac arrhythmia, but the molecular mechanisms underlying AF are not known. We aimed to identify the pivotal genes and pathways involved in AF pathogenesis because they could become potential biomarkers and therapeutic targets of AF.

Methods

The microarray datasets of GSE31821 and GSE41177 were downloaded from the Gene Expression Omnibus database. After combining the two datasets, differentially expressed genes (DEGs) were screened by the Limma package. MicroRNAs (miRNAs) confirmed experimentally to have an interaction with AF were screened through the miRTarBase database. Target genes of miRNAs were predicted using the miRNet database, and the intersection between DEGs and target genes of miRNAs, which were defined as common genes (CGs), were analyzed. Functional and pathway-enrichment analyses of DEGs and CGs were performed using the databases DAVID and KOBAS. Protein–protein interaction (PPI) network, miRNA- messenger(m) RNA network, and drug-gene network was visualized. Finally, reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR) was used to validate the expression of hub genes in the miRNA–mRNA network.

Results

Thirty-three CGs were acquired from the intersection of 65 DEGs from the integrated dataset and 9777 target genes of miRNAs. Fifteen “hub” genes were selected from the PPI network, and the miRNA-mRNA network, including 82 miRNAs and 9 target mRNAs, was constructed. Furthermore, with the validation by RT-qPCR, macrophage migration inhibitory factor (MIF), MYC proto-oncogene, bHLH transcription factor (MYC), inhibitor of differentiation 1 (ID1), and C-X-C Motif Chemokine Receptor 4 (CXCR4) were upregulated and superoxide Dismutase 2 (SOD2) was downregulated in patients with AF compared with healthy controls. We also found MIF, MYC, and ID1 were enriched in the transforming growth factor (TGF)-β and Hippo signaling pathway.

Conclusion

We identified several pivotal genes and pathways involved in AF pathogenesis. MIF, MYC, and ID1 might participate in AF progression through the TGF-β and Hippo signaling pathways. Our study provided new insights into the mechanisms of action of AF.

Keywords:

• atrial fibrillation
• bioinformatics analysis
• differentially expressed genes
• pathways
• protein–protein interaction network

Introduction

Atrial fibrillation (AF) is a type of sustained arrhythmia characterized by rapid and disordered atrial electrical activity. AF is one of the leading causes of stroke, heart failure, cardiovascular disease, and sudden death.^1,2 With the increasing age of populations, AF incidence is increasing worldwide.^3,4 However, the pathophysiologic mechanisms of AF are incompletely understood.

Bioinformatic analysis and studies using microarrays to measure gene expression can be employed to screen molecular markers in patients and healthy individuals. This strategy can help to provide novel insights into diseases at multiple levels, ranging from alterations in the copy number at the genome level to gene expression at the transcriptome level (and even epigenetic alterations).^5,6 Microarray studies are used commonly to obtain gene-expression profiles to uncover the pathogenesis of complicated diseases and for biomarker identification.

Using bioinformatics analysis, Wang et al discovered that Zinc finger and BTB domain containing 20 (ZBTB20), Erb-B2 receptor tyrosine kinase 4 (ERBB4), yin yang 1 (YY1), CAMP responsive element binding protein 1 (CREB1), BCL2 interacting protein 3 like (BNIP3L), AT-rich interactive domain-containing protein 1A (ARID1A), protein inhibitor of activated STAT 1 (PIAS1) and forkhead box O3 (FOXO3) may be beneficial for the early diagnosis and future treatment of AF.⁷ Li and coworkers found that cyclin dependent kinase 1 (CDK1), CDK6, and cyclin D3 (CCND3) might be potential biomarkers for the diagnosis and targets for treatment of AF.⁸ Nevertheless, the genes involved in AF pathogenesis have not been reported fully.

In the present study, datasets on the gene-expression profile of AF were downloaded from the Gene Expression Omnibus (GEO) database. Several pivotal differentially expressed genes (DEGs) and pathways were revealed. This information could help to identify biomarkers with high sensitivity in the diagnosis and treatment of AF.

Materials and Methods

Microarray Data

Gene-expression profiles of AF were collected from the GEO database (www.ncbi.nlm.nih.gov/geo).⁹ The GSE31821 dataset includes the atrial tissues of two samples from healthy controls (HCs) and four patients with AF. The GSE41177 dataset contains 38 samples (32 patients with AF and six samples from HCs). The two datasets were both generated using GPL570 (HG-U133_Plus_2) Human Genome U133 Plus 2.0 Array (Affymetrix, Santa Clara, CA, USA).¹⁰ The search term “atrial fibrillation” and the microRNAs (miRNAs) associated with AF were used in the miRTarBase database (http://mirtarbase.mbc.nctu.edu.tw/php/index.php). This is a database on interactions between miRNAs and targets that has been validated experimentally.¹¹ All data in the present study were collected from public databases, so ethical approval from our institution was not required. The study flowchart was presented in Figure 1.

Figure 1 Study flowchart.

Abbreviations: GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Identification of DEGs

We employed the robust multiarray average (RMA) algorithm in R 3.6.3 (www.R-project.org/). Datasets were analyzed with the Affymetrix platform to convert raw microarray data into expression values. This strategy was followed by background correction, quintile normalization, and probe summarization.¹² After merging the microarray datasets (GSE31821 and GSE41177), batch effects were adjusted by the “combat” function of the “sva” package of R using empirical Bayes frameworks.¹³ The DEGs between patients with AF and HCs were screened using the Linear Models for Microarray Data (LIMMA, http://www.bioconductor.org/packages/release/bioc/html/limma.html) package in R Bioconductor.¹⁴ Gene-expression values of | log₂ (fold change, FC) | > 1 and p-value < 0.05 were chosen as cut-off criteria. Volcano maps and heatmaps of DEGs were constructed, respectively.

Analyses of Pathway Enrichment Using the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Databases

Database for Annotation, Visualization and Integration Discovery (DAVID; https://david.ncifcrf.gov/) is an online database that integrates considerable biological data and analytical tools.¹⁵ KEGG Orthology-Based Annotation System (KOBAS, http://kobas.cbi.pku.edu.cn) is an online tool for functional annotation of genes/proteins and functional enrichment of gene sets.¹⁶ For enrichment analyses, KOBAS can accept a gene list or gene-expression data as the input. Then, it generates enriched gene sets, the corresponding name, p-value, or a probability of enrichment and enrichment score based on the results of multiple methods. The GO database is used widely in bioinformatics analysis, and covers three aspects of biology: Biological Process (BP), Cellular Components (CC) and Molecular Function (MF).¹⁷ GO annotations were made using DAVID on DEGs. Pathway analyses of DEGs were done using KOBAS. GO terms and KEGG pathways of biological functions associated with p-value < 0.05 and count number ≥ 2 were considered to be enriched significantly.

Construction of Protein-Protein Interaction (PPI) Network

The miRTarBase database (http://mirtarbase.mbc.nctu.edu.tw/php/index.php) is an experimentally validated database on miRNA–target interactions.¹⁸ The miRNet database (http://www.mirnet.ca) contains information about miRNA–target interactions and displays the association in a visual network.¹⁹ Target genes of miRNAs obtained from the miRTarBase database were selected within the miRNet database. An online tool (https://bioinfogp.cnb.csic.es/tools/venny/index.html) was used to present the intersection of DEGs and the target genes of miRNAs in the form of a Venn diagram. The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (http://string-db.org/) consolidates known data on predicted PPIs for many organisms. This strategy helps to uncover the direct (physical) and indirect (functional) relationships of DEGs.²⁰ Commonly expressed genes were uploaded to the STRING database. Then, a PPI network was established with the minimum required interaction score set as the medium confidence (>0.4). Subsequently, the PPI network was visualized by CytoScape 3.2.1 (https://cytoscape.org/).²¹ The intersection of DEGs and target genes of miRNAs was selected for analyses of the regulation of the miRNA–mRNA network.

Construction of miRNA–mRNA Network in AF

The target miRNAs of commonly expressed genes (CGs) were predicted using the miRTarBase database. The miRNAs of CGs were included once they were found in the reporter assay. Western blotting or real-time reverse transcription-quantitative polymerase chain reaction (RT-qPCR) can also provide strong experimental evidence. Following screening, the miRNA–mRNA network was established further by CytoScape.

Drug-Gene Network Analysis

Consolidating data on gene druggability and drug-gene interactions is the purpose of the Drug-Gene Interaction database. The DGIdb database (http://www.dgidb.org/) was used to predict drug-gene pairings for the co-regulated genes.²² Finally, the drug-gene network was developed using Cytoscape software.

Clinical Experimental Design

We recruited 5 patients with AF and 5 HCs between October 15, 2021, and October 30, 2021 to further confirm the reliability of 9 hub genes selected from the miRNA–mRNA network. The inclusion and exclusion criteria for AF patients and HCs were listed in Appendix 1. The basic characteristic information of patients and HCs were displayed in Supplementary Table 1. Blood samples of all the participants were taken on an empty stomach in the morning on the first day of admission. The clinical blood collection scheme was approved by the Ethics Committee of the Affiliated Hospital of Xuzhou Medical University [Approval number: XYFY2021-KL302]. All volunteers were recruited by the Affiliated Hospital of Xuzhou Medical University and provided written informed consent to the blood samples used in the experiment. The informed consent details are available in Appendix 2.

Reverse Transcription-Quantitative Real-Time Polymerase Chain Reaction (RT-qPCR)

Total RNA was extracted from peripheral blood using TRIzol^® Reagent according to the manufacturer’s instructions (TaKaRa Biotechnology, Shiga, Japan). The iScriptTM cDNA Synthesis Kit (Bio-Rad Laboratories, Hercules, CA, USA) was then used to reverse transcription 600ng of extracted total RNA to synthesis complementary DNA (cDNA).

The ABI-384 Real-Time PCR Detection System (Applied Biosystems, Foster City, CA, USA) was used to quantify mRNA expression with the appropriate proportions of SYBR™ Green qPCR Master Mix kit (Bio-Rad Laboratories), cDNA and deionized water using the comparative quantification method (2^−ΔΔCT). β-actin was used as the internal reference. A list of the primer sequences for the genes are shown in Supplementary Table 2.

Statistical Analyses

All data were performed using SPSS 22.0 (IBM Corporation, Armonk, NY) and GraphPad Prism 7 (GraphPad, San Diego, CA). Continuous variables were expressed as mean ± standard deviation, and classified variables were represented by N (%). The chi-square test was used to calculate the differences in the rates between patients and HCs. Independent samples t test was used to analyze differences in general characteristics and the hub gene expression levels between patients and HCs. P-value < 0.05 was selected to indicate a statistically.

Results

Identification of DEGs

A total of 65 DEGs were screened from the integrated dataset: expression of 49 genes was upregulated and expression of 16 genes was downregulated between AF patients and HCs (Figure 2A). The cluster heatmaps of all DEGs are shown in Figure 2B, and all DEGs are displayed in Supplementary Table 3.

Figure 2 Volcano plots and Heatmap for the DEGs identified from the integrated dataset.

Notes: (A) Red points represent upregulated DEGs and green points denote downregulated DEGs. (B) Each row represents the DEGs and each column represents one of the samples of normal samples or AF samples. Red and blue represent upregulated and downregulated DEGs, respectively.

Abbreviations: DEG, differentially expressed gene; AF, atrial fibrillation.

Pathway Enrichment of DEGs Using GO and KEGG Databases

Analyses using the GO database showed that the DEGs were enriched mainly in BPs, including “innate immune response”, “neutrophil chemotaxis”, and “negative regulation of the apoptotic process”. With regard to CC, DEGs were enriched in “extracellular exosome”, “extracellular region”, and “cell surface”. With respect to MF, the DEGs were enriched significantly in “RAGE receptor binding”, “arachidonic acid binding”, and “protein binding” (Supplementary Table 4). Pathway enrichment-analyses using the KEGG database showed that DEGs were expressed mainly in the “TGF-β signaling pathway” and “Hippo signaling pathway” (Figure 3, Supplementary Table 5).

Figure 3 Terms of BP, CC and MF in pathway-enrichment analyses for DEGs using the GO database.

Notes: Red lines represent biological process. Brown lines represent cellular components. Yellow lines represent molecular function. Green lines represent the KEGG database.

Abbreviations: DEGs, differentially expressed genes; BP, biological process; CC, cellular component; MF, molecular function; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Analyses of PPI Networks

To further explore the molecular mechanism of AF, target genes of mRNAs were predicted using the miRNet database. A total of 9777 genes were predicted as target genes of miRNAs, of which 33 genes overlapped with DEGs (Figure 4A). For a combined score > 0.4, the STRING database was utilized to construct a PPI network with the intersection of 33 DEGs. After removal of isolated nodes, the PPI network was constructed with 15 nodes and 13 edges using CytoScape (Figure 4B). Figure 4C–F and Supplementary Tables 6 and 7 reveals that these 33 common genes were enriched mainly in the TGF-β signaling pathway and Hippo signaling pathway.

Figure 4 Venn diagram, PPI network and the functional terms and pathways enriched for common genes.

Notes: (A) Venn diagram showing common genes in the target genes of DEGs and mRNAs. (B) Red nodes represent upregulated genes, and blue nodes represent downregulated genes. (C–F) Enrichment analyses of common genes using GO and KEGG databases.

Abbreviations: PPI, Protein–protein interaction; DEGs, differentially expressed genes; mRNA, messenger RNA; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Analyses of Regulatory miRNA–mRNA Networks

A regulatory miRNA–mRNA network was constructed. It consisted of 91 nodes and 120 edges, including 9 target mRNAs and 82 miRNAs (Figure 5A). Nodes of degree ≥ 2 in the mRNA–miRNA network are shown in Table 1. Also, hsa-miRNA-212-3p, hsa-miRNA-494-3, hsa-miRNA-335-5, hsa-miRNA-451a and hsa-miRNA-146a-5p interacted with ≥2 mRNAs simultaneously (Figure 5B).

Table 1 The Nodes with Degree of ≥ 2 in the mRNA-miRNA Network

mRNA		miRNA
Symbol	Degree	Symbol	Degree
MYC	55	hsa-miR-146a-5p	3
CXCR4	15	hsa-miR-494-3p	2
SOD2	4	hsa-miR-451a	2
MIF	3	hsa-miR-335-5p	2
ID1	3	hsa-miR-212-3p	2
DDIT4	3
S100A12	2
PDK4	2

Abbreviations: mRNA, messenger RNA; miRNA, microRNA.

Figure 5 miRNA–mRNA integrated network.

Notes: Red balls represent upregulated mRNAs. Green balls represent downregulated mRNA. Blue triangles are predicted miRNAs.

Abbreviations: mRNA, messenger RNA; miRNA, microRNA.

Drug-Gene Network Analysis

For the 9 common genes that interact with predicted miRNAs, 79 drug-gene pairs were acquired, including 7 common genes (including Macrophage migration inhibitory factor (MIF), MYC proto-oncogene, bHLH transcription factor (MYC), Pyruvate Dehydrogenase Kinase 4 (PDK4), C-X-C Motif Chemokine Receptor 4 (CXCR4), S100 Calcium Binding Protein A12 (S100A12), DNA Damage Inducible Transcript 4 (DDIT4) and Superoxide Dismutase 2 (SOD2) and 77 drugs (Figure 6).

Figure 6 The drug-gene network.

Notes: Red balls, green ball, and grey squares represent upregulated mRNAs, downregulated gene, and drugs, respectively.

Abbreviation: mRNAs, messenger RNA.

Validation of Hub Genes

The RT-qPCR was used to confirm the expression levels of MIF, MYC, CXCR4, SOD2, inhibitor of differentiation 1 (ID1), DDIT4, S100A12, Ras related dexamethasone induced 1 (RASD1), and PDK4. As demonstrated in Figure 7, the expression levels of MIF, MYC, ID1, and CXCR4 were significantly higher in AF samples, as compared to HCs (all p-value <0.05). The expression of SOD2 was significantly downregulated in AF patients when compared to HCs (p-value <0.05). The findings of the RT-qPCR analysis in our research were consistent with the results of bioinformatics analysis. However, there was no difference in the expression of PDK4, S100A12, RASD1, or DDIT4 between AF patients and HCs.

Figure 7 Validation of hub genes in AF patients.

Notes: The expression level of MIF, MYC, CXCR4, ID1, SOD2, RASD1, PDK4, DDIT4 and S100A12 determined by RT-qPCR.

Abbreviations: AF, atrial fibrillation; MYC, MYC proto-oncogene, bHLH transcription factor; CXCR4, C-X-C Motif Chemokine Receptor 4; SOD2, Superoxide Dismutase 2; MIF, Macrophage migration inhibitory factor; ID1, inhibitor of differentiation 1; DDIT4, DNA Damage Inducible Transcript 4; S100A12, S100 Calcium Binding Protein A12; RASD1, Ras Related Dexamethasone Induced 1; PDK4, Pyruvate Dehydrogenase Kinase 4.

Discussion

The mechanism of AF is incompletely understood. Therefore, the early diagnosis of AF is the focus of basic scientific and clinical research.¹ Bioinformatics analysis of AF has focused mainly on DEGs using GEO datasets. However, there are few datasets on AF. We concentrated on the miRNAs related to AF as demonstrated by clinical trials, RT-qPCR, Western blotting, and luciferase assays. Moreover, the samples we used were from human atrial tissues. To research the biological function of common genes, we conducted enrichment analysis, constructed PPI, miRNA–mRNA and drug-gene networks. We speculated that the TGF-β and the Hippo signaling pathway might have vital roles in AF progression revealing from the KOBAS database. The genes associated with the TGF-β signaling pathway and the Hippo signaling pathway included MIF, MYC, and ID1, which were illustrated in the PPI network as well.

The TGF-β signaling pathway regulates the proliferation and differentiation of cells as well as fibrosis development.²³ Studies have demonstrated AF pathogenesis is related mainly to structural remodeling of atria involving cardiomyocyte apoptosis, inflammation, and activation of fibrotic pathways via fibroblasts.^24,25 Moreover, the relationship among TGF-β, fibrosis and AF, as well as partial regulatory mechanisms, have been shown. Wang et al²⁶ demonstrated that long non-coding RNA predicting cardiac remodeling (lnc LIPCAR) might regulate atrial fibrosis via the TGF-β signaling pathway in samples from AF patients. Besides, Qiu and collaborators²⁷ found that salvianolate could suppress atrial interstitial fibrosis by inhibiting the TGF-β signaling pathway in AF rats. Our results confirmed the function of the TGF-β signaling pathway in AF pathogenesis.

The Hippo signaling pathway regulates the survival, proliferation, and apoptosis of cells. It has an important role in development/regeneration and diseases of the heart.^28,29 However, few researchers have focused on the function of the Hippo signaling pathway in AF. Zhang et al³⁰ and Li and coworkers³¹ suggested the role of the Hippo signaling pathway in AF regulation by bioinformatics analysis. Dysregulation of the Hippo signaling pathway in the cardiovascular system has been reported to be associated with myocardial infarction (MI), cardiac hypertrophy and vascular remodeling.³² Ramjee et al³³ found that deficiency of two core effectors of the Hippo signaling pathway, Yes-associated protein (YAP) and transcription regulator protein 1 (TAZ), could lead to the post-MI pericardial inflammation and myocardial fibrosis in mice, which suggested the profound effect of the Hippo signaling pathway in MI. Yang and collaborators³⁴ found that YAP could promote cardiomyocyte growth, and cause hypertrophy in postnatal hearts. He et al³⁵ discovered that the Hippo/YAP signaling pathway alters the growth, death, and migration of vascular smooth muscle cells and endothelial cells, which contributes to vascular remodeling in cardiovascular diseases. In the progression of cardiac hypertrophy and vascular remodeling, fibrosis is vital.^36,37 In this study, we showed that the Hippo signaling pathway participates in AF. Upon consideration of PPI networks, we found that MYC, ID1, and MIF were involved in the pathways mentioned above, thereby suggesting their roles in AF pathogenesis.

MIF, MYC, ID1, and CXCR4 had substantially higher levels of expression and SOD2 was downregulated in patients with AF compared to HCs. MIF is a pleiotropic inflammatory cytokine which shows high expression in patients with AF. Several studies have demonstrated its detailed mechanism of action in AF progression. Cheng et al³⁸ demonstrated that MIF might promote cardiac fibrosis in AF patients by activating the TGF-β signaling pathway. The relationship between MIF and Hippo is not known. We speculate that, because Hippo might be related to AF, and MIF was found to take part in the Hippo signaling pathway, then MIF might function in AF through the Hippo signaling pathway. MYC is a proto-oncogene implicated in regulation of the formation and metastasis of tumors. Yee and collaborators³⁹ discovered that polymorphisms in MYC-related genes lead to various bleeding complications in patients with a stable response to warfarin (which is used commonly in AF). Zhang and coworkers⁴⁰ found that TGF-β expression could be regulated by long noncoding RNA regulator of reprogramming through regulation of c-MYC expression. Xiao et al⁴¹ found that Yap (a pivotal molecule in the Hippo signaling pathway) activates MYC directly in response to cardiac injury. Based on the relationship among genes and pathways in cardiac models, genes and AF, AF and pathways, we postulate that MYC might participate in AF pathogenesis. ID1 is a regulator of the cell cycle and cell differentiation.⁴² ID1 expression in the heart is confined to non-myocardial layers. ID1 is essential for development of the neonatal heart and has important roles in embryogenesis.⁴³ Sato and collaborators⁴⁴ found that ID1 blocks TGF-β-induced apoptosis of mesangial cells. Lin et al⁴⁵ revealed that ID1 could mediate the fibroblast response to TGF-β in pulmonary fibrosis. Fibrosis is related to TGF-β and Hippo signaling pathways. Therefore, we speculate that ID1 might participate in TGF-β and Hippo signaling pathways in AF. Chen M found that by increasing antioxidant genes such as SOD2 and UCP3, fibroblast growth factor 21 successfully reduced atrial remodeling by lowering oxidative stress.⁴⁶ In addition, Wang et al revealed that CXCR4 expression is increased in chronic AF patients with mitral valve dysfunction, that it is related with atrial remodeling.⁴⁷

However, no difference in the expression of PDK4, S100A12, RASD1, or DDIT4 was seen between AF patients and HCs. Raman K demonstrated that PDK4 is a useful biomarker for discriminating between atrial fibrillation and sinus rhythm.⁴⁸ Via the least absolute shrinkage and selection operator model, Liu L identified that S100A12 is potentially involved in the molecular mechanisms of AF development.⁴⁹ The relationship between RASD1, DDIT4 and AF has not been reported. However, the involvement of the RASD1 and DDIT4 genes in AF requires more investigation.

Our study had three main limitations. First, we created the PPI network with the minimum required interaction score of medium confidence (>0.4) in order to obtain more proteins related with AF. Second, we concentrated only on reverse regulation of miRNA–mRNA pairs and ignored the more complicated mechanisms of miRNA–mRNA pairs. Finally, the miRNA and mRNA we obtained were not from identical samples.

Conclusions

In the study, we identified 9 pivotal genes and several related pathways involved in AF pathogenesis. We further verified experimentally that MIF, MYC, CXCR4, and ID1 were up-regulated and SOD2 was down-regulated in patients with AF. Based on the miRNA-mRNA network, 11 pairs of miRNA-mRNA pairs associated with AF were identified, including miR-335-5p/CXCR4, miR-335-5p/MYC, miR-494-3p/CXCR4, miR-494-3p/MYC, miR-451a/MYC, miR-451a/MIF, miR-146-5p/CXCR4, miR-146-5p/MIF, miR-146-5p/S100A12, miR-212-3p/SOD2, miR-212-3p/MYC, that might regulate development and progression of AF. In addition, MIF, MYC, and ID1 might participate in AF progression through TGF-β and Hippo signaling pathways. Our study provides new insights into AF mechanisms.

Data Sharing Statement

The data used to support the findings of this study are available from the corresponding author upon request.

Author Contributions

All authors contributed significantly to the conception and design, data acquisition, and data analysis and interpretation; participated in the drafting of the article or critically revised it for important intellectual content; agreed to submit to the current journal; gave final approval of the version to be published; and agreed to be accountable for all aspects of the work.

Disclosure

The authors declare that they have no competing interests.

Additional information

Funding

This work was supported by grant from the Xuzhou Science and Technology Bureau to Defeng Pan [grant number: KC20097].