Identification of key biomarkers and functional pathways in osteosarcomas with lung metastasis: Evidence from bioinformatics analysis.

BACKGROUND: In osteosarcoma, the lung is the most common metastatic organ. Intensive work has been made to illuminate the pathogeny, but the specific metastatic mechanism remains unclear. Thus, we conducted the study to seek to find the key genes and critical functional pathways associated with progression and treatment in lung metastasis originating from osteosarcoma.
METHODS: Two independent datasets (GSE14359 and GSE85537) were screened out from the Gene Expression Omnibus (GEO) database and the overlapping differentially expressed genes (DEGs) were identified using GEO2R online platform. Subsequently, the Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment analysis of DEGs were conducted using DAVID. Meanwhile, the protein-protein interaction (PPI) network constructed by STRING was visualized using Cytoscape. Afterwards, the key module and hub genes were extracted from the PPI network using the MCODE and cytoHubba plugin. Moreover, the raw data obtained from GSE73166 and GSE21257 were applied to verify the expression differences and conduct the survival analyses of hub genes, respectively. Finally, the interaction network of miRNAs and hub genes constructed by ENCORI was visualized using Cytoscape.
RESULTS: A total of 364 DEGs were identified, comprising 96 downregulated genes and 268 upregulated genes, which were mainly involved in cancer-associated pathways, adherens junction, ECM-receptor interaction, focal adhesion, MAPK signaling pathway. Subsequently, 10 hub genes were obtained and survival analysis demonstrated SKP2 and ASPM were closely related to poor prognosis of patients with osteosarcoma. Finally, hsa-miR-340-5p, has-miR-495-3p, and hsa-miR-96-5p were found to be most closely associated with these hub genes according to the interaction network of miRNAs and hub genes.
CONCLUSION: The key genes and functional pathways identified in the study may contribute to understanding the molecular mechanisms involved in the carcinogenesis and progression of lung metastasis originating from osteosarcoma, and provide potential diagnostic and therapeutic targets.

Introduction

