Network and pathway-based analyses of genes associated with osteoporosis.

Osteoporosis (OP) is a disease characterized by bone mass loss, bone microstructure damage, increased bone fragility, and easy fracture. The molecular mechanism underlying OP remains unclear.In this study, we identified 217 genes associated with OP, and formed a gene set [OP-related genes gene set (OPgset)].The highly enriched GOs and pathways showed OPgset genes were significantly involved in multiple biological processes (skeletal system development, ossification, and osteoblast differentiation), and several OP-related pathways (Wnt signaling pathway, osteoclast differentiation, steroid hormone biosynthesis, and adipocytokine signaling pathway). Besides, pathway crosstalk analysis indicated three major modules, with first module consisted of pathways mainly involved in bone development-related signaling pathways, second module in Wnt-related signaling pathway and third module in metabolic pathways. Further, we calculated degree centrality of a node and selected ten key genes/proteins, including TGFB1, IL6, WNT3A, TNF, PTH, TP53, WNT1, IGF1, IL10, and SERPINE1. We analyze the K-core and construct three k-core sub-networks of OPgset genes.In summary, we for the first time explored the molecular mechanism underlying OP via network- and pathway-based methods, results from our study will improve our understanding of the pathogenesis of OP. In addition, these methods performed in this study can be used to explore pathogenesis and genes related to a specific disease.

Introduction