Osteosarcoma is the most common malignant bone tumor that mainly originates from the tubular long bones with uneven crowd distribution and high prevalence in children and adolescents.[ As reported, osteosarcoma is prone to develop local infiltration and distant metastasis.[ The most common metastatic organ is lung, with an incidence of approximately 40% to 50%.[ In recent years, due to the development of surgical techniques, multidrug system chemotherapy, precise radiotherapy, and immunotherapy, significant progress has been made in the treatment of osteosarcoma. However, even with aggressive treatment, the prognosis of patients with lung metastasis remains poor, with a 5-year survival rate of only 20% to 30%.[ Therefore, it is essential to understand the precise molecular mechanisms involved in lung metastasis of osteosarcoma and thus develop effective diagnostic and therapeutic strategies.

In recent years, with the constant development of high-throughput technologies and improvements in data analysis techniques, microarray techniques combined with bioinformatics analyses have been extensively applied to identify genetic alterations at the genome level of tumors, particularly in seeking the differentially expressed genes (DEGs). At present, Shi et al[ have analyzed the enrichment pathway and key genes of lung metastasis of osteosarcoma using a single microarray. However, due to the inherent differences in experimental samples, detection platforms, and experimental methods, false-positive rates in independent microarray analysis reduce the reliability of the results obtained. Therefore, in the present study, several original microarray datasets of the Gene Expression Omnibus (GEO) database were downloaded and analyzed to obtain the DEGs between lung metastasis originating from osteosarcoma and non-metastatic osteosarcoma samples. Subsequently, we performed Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses to further investigate the potential functions of these DEGs. Then, the hub genes were identified by the cytoHubba and the other important bioinformatics approaches, including the protein-protein interactions (PPI) network analysis, co-expression network construction, survival analysis, and miRNA-hub genes network construction were also conducted, which may contribute to seeking the key gene signatures and functional pathways involved in lung metastasis of osteosarcomas.

Materials and methods

In the current study, all analyses were based on published datasets that are available in the public repository (http://www.ncbi.nlm.nih.gov/geo/). Thus, no ethical approval or patient consent was required.

Microarray data

The Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database, containing high-throughput gene expression data, chips, and microarrays, is an international public repository.[ After rigorous screening of all relevant datasets that analyze gene expression differences between osteosarcomas with lung metastasis group and non-metastasis group, 2 independent datasets GSE14359[ and GSE85537 were identified from the GEO database. Among them, the GSE14359 dataset was based on the GPL96 platform ([HG-U133A] Affymetrix Human Genome U133A Array) containing 8 lung metastasis samples and 10 non-metastasis osteosarcoma samples. The other dataset utilized the GPL570 platform ([HGU133_Plus_2] consisting of 3 lung metastasis samples and 3 non-metastasis osteosarcoma samples.

Identification of DEGs

GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r), a web-based application, contributes to filtrating DEGs through the Bioconductor's GEO Query and limma R software.[ Therefore, in the present study, the DEGs between the lung metastasis group and the non-metastasis group were screened out by the GEO2R online platform with cut-off criteria of P < .05 and |log fold-change| (|logFC|) >1. Then, a Venn diagram was drawn to obtain the overlapping DEGs from the intersection of the two independent datasets.

GO and KEGG enrichment analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID, http://david.ncifcrf.gov; version 6.8)[ is an online bioinformatics resource website that provides a series of functional annotation tools to study the biological function of genes or proteins. GO, including biological process (BP), cytological component (CC), and molecular function (MF), is a standard semantic vocabulary to define and describe the function of genes and proteins.[ While KEGG is an integrated database resource containing biological interpretation of genome sequences and other high-throughput data.[ Hence, DAVID was utilized to conduct GO annotation and KEGG pathway enrichment analyses to further analyze the characteristic biological functions and potential pathways of DEGs in the present study. P-value < .05 was considered statistically significant.

PPI networks construction and module analysis

The PPIs, existing in the life activities of each cell of the body, regulate many biological phenomena such as replication, transcription, translation, splicing, secretion, cell cycle regulation, signal transduction, and intermediate metabolism. Therefore, it may provide important insights into the carcinogenesis and progression of tumors. In this study, the PPI network of DEGs was constructed using the Search Tool for the Retrieval of Interacting Gene (STRING, http://string-db.org; version 11.0)[ with an interaction score ≥0.4. Afterwards, the PPI network was visualized using Cytoscape 3.7.2 software,[ which is a public software platform for visualizing, integrating, and analyzing molecular interaction networks. And the most key module in the PPI networks was identified using the MCODE plugin of Cytoscape with MCODE score >5, degree cutoff = 2, node score cutoff = 0.2, k-core = 2, and max. depth = 100. Finally, in order to further analyze the functions of genes in the most key module, GO annotation and KEGG enrichment analysis were performed using the DAVID online website.

Hub genes selection and analysis

The top 10 genes of DEGs were defined as the hub genes based on the Maximum Correlation Criteria of the cytoHubba plugin of Cytoscape.[ And the co-expression network of hub genes was constructed by STRING. Meanwhile, in order to validate the expression differences of hub genes between lung metastasis group and non-metastatic group, the expression levels of hub genes obtained from GSE73166 datasets were analyzed and shown by box-plots. In addition, we obtained the online GSE21257 dataset,[ which provided survival data of patients with metastasis and non-metastasis osteosarcoma. According to the expression level of hub genes, these patients were divided into high expression groups (≥50%) and low expression groups. And the prognostic value of the hub genes was confirmed by Kaplan--Meier method using SPSS version 22.0 (SPSS Inc, Chicago, IL). P-value < .05 was considered statistically significant.

Construction of the miRNA-hub gene regulatory network

miRNAs are a type of endogenous noncoding RNA with regulatory functions, which can bind to target genes through base complementary pairing to inhibit translation or induce degradation of target mRNA. The Encyclopedia of RNA Interactomes (ENCORI, http://starbase.sysu.edu.cn/index.php) is a public platform that provides a series of miRNA-hub genes prediction databases to explore miRNA–target gene interactions. In the present study, the interaction relationship between miRNAs and hub genes is established based on the positive results of at least 2 databases in the following miRNA-hub genes prediction databases: miRanda, PITA, PicTar, and TargetScan. Finally, Cytoscape was applied to visualize the miRNA-mRNA regulatory network.[

Results

The raw data of 2 gene expression microarray datasets (GSE14359, GSE85537) were obtained from the GEO database and then subjected to differential expression analysis by the GEO2R online platform. On the basis of the predefined cut-off values, 1742 DEGs were screened from the GSE14359 dataset and 3133 DEGs were identified from the GSE85537 dataset, respectively (Fig. 1A,B). The overlap among the two datasets was identified by drawing a Venn diagram (Fig. 1C), consisting of 96 downregulated genes and 268 upregulated genes.

Figure 1

Identification of differentially expressed genes (DEGs) in 2 gene expression omnibus datasets. (A) Volcano plot of DEGs in GSE14359. (B) Volcano plot of DEGs in GSE85537. Red, blue, and gray color represents the relatively high, low and equal expression of genes in the corresponding dataset, respectively. (C) Venn diagram of overlapping DEGs from the intersection of the two independent datasets.

KEGG and GO enrichment analyses of DEGs

In order to explore the biological functions of DEGs, DAVID was utilized to perform GO annotation and KEGG pathways enrichment analysis. With regard to the biological processes, the DEGs were mainly involved in the regulation of apoptosis, lymphocyte mediated immunity, adaptive immune response, gland development, cell proliferation, enzyme linked receptor protein signaling pathways, and DNA-dependent transcription (Fig. 2A). With respect to cell components, the DEGs were significantly enriched in the extracellular region, extracellular matrix (ECM), extracellular space, plasma membrane part, and cell surface (Fig. 2B). As for the molecular functions, the DEGs were mainly enriched in transcription activator activity, identical protein binding, sequence-specific DNA binding, polysaccharide binding, and RNA polymerase II transcription factor activity (Fig. 2C). Moreover, the results of KEGG pathway enrichment analyses demonstrated that the DEGs were significantly enriched in cancer-associated pathways, adherens junction, ECM-receptor interaction, focal adhesion (FA), and MAPK signaling pathway (Fig. 2D).

Figure 2

GO annotation and KEGG pathways enrichment analyses of the differentially expressed genes (DEGs). (A) Top 20 of the biological process of the DEGs. (B) Top 20 of the cytological components of the DEGs. (C) Top 20 of the molecular function of the DEGs. (D) KEGG signaling pathways of the DEGs. GO = Gene Ontology; KEGG = Kyoto Encyclopedia of Genes and Genomes.

PPI network construction and module analysis

To analyze the functional interactions between DEGs, the PPI network constructed by STRING was visualized using Cytoscape, consisting of 282 nodes and 673 edges (Fig. 3). Subsequently, a total of 8 modules were obtained from the PPI network through the MCODE plugin and the most key module was identified based on MCODE score >5 (Fig. 4A). And the functional analyses of genes in the most key module were conducted by DAVID. As shown in Figure 4B, these genes were significantly enriched in the modification-dependent macromolecule catabolic process, cell cycle, DNA replication, DNA metabolic process, proteolysis, cell division, protein ubiquitination, regulation of programmed cell death, cell proliferation, and protein modification.

Figure 3

The protein-protein interaction network of DEGs. Upregulated genes are marked in red; downregulated genes are marked in blue.

Figure 4

(A) The most key module of the protein-protein interaction network. (B) The biological process analyses of those genes involved in the key module.

Hub gene selection and analysis

The top 10 genes (SKP2, SMC2, CDC25C, SOCS3, PARPBP, DTL, ASPM, ECT2, CENPN, TOP2A) were defined as hub genes after ranking by the Maximum Correlation Criteria of the cytoHubba plugin. The abbreviations, full names, and functions of hub genes are summarized in Table 1. Meanwhile, the coexpression network of these hub genes showed that ASPM was more closely related to other genes (Fig. 5A). Subsequently, we utilized the microarray profile GSE73166 to verify the expression differences of hub genes between lung metastasis group and non-metastatic group. As shown in Figure 5B, the expression levels of hub genes in the lung metastasis group were significantly elevated compared with the non-metastatic group, which was consistent with our results. Next, the original survival data obtained from GSE21257 was used to analyze the prognostic value of hub genes in osteosarcoma by Kaplan--Meier method. And the results demonstrated SKP2 and ASPM were closely related to the poor prognosis of patients with osteosarcoma (Fig. 6). However, it is regrettable that the prognostic data of PARPBP, CENPN, and SMC2 cannot be obtained from the GSE21257 datasets.

Table 1

The detailed information of hub genes.

Gene symbol	Full name	Function
SKP2	S-Phase kinase associated protein 2	SKP2 encodes a substrate recognition component of a SCF (SKP1-CUL1-F-box protein) E3 ubiquitin-protein ligase complex, which mediates the ubiquitination and subsequent proteasomal degradation of target proteins involved in cell cycle progression, signal transduction, and transcription.
ASPM	Abnormal spindle microtubule assembly	ASPM encodes a centrosomal protein that plays an essential part in maintaining the normal function of the mitotic spindle and regulating neurogenesis.
TOP2A	Topoisomerase (DNA) II alpha	TOP2A encodes a DNA topoisomerase that is essential in the regulation of DNA structure and transcription, and it is the direct molecular target of anthacyclines.
SMC2	Structural maintenance of chromosomes 2	SMC2 encodes the central component of the condensin complex, which is required for conversion of interphase chromatin into mitotic-like condense chromosomes.
CDC25C	Cell division cycle 25C	CDC25C encoded a tyrosine protein phosphatase that participated in regulating G2/M progression and mediating DNA damage repair
SOCS3	Suppressor of cytokine signaling 3	SOCS family proteins form part of a classical negative feedback system that mainly regulates cytokine signal transduction.
DTL	Denticleless E3 ubiquitin protein ligase homolog	DTL encodes a substrate-specific adapter of a DCX (DDB1-CUL4-X-box) E3 ubiquitin-protein ligase complex required for cell cycle control, DNA damage response, and translesion DNA synthesis
ECT2	Epithelial cell transforming 2	ECT2 is a guanine nucleotide exchange factor (GEF) of the Rho family members of small GTPases and essential for signal transduction pathways involved in the regulation of cytokinesis.
CENPN	Centromere protein N	CENPN is the vital component of the CENPA-NAC (nucleosome-associated) complex, which plays a central role in assembly of kinetochore proteins, mitotic progression and chromosome segregation
PARPBP	PARP1 binding protein	PARPBP contributes to suppressing inappropriate homologous recombination, thereby playing a central role in DNA repair and the maintenance of genomic stability.

Figure 5

(A) The coexpression network of hub genes. (B) The expression differences of hub genes between lung metastasis group and nonmetastatic group.

Figure 6

Survival analyses of hub genes in patients with osteosarcoma. P < .05 was considered statistically significant.

MiRNA-hub gene regulatory network

To illustrate the regulatory relationships between miRNA and hub genes, the interaction network of miRNAs and hub genes constructed by ENCORI was visualized using Cytoscape. As illustrated in Figure 7, the interaction network consists of 10 hub genes and 227 miRNAs. Through the interactive network degree analysis, ECT2, SOCS3, and PARPBP were apparently the interactive hub genes targeted by most miRNAs, while -miR-495-3p, hsa-miR-340-5p, and hsa-miR-96-5p were the miRNAs that interacted with the most target genes.

Figure 7

Interaction network of hub genes and targeted miRNAs. Hub genes are presented in red circles, whereas targeted miRNAs are shown in blue circles.

Discussion

Due to its malignant properties and invasion behavior, approximately 20% to 30% of osteosarcoma patients present with lung metastasis at the initial diagnosis.[ Even with aggressive treatment, the prognosis of these patients remains poor.[ Thus, exploring the specific molecular mechanisms involved in the progression of osteosarcoma, and then guiding the development of promising therapeutic strategies targeting the molecular pathway will contribute to ameliorating the prognosis of these patients.

In recent years, owing to the wide application of gene chips and high-throughput sequencing technology, it paves the way for our biological research. In the current study, raw data of 2 gene expression datasets were obtained to conduct bioinformatics analysis. As a result, a total of 364 DEGs were identified between lung metastasis originating from osteosarcoma and non-metastatic osteosarcoma. And the results of the functional enrichment analyses suggest that the DEGs mainly involved in some important signaling pathways, including MAPK signaling pathway, adherens junction, ECM-receptor interaction, FA. MAPK signal pathway, consisting of RAF1 (MAP kinase-kinase-kinase), MEK (MAP kinase-kinase), and ERK (extracellular signal-related kinase), play a critical role in promoting tumor proliferation, migration, invasion and angiogenesis.[ Wang et al[ previously reported that the gene mutations associated with lung metastasis of osteosarcoma are mainly concentrated in the MAPK signaling pathway and another study also demonstrated macrophage migration inhibitory factor can promote osteosarcoma growth and lung metastasis through activating the MAPK signaling pathway.[ The ECM, composed of polysaccharides and fibrous proteins, is an important structural component of the tumor microenvironment. Previous studies have shown that dysregulation of the ECM-receptor interaction signaling pathway plays a vital role in regulating tumor invasion and metastasis.[ Not only that, FA promotes cell traction by associating the ECM bound to the transmembrane integrin molecule with the actin cytoskeleton, which also plays an important role in cell migration dynamics.[

To conduct deeper research, the top 10 genes (SKP2, SMC2, CDC25C, SOCS3, PARPBP, DTL, ASPM, ECT2, CENPN, TOP2A) of DEGs were regarded as hub genes based on the predefined criteria. And the expression levels of hub genes were verified to be significantly elevated in the patients with lung metastasis from osteosarcoma. Meanwhile, the results of survival analysis indicated that SKP2 and ASPM were closely associated with poor prognosis of patients with osteosarcoma, which means SKP2 and ASPM may play vital roles in the carcinogenesis and progression of patients with osteosarcoma. SKP2 encodes a substrate recognition component of a SCF (SKP1-CUL1-F-box protein) E3 ubiquitin-protein ligase complex which mainly mediates the ubiquitination and subsequent proteasomal degradation of target proteins involved in cell cycle progression and signal transduction.[ At present, accumulating evidence demonstrates that SKP2 plays a vital role in promoting tumor invasion and metastasis. On one hand, SKP2 could take part in the ubiquitination and degradation of E-cadherin to positively regulate cell migration.[ On the other hand, Skp2 could combine with Myc, Miz1, and p300 to form a transcriptional complex mediating the transcription of RhoA to participate in the invasion and metastasis of tumors.[ Furthermore, SKP2 has also been proven to enhance tumor invasion ability by promoting epithelial-mesenchymal transition (EMT) process and matrix metalloproteinase (MMP) expression.[ On the basis of these encouraging results, Zhang et al[ found knockdown of Skp2 could result in RhoA and MMP-9 inhibition, thereby significantly reducing tumor cell invasion and lung metastasis in an orthotopic mouse model of osteosarcoma. Not only that, SKP2 has also been reported to exert an effect in DNA damage response, DNA repair, tumor angiogenesis, and drug resistance.[ For example, Ding et al[ have demonstrated that SKP2 overexpression was also closely related to the acquisition of methotrexate-resistant and EMT properties in osteosarcoma. Due to the extensive carcinogenic effects of SKP2, research on the anti-tumor activity of SKP2 small molecule inhibitors is in full swing. Among them, flavokawain A (FKA), an effective apoptotic inducer and antiproliferative agent in the tumor, which could lead to downregulation of SKP2, results in osteosarcoma growth retardation and metastasis inhibition.[ However, more preclinical and clinical trials are necessary to explore the application of the specific SKP2 inhibitors in the treatment of osteosarcoma with lung metastasis.

ASPM encodes a centrosomal protein that plays an essential part in maintaining the normal function of the mitotic spindle and regulating neurogenesis.[ Previous studies have prevailingly focused on the relationship between ASPM gene mutations and autosomal recessive primary microcephaly and abnormal neuronal differentiation.[ Nevertheless, except for its role in normal physiological regulatory functions, recent studies have also found that the upregulation of ASPM promotes the metastasis of tumors, including prostate cancer, endometrial cancer, melanomas, bladder cancer, and clear cell renal cell carcinoma.[ Pai et al[ previously pointed out that overexpression of ASPM has energetically influence on prostate cancer cell proliferation and invasion, as it could enhance the Wnt-induced β-catenin transcriptional activity through interacting with disheveled-3, which is a cardinal upstream regulator of Wnt signaling. Likewise, Zhou et al[ found that knockout of ASPM could inhibit tumor proliferation, migration, and metastasis by restraining the Wnt/β-catenin pathway in endometrial cancer. Interestingly, ASPM was situated on chromosome 1q31, an enhancement area that was significantly correlated with invasion and metastasis of tumors.[ The above studies suggested that ASPM could give tumor an impetus to metastasize, but the role of ASPM in lung metastasis of osteosarcoma needs to be further explored.

Although no significant correlation was found between the other hub genes and prognosis in our study, they are still proven to be associated with the invasion and metastasis of tumors. TOP2A gene encodes a DNA topoisomerase that is essential in the regulation of DNA structure and transcription, and it is the direct molecular target of anthracycline drugs. Previous investigations have demonstrated that TOP2A overexpression was involved in tumor proliferation, invasion, and angiogenesis in multiple malignancies, including prostate cancer, pancreatic cancer, soft tissue sarcoma, and osteosarcoma.[ And Nguyen et al[ found that TOP2A overexpression was inclined to poor prognosis of osteosarcoma patients without ERBB2 coamplification. Meanwhile, the therapeutic effect of pirarubicin for osteosarcoma patients with lung metastasis is also closely related to the expression level of TOP2A.[ The oncogene ECT2 is a guanine nucleotide exchange factor (GEF)of the Rho family members of small GTPases and essential for signal transduction pathways involved in the regulation of cytokinesis. Chen et al[ showed that overexpression of ECT2 promoted the metastasis of osteosarcoma by regulating the EMT process. And the downregulation of ECT2 by siRNA and miR-223 could suppress osteosarcoma cell migration and invasion.[ Moreover, Qiu et al[ found that Calycosin could inhibit osteosarcoma cell migration through restraining metastasis-associated IκBα/ECT2 signal pathway, but further research is needed to confirm the results. SOCS family proteins form part of a classical negative feedback system that regulates cytokine signal transduction through the JAK/STAT signaling pathway. Recent studies suggested that Circ_ANKIB1 inhibited the expression of SOCS3 and activated the JAK/STAT signaling pathway by regulating miR-19b, thereby promoting osteosarcoma growth and invasion.[ CDC25C encoded cell cycle regulatory protein that participated in regulating G2/M progression and mediating DNA damage repair. It has been reported that Evodiamine, Flavokawain B, and Ludartin could inhibit proliferation, migration, and apoptosis of osteosarcoma cells by downregulating the expression level of CDC25C.[ At present, the relationship between the rest of hub genes (DTL, CENPN, PARPBP, SMC2) and osteosarcoma remains unclear and the specific role in lung metastasis deserves further exploration.

On the basis of the established miRNA-hub genes interaction network, we found that hsa-miR-340-5p, has-miR-495-3p, and hsa-miR-96-5p were the miRNAs that interacted with the most hub genes. And the role of those miRNAs in the carcinogenesis and progression of osteosarcoma has been reported. For example, miR-340-5p could negatively regulate the Wnt/β-catenin signaling pathway by targeting the STAT3 gene, thereby inhibiting the development of osteosarcoma.[ And overexpression of miR-495-3p could suppress the proliferation, migration, and invasion of osteosarcoma cells by downregulating the CTRP3 gene.[ However, the specific regulation mechanism between these miRNAs and hub genes in patients with lung metastasis of osteosarcoma remained unclear, which need to be explored by further studies.

In the present study, key genes and critical signaling pathways that may be involved in the development of lung metastasis of osteosarcoma were identified via a comprehensive analysis of the 2 gene expression datasets, which may serve as diagnostic and therapeutic targets. Moreover, SKP2 and ASPM were found to be closely related to the poor prognosis of patients with osteosarcoma, which may present as potential prognostic biomarkers in the future. Meanwhile, the interaction network of miRNAs and hub genes illustrates the regulatory relationships of the hub genes and miRNA. And hsa-miR-340-5p, has-miR-495-3p, and hsa-miR-96-5p, which were most closely associated with these hub genes are reported to be involved in the progression of osteosarcoma. These results provide important ideas for a comprehensive understanding of lung metastasis originating from osteosarcoma; however, further studies are needed to validate the current findings and elucidate the specific molecular mechanisms of these genes in lung metastasis originating from osteosarcoma.

Author contributions

Conceptualization: Denghua Liu, Aihua Zhou.

Data curation: Denghua Liu, Rui Zhou, Aihua Zhou.

Formal analysis: Denghua Liu, Rui Zhou, Aihua Zhou.

Investigation: Denghua Liu, Rui Zhou, Aihua Zhou.

Methodology: Denghua Liu, Rui Zhou, Aihua Zhou.

Project administration: Aihua Zhou.

Resources: Denghua Liu, Rui Zhou, Aihua Zhou.

Software: Denghua Liu, Rui Zhou, Aihua Zhou.

Supervision: Denghua Liu, Rui Zhou, Aihua Zhou.

Validation: Denghua Liu, Rui Zhou, Aihua Zhou.

Visualization: Denghua Liu, Rui Zhou, Aihua Zhou.

Writing – original draft: Denghua Liu.

Writing – review & editing: Denghua Liu.