Osteoporosis (OP) is a systemic skeletal disease characterized by bone mass loss and bone microstructure damage which results in bone fragility and increased fracture risk.[ As a major threat to elderly people, the prevalence of OP increases with continued aging of the population.[ It is estimated that in the People's Republic of China, the mean prevalence of OP in older adults is 15.7%, and there were about 202.43 million people aged 60 years and older at the end of 2013, therefore, OP will cost a large financial burden.[

OP is a chronic condition. Osteoporotic fracture, which is the most serious clinical consequence of OP, increases gradually with the increasing age and occur at trabecular bone. Osteoporotic fractures majorly occur in spine, hip, distal forearm and proximal humerus,[ while osteoporotic fractures also occur in other sites, such as the humerus, ribs, tibia, and pelvis femoral fractures.[ Osteoporotic fractures, particularly hip fracture, cause a large economic burden.[ Therefore, OP has attracted much attention in many countries around the world.

OP is caused by the imbalance between bone formation and bone resorption,[ and considered to be resulted from a complicated interaction of environmental and genetic factors. Many studies point out that the environmental factors, such as diet and lifestyle factors, was the risk of OP,[ while other studies revealed that the genetic factor determined the variation in BMD and OP. The identification of the possible essential proteins, such as vitamin D receptor gene,[ transforming growth factor-β1 (TGF-β1) gene,[ osteoprotegerin (OPG) gene,[ collagen type I α1 gene,[ and the estrogen receptor-α gene,[ has improved our understanding of the pathology of OP. Ahmad et al showed that VDR gene FokI and BsmI polymorphism is significantly associated with low-bone mineral density, and f allele of VDR FokI gene may be used as an important risk factor for OP.[ Several studies have already assessed the connection between OPG gene polymorphism, which might be a candidate gene for OP predisposition, and bone mineral density (BMD).[ With the development of biomedical science, a number of studies have accelerated the process of understanding of OP, from the perspective of animal models,[ gene analysis,[ genome-wide association studies,[ and systems biology.[ However, the molecular mechanisms of OP were still far to understand.

Most studies have focused single gene and/or pathway to elucidate the molecular mechanisms of OP, few studies are made to implement systematic bioinformation analyses to elucidate the molecular mechanisms. For example, miR-217 has been reported to promote proliferation and osteogenic differentiation of BMSC in steroid-associated osteonecrosis via inhibiting Dikkopf-1 (DKK1) and activating Wnt signaling pathway.[ MiRNA-10b promotes osteogenic differentiation of MSCs and inhibits adipogenic differentiation by repressing SMAD2, as well as stimulating transforming growth factor beta (TGF-β) signaling pathway.[ The GDF11-FTO-PPARγ axis controlled the shift from osteoporotic MSC to adipocyte and inhibited bone formation during OP.[ PI3K-Akt signaling pathway is reported to promote osteoblast proliferation, differentiation and bone formation and inhibit process of OP.[ It is now recognized that most complicated psychiatric phenotypes are influenced by numbers of genes with small effects, instead of single genes with large effects.[ Therefore, it will provide more useful insights beyond the conventional single-gene analyses to analysis multiple genes in a system biological level.[ In current study, we firstly collected OP-associated genes comprehensive, which were reported in PUBMED. Second, we identified the significant biological terms and pathways within the genetic factors using GO and KEGG enrichment analyses. Third, we analyzed the OP-related genes in the context of human protein-protein interaction (PPI) network. At last, OP-specific genes network, which play crucial roles in OP, was inferred using note degree of the network and the k-core algorithm. This study will advance our understanding of the pathological mechanisms of OP from the systems biological level.

Materials and methods

Collection of OP-associated genes

We retrieved the human genetic association studies deposited in PUBMED (http://www.ncbi.nlm.nih.gov/pubmed/), and collected OP-associated genes. The collection of genes was performed according to the method described in previous studies.[ Briefly, the reports related to OP were queried with the term (Osteoporosis [MeSH]) and (polymorphism [MeSH] or genotype [MeSH] or alleles [MeSH]) not (neoplasms [MeSH]), and a total of 1370 articles were retrieved associated with OP by November 30, 2018. Then, we reviewed the abstracts of initial publications and collected the genetic association studies of OP. We narrowed our selection via focusing on the selected publications, which reported significant associations between genes and OP. We reduced the number of false-positive finding by excluding the publications, which reported negative or insignificant associations. We reviewed the full texts of the selected publications and ensured that the content supported the conclusions. The genes, which were reported to be significantly associated with OP in these studies, were selected for this study. An ethical approval was not needed, since it was not involved with human or animals.

GO and pathway enrichment analysis of OP-related genes

We performed GO and pathway enrichment analysis of the OP-related genes using WebGestalt[ and ToppGene,[ respectively. WebGestalt is used to evaluate the significantly enriched GO terms in the OP-related genes gene set (OPgset). Pathway analysis of OPgset genes was carried out to find out the important pathways using ToppGene. Pathways with one or more genes overlapping the OPgset were selected. The significant GO terms and pathways were identified by Fisher exact test, and FDR value of P < .05 were considered to be significantly enriched.

Pathway crosstalk analysis

Pathway crosstalk analysis was performed to explore the interactions between any pair of significantly enriched pathways. Similar to Ref.,[ we computed two measurement and evaluated the overlap between any pair of pathways, that is, the Jaccard Coefficient and the Overlap Coefficient , where A and B are the lists of genes included in the two tested pathways.[ The significant enriched pathways containing more than two candidate genes were selected. Then, the number of shared candidate genes between any pair of pathways was counted and pathway containing more than two overlapped genes was included. JC and OC value of pathway pairs were calculated and ranked. The software Cytoscape was used to visualize the selected pathway crosstalk.[

Network analysis

Network construction: In order to investigate the function of OP-related genes at the protein level, we constructed OPgset-related network via a biological database, Search Tool for the Retrieval of Interacting Genes/Proteins (STRING). The PPI data were visualized via Cytoscape (Version 3.6.1).

Defining network topological feature set: We computed three measures, such as Degree, Edge betweenness, and K-core, for assessing topological feature of each node i in OPgset-related network. Degree, which is defined as the number of linkers to node, was used to measure the topological importance of protein in the network.[ Edge Betweenness is defined as the frequency of an edge that places on the shortest paths between all pairs of vertices in network.[ The edges with highest betweenness values are most likely to lie between sub-graphs. K-core analysis, which is used to measure the centrality of node i, is an iterative process in which the nodes are removed from the networks in order of least-connected.[ We generated the k-core sub-network of OPgset-related genes using Cytoscape plug-in MCODE. GO analysis of sub-networks was reconstructed using Cytoscape plug-in BinGO.

Results

Identification of OP-related genes

There were more than 500 studies were collected. In these publications, 217 genes were reported to be significantly associated with OP and formed a gene set (OPgset) for following analysis (Table S1). Among them were five hormones, that is, GH1, GHRH, PTH, GNRH1, and PTHLH, several hormone receptors, that is, ESR1, ESR2, vitamin D receptor, ESRRG, FSHR, GHR, NR3C1, CALCR, and CASR, three bone morphogenetic proteins, that is, BMP2, BMP4, and BMP15, five Wnt family members, that is, WNT1, WNT3A, WNT4, WNT5B, and WNT16, four Wnt antagonists, that is, SOST, SFRP1, SFRP4, and DKK1, several TNF members, that is, TNF, TNFRSF11A, TNFRSF11B, THFRSF1B, and TNFSF11, four transcription factors, that is, RUNX2, MRTF, SP7, and TWIST1. Several genes were involved in the functions associated with nitric oxide synthesis, adipo-metabolism-related genes (e.g., ADIPOQ, ADIPOR1, and APOE), glutathione metabolism relate genes (e.g., GPX1, GSR, GSTM1, GSTM3, and GSTP1), as well as immune system (e.g., IL1A, IL-1B, IL6, IL10, IL-15, IL-16, IL-17A, IL17B, and IL 21R). These data suggested that the genes significantly associated with OP were multifarious.

Biological functions enriched in OPgset

GO enrichment analysis was performed to investigate the biological function of 217 genes in OPgset (Table S2). Significantly enriched GO terms, in the OPgset genes, include those associated with bone development, osteoblast differentiation, osteoclast differentiation, adipogenesis, and vitamin D function. GO terms associated with bone development (e.g., Skeletal system development, bone morphogenesis, bone maturation, bone mineralization involved in bone maturation, bone growth, bone resorption, bone cell development, ossification, regulation of ossification, bone mineralization, and bone development) were enriched in genes in OPgset. Terms directly related to osteoblast differentiation (e.g., osteoblast differentiation, osteoblast development, osteoblast proliferation, regulation of osteoblast differentiation, and regulation of osteoblast proliferation) and osteoclast differentiation (osteoclast differentiation, macrophage cytokine production, and positive regulation of osteoclast differentiation) were included. These data were accord with the previous reports that an imbalance between bone resorption and formation is the pathophysiological of OP. In addition, GO terms associated with adipogenesis (e.g., regulation of fat cell differentiation, fat cell differentiation), vitamin D function (response of vitamin D, vitamin D metabolic process, cellular response to vitamin D, vitamin D biosynthetic process, vitamin D receptor signaling pathway), signaling pathways (e.g., MAPK, BMP signaling pathway, Wnt signaling pathway), and immune function (e.g., immune system development, leukocyte proliferation, leukocyte differentiation, interleukin-6 production, lymphocyte differentiation, T cell differentiation, and interleukin-17 production) were also enriched in these genes.

Pathway enrichment analysis in OPgset

Fifty-seven significantly enriched pathways related to OP were identified (Table S3). Consistent with previous studies,[ a number of pathways related to signaling pathway, for example, Wnt signaling pathway (ranked 2nd in Table S3), Hippo signaling pathway, TNF signaling pathway, MAPK signaling pathway, and NF-kappa B signaling pathway were enriched in OPgset. In addition, bone development-related pathways were identified, for example, osteoclast differentiation, ECM-receptor interaction, all of which were closely involved in bone development. Furthermore, immune processes including cytokines and inflammatory response, IL-10 anti-inflammatory signaling pathway, signal transduction through IL-1R, TH17 cell differentiation, IL-17 signaling pathway, and IL-5 signaling pathway were also enriched, indicating that the immune system play an important role in the pathological process of OP. Further, pathways related to adipogenesis, such as adipocytokin signaling pathway and PPAR signaling pathway, were also enriched in the OPgset genes, in accordance with previous studies.[

Crosstalk analysis of significantly enriched pathways

A pathway crosstalk analysis was performed to understand how lists of significantly enriched pathways interact with each other. This analysis was on account of the assumption that two pathways were considered to crosstalk if they shared part of OPgset.[ Total of 48 pathways was selected according to the criterion: pathway contain at least three members in OPgset and pathway shared ≥2 genes with one or more other pathways. Similar to Ref.,[ pathway crosstalk was formed using the pathway pairs (edges) of selected pathways, and the average scores of coefficients JC and OC was used to measure the overlapping level between two pathways (Fig. 1). We calculated all node degrees of the signaling pathways. According to the previous study by Han et al,[ in which they defined a hub as a node degree exceeding 5, we found that 32 nodes could be chosen as hub nodes and the results are shown in Table S4. As shown in Figure 1, there are three major modules in the pathway network, in which pathways shared more interactions with each other. The first module mainly included bone development-related signaling pathways, such as osteoclast differentiation, MAPK signaling pathway, TGF-β signaling pathway, PI3K-Akt signaling pathway and HIF-1 signaling pathway, as well as immune system-related pathways, such as cytokine-cytokine receptor interaction, IL-10 anti-inflammatory signaling pathway, IL-17 signaling pathway, TH17 cell differentiation, and signal transduction through IL-1R. The second module included Hippo signaling pathway, signaling pathways regulating pluripotency of stem cells, Wnt signaling pathway, Melanogenesis, and mTOR signaling pathway. The third module was primarily dominated by the metabolic pathways of hormone or drug, including Drug metabolism-cytochrome P450, Glutathione metabolism, Steroid hormone biosynthesis, Ovarian steroidogenesis, Free Radical Induced Apoptosis, and Metabolism of xenobiotics by cytochrome P450. At the same time, the first and second modules were connected via several pathway interactions.

Figure 1

Pathway crosstalk among OSgset-enriched pathways. Nodes represent pathways, and edges represent crosstalk between pathways. Edge-width corresponds to the score of specific pathway pair. Larger edge-width indicated higher score. Node-size corresponds to the degree of pathway. Larger node-size indicated higher degree.

Functional analysis of differentially expressed genes

In order to investigate the function of OP-related genes at the protein level, we constructed OPgset-related network via a biological database, Search Tool for the Retrieval of Interacting Genes/Proteins (STRING). The PPI data were visualized via Cytoscape (Version 3.6.1). The highest confidence score (0.9) was adopted to evaluate the protein interactions for OPgset. The network included 143 nodes and 336 edges (Fig. 2). We calculated all node degrees of genes and the results are shown in Table S5. The top 10 note degree genes were TGFB1, IL6, WNT3A, TNF, PTH, TP53, WNT1, IGF1, IL10, and SERPINE1.

Figure 2

Protein-protein interaction networks of OSgset. The PPI data were obtained from STRING. A confidence score that calculated for all protein interactions based on experimentally and computationally interaction was set as the highest (>0.9). Then, we applied Cytoscape (Version 3.6.1) to visualize the networks. Nodes represent genes, and edges represent interaction between genes. Edge-width corresponds to the combined score of specific genes pair. Larger edge-width indicated higher score. Node-size corresponds to the degree of pathway. Larger node-size indicated higher degree. PPI = protein-protein interaction.

Functional prediction of sub-networks

Via K-core analysis, we extracted three sub-networks from the global network and analysis the function of genes. As shown in Table 1, the “1” sub-network was composed of 9 genes nodes and 36 edges. The “2” sub-network was composed of 8 genes nodes and 28 edges. The “3” sub-network was composed of 11 genes nodes and 36 edges. We evaluated the enrichment significance of GO terms in the subnets. The GO term with FDR value less than 0.01 were considered statistically significant. The results of GO and pathway analysis revealed that 157 GO terms and 2 pathways were enriched in the “1” sub-network (Table S6 and S7), 155 GO terms and 12 pathways were enriched in the “2” sub-network (Table S8 and S9), 264 GO terms and 11 pathways were enriched in the “3” sub-network (Table S10 and S11). The GO terms interaction networks are shown in Figure 3A, 4A, and 5A. The top 10 significant GO terms are shown in Figures3B, 4B, and 5B, and the significant enriched pathways are shown in Figures3C, 4C and 5CB.

Table 1

Information of sub-networks.

Figure 3

The GO terms interaction network of “1” sub-network (A). GO terms displayed as an interaction network using Cytoscape plug-in BinGO. Yellow nodes: nodes with P-value < 0.01 and Benjamini corrected P-value < 0.01. Biological function and pathway analysis of “1” sub-network. (B) The significant changes in the GO biological process. (C) The significant changes in the pathway.

Figure 4

The GO terms interaction network of “2” sub-network (A). GO terms displayed as an interaction network using Cytoscape plug-in BinGO. Yellow nodes: nodes with P-value < 0.01 and Benjamini corrected P-value < .01. Biological function and pathway analysis of “2” sub-network. (B) The significant changes in the GO biological process. (C) The significant changes in the pathway.

Figure 5

The GO terms interaction network of “3” sub-network (A). GO terms displayed as an interaction network using Cytoscape plug-in BinGO. Yellow nodes: nodes with P-value < .01 and Benjamini corrected P-value < .01. Biological function and pathway analysis of “3” sub-network. (B) The significant changes in the GO biological process. (C) The significant changes in the pathway.

Discussion

OP is a systemic skeletal disease with characters of low bone mass, impaired bone microstructure and increase skeletal fragility, with a subsequent increase incidence of fractures. OP results from imbalance of bone formation and resorption mediated by interactions of several genetic, epigenetic and environmental factors. Over the past few years, great efforts have been made to explore the molecular mechanisms of OP. Numbers of genes/proteins have been identified to involved in OP, but it is still far from complete to thoroughly understand the molecular mechanism of OP. Therefore, it is necessary to decode the pathological mechanisms of OP at systems biology level. In current study, we first collected the OP-related genes, explored the genes’ interaction in system using network analyses, and showed a systematic framework to delineate the biochemical processed involving in OP.

GO analysis and pathway analysis have been used to assess biological functions enriched among OPgset genes.[ Our GO analysis indicated that OPgset genes enriched in biological regulation, metabolic process, developmental process and cell proliferation. GO terms, such as skeletal system development, ossification and osteoblast differentiation, were most significantly enriched in OPgset genes, suggesting the importance of these terms in the pathologic processes of OP. In addition, we found that the GO terms involved in MAPK cascade, blood vessel development, immune system development and apoptotic process were also in the most ten enriched GO terms, which consistent with the previous studies.[ At the same time, the pathway analysis showed that 57 pathways were enriched and involved in Wnt signaling pathway, osteoclast differentiation, cytokines and inflammatory response, steroid hormone biosynthesis and adipocytokine signaling pathway, all of which have been considered to play important roles in OP.[ Wnt can promote osteoblast differentiation and bone formation, while inhibit osteoclast differentiation and bone resorption.[ Inhibition of Wnt antagonists, such as sclerotin (SOST) and DKK1, stimulated bone formation and increase bone mineral density.[ Osteoclast differentiation, which contributes to pathological bone resorption, plays an important role in the pathology of OP, and is considered as a viable therapeutic target for OP.[ Cytokines and inflammatory response play key roles in regulating the bone regeneration and bone resorption.[ Typically, pro-inflammatory cytokines can suppress the bone formation and/or promote bone resorption activity[ while anti-inflammatory cytokines IL-10 and IL-13 showed the opposite effect.[ Steroid hormones, including glucocorticoids, androgens and estrogens, regulate the calcium and phosphorus homeostasis and bone mineralization via endocrine effects on bone, intestine, parathyroid glands and kidney.[ The imbalance between adipogenesis and osteogenesis of MSCs is involved in OP. A reduction in the osteogenesis of MSCs results in the impairment of bone formation and an increase of adipose tissue in bone marrow is found in OP.[

Three main modules were identified in pathway crosstalk analysis. The first module was mainly consisted with the pathways related to the bone resorption and formation. Among these pathways, osteoclast differentiation, MAPK signaling pathway, NF-kappa B signaling pathway, and HIF-1 signaling pathway, have been showed to be involved in bone tissue, as well as the development of osteogenesis.[ The second module included five pathways associated with Wnt signaling pathway, which indicated the important role of Wnt signaling pathway in OP. The third module was primarily dominated by the metabolic pathways of hormone or drug, including Drug metabolism-cytochrome P450, Glutathione metabolism, Steroid hormone biosynthesis, Ovarian steroidogenesis, Free Radical Induced Apoptosis, and Metabolism of xenobiotics by cytochrome P450. Furthermore, we calculated sharing frequency of OPgset genes in the network, and found that the most frequently shared genes included ligand for members of the frizzled family of seven transmembrane receptors (Wnt), TNF receptor superfamily member (TNFRSF11A, 11B, 1B), interleukins (IL-1), and Collagen type I (COL1A1, COL1A2), suggesting these genes might be more important in the pathogenesis of OP. The first and second modules connected through several edges, suggesting that these pathways/modules play a concerted role in pathogenesis of OP, rather than separate role.

According to the network theory, genes and proteins work in interconnected networks. We further constructed the PPI network of OPgset and screened the candidate key gene of OP. Hub nodes, which is high degree of connectivity to other nodes, are reported to be used as topological properties of the network to evaluate the importance of genes.[ Therefore, we evaluated node degree centrality, and also extracted the important sub-network via K-core analysis. In this study, we identified 143 notes from 217 OPgset genes/proteins and calculated degree centrality and K-coreness to identify hub nodes. We found that the genes/proteins including TGFB1, IL6, WNT3A, TNF, PTH, TP53, WNT1, IGF1, IL10, and SERPINE1, were in the most ten note degree. TGFB1 is a well-known molecular of the transforming grown factor B, which play key roles in bone development.[ The previous studies showed that TGF induces commitment to osteoblastic cell lineage, and prevents terminal osteoblastic differentiation.[ TGFB activated receptors form a complex with Smad4, then Smad4 translocated into the nucleus and interacted with RunX2, which was considered to be a key transcription factors of osteogenesis.[ IL-6 was reported to stimulate osteogenic differentiation of MSCs, and enhance extracellular matrix mineralization through increasing the expression of ALP, Runx2, and Collagen II.[ Infection-stimulated bone resorption can be inhibited by IL 10 in vivo and IL-10-deficient mice exhibited impaired bone formation, and bone fragility.[ TNF-α played a suppressed or promoted role in osteogenesis depending on its concentration, cell type, and exposure time.[ Wnt1 and Wnt3a are members of Wnt signaling pathway, which regulates bone formation.[ PTH is indicated to regulate calcium homeostasis via effects on kidney and bone.[ PTH was confirmed to directly activate survival in osteoblasts and increase osteoblast number in vitro and in vivo.[ IGF1 was reported to be produced by osteoblasts and regulated bone metabolism.[ The low serum IGF1 concentrations may reflect the low bone formation in humans.[

Aimed to extract the central modules, which are embedded in OPgset network, we analyze the K-core and construct k-core sub-networks of OPgset genes. Our data showed that there are three sub-networks from whole networks. The first sub-network is consisted with PTH, CALCA, GIPR, PTHLH, FSHR, ADRB3, PTH1R, and CALCR. Based on the results of GO and pathway analysis, we found that 156 GO terms and two pathways were enriched. The second sub-network is consisted with TGFB1, IGF1, SERPINE1, CFD, SPARC, MMRN1, ACTN1, and AHSG. We found that 114 GO terms and 12 pathways were enriched. The third sub-network was consisted with SFRP1, WNT16, SFRP4, WNT1, WNT4, WNT5B, WNT3A, LRP6, WLS, DKK1, and SOST. We found that 219 GO terms and 11 pathways were enriched. It is especially mentioned that the Wnt-related pathways were enriched in pathway analysis of whole network and sub-networks. The Wnt signaling pathway, is indicated to play important roles in bone formation, skeletal development and adult skeletal homeostasis.[ Recently, accumulated evidence suggested that the Wnt signaling pathway also plays central roles in the development of OP, possibly through inflammatory regulation, bone resorption, bone remodeling and joint destruction.[

In this study, based on the network theory, we reconstructed a gene network using OP associated genes compiled from selective literatures deposited in PUBMED, which for the first time enables analysis of the genes related to OP at a systematic level. Go analysis revealed biological processes related to skeletal system development, ossification and osteoblast differentiation were most significantly enriched in OPgset. Pathways analysis showed that 57 pathways were enriched and primarily involved Wnt signaling pathway, osteoclast differentiation, cytokines and inflammatory response, steroid hormone biosynthesis and adipocytokine signaling pathway, all of which have been shown to play important roles in OP. Besides, pathway crosstalk analysis indicated three major modules, with first module consisted of pathways mainly involved in bone development-related signaling pathways, second module in Wnt-related signaling pathway and third module in metabolic pathways. In addition, we calculated degree centrality to identify hub nodes, including TGFB1, IL6, WNT3A, TNF, PTH, TP53, WNT1, IGF1, IL10, and SERPINE1 in the most 10-note degree. We analyze the K-core and construct three k-core sub-networks of OPgset genes. Our study will improve our understanding of the pathogenesis of OP from the systems biological level and provide potential biomarkers for further study.

Author contributions

Conceptualization: Huijie Gu, Zhongyue Huang, Guangnan Chen, Xiaofan Yin.

Data curation: Huijie Gu.

Formal analysis: Huijie Gu, Zhongyue Huang, Guangnan Chen.

Funding acquisition: Xiaofan Yin.

Methodology: Huijie Gu, Zhongyue Huang, Guangnan Chen, Kaifeng Zhou, Yiming Zhang, Jiong Chen, Jun Xu.

Writing – original draft: Huijie Gu.

Writing – review & editing: Xiaofan Yin.