Transcriptome analyses identify hub genes and potential mechanisms in adenoid cystic carcinoma.

Adenoid cystic carcinoma (ACC) is one of the most frequent malignancies of salivary glands. The objective of this study was to identify key genes and potential mechanisms during ACC samples.The gene expression profiles of GSE88804 data set were downloaded from Gene Expression Omnibus. The GSE88804 data set contained 22 samples, including 15 ACC samples and 7 normal salivary gland tissues. The gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were constructed, and protein-protein interaction network of differentially expressed genes (DEGs) was performed by Cytoscape. The top 10 hub genes were analyzed based on Gene Expression Profiling Interactive Analysis. Then, DEGs between ACC samples and normal salivary gland samples were analyzed by gene set enrichment analysis. Furthermore, miRTarBase and Cytoscape were used for visualization of miRNA-mRNA regulatory network. KEGG pathway analysis was undertaken using DIANA-miRPath v3.0.In total, 382 DEGs were identified, including 119 upregulated genes and 263 downregulated genes. GO analysis showed that DEGs were mainly enriched in extracellular matrix organization, extracellular matrix, and calcium ion binding. KEGG pathway analysis showed that DEGs were mainly enriched in p53 signaling pathway and salivary secretion. Expression analysis and survival analysis showed that ANLN, CCNB2, CDK1, CENPF, DTL, KIF11, and TOP2A are all highly expressed, which all may be related to poor overall survival. Predicted miRNAs of 7 hub DEGs mainly enriched in proteoglycans in cancer and pathways in cancer.This study indicated that identified DEGs and hub genes might promote our understanding of molecular mechanisms, which might be used as molecular targets or diagnostic biomarkers for ACC.

Introduction

Adenoid cystic carcinoma (ACC) is one of the most frequent malignancies of the minor and major salivary glands and has poor long-term prognosis.[ ACC displays heterogeneous morphology because of their slow growth and tendency for perineural invasion, which makes it difficult to be diagnosed and characterized. After primary tumor resection, ACC can recur loco-regionally or with distant metastases in decades, which would require the long-term surveillance of all patients with ACC. Due to the resistance of ACC to chemotherapy or radiation therapy, nonresectable cases would be usually fatal.[ Therefore, the understanding of the molecular mechanism involved in proliferation, apoptosis, and invasion of ACC would be extraordinarily important for more effective diagnostic and therapeutic strategies.[

Microarrays are increasingly valued as a promising tool with great clinical applications in medical oncology: from molecular diagnosis to molecular classification of tumors, from new drug targets discovery to tumor response prediction, from patients’ stratification to prognosis prediction, and so on.[ The gene expression profiling study on ACC samples has been performed using microarray technology, which showed differentially expressed genes (DEGs) involved in different pathways, biological processes (BPs), or molecular functions (MFs). Now, microarray technology made it able to analyze the expression changes of mRNA comprehensively in the development and progression of ACC. Andersson et al[ collected tissue samples and investigated differences in gene expression between ACC and NSG. However, the interactions among the DEGs remain to be elucidated.

In this study, we downloaded the original data (GSE88804) from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/).[ Gene expression profiles of tumor cells in patients with ACC were compared with those in normal salivary gland (NSG) to identify DEGs. Whereafter, the DEGs were screened using R and Morpheus, followed by gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, integration of protein–protein interaction (PPI) network, module analysis, expression analysis, and survival analysis based on The Gene Expression Profiling Interactive Analysis (GEPIA), and gene set enrichment analysis (GSEA). Furthermore, miRTarBase and Cytoscape v3.6.0 were used for visualization of the miRNA-mRNA regulatory network. KEGG pathway analysis for predicted miRNAs was undertaken using DIANA-miRPath v3.0. By the method of analyzing the biological functions and pathways, we may explore the potential biomarkers for diagnosis, drug targets, and prognosis of ACC.

Methods

Microarray data

The gene expression profiles of GSE88804 were downloaded from the GEO database, which were the only available data set of ACC samples. And the original authors neither reported DEGs data nor performed bioinformatic analysis based on this data set. GSE88804, which was based on Affymetrix GPL6244 platform (Affymetrix Human Gene 1.0 ST Array), was submitted by Andersson et al.[ The GSE88804 data set contained 22 samples, including 15 ACC samples (13 surgical samples of ACC and 2 ACC xenografts) and 7 NSG tissues.

Differential expression analysis and identification of DEGs

Gene microarray analyses were all conducted through R software (version 3.5.1, https://www.r-project.org/; The R Foundation). Raw CEL data were imported into R, and we performed relative log expression (RLE) plots for detecting and visualizing the unwanted variation in high dimensional microarray data among all tissue samples through affyPLM package and RcolorBrewer package. DEGs between ACC samples and NSGs tissues were identified through GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/) with the cut-off criterion (P-value < .05 and Log|FC| > 2).[ Then, we clicked “Save all results” and the results of DEGs were shown after entering into GSE88804 in GEO accession. Then, the heat map of the DEG expression (top 50 upregulated and downregulated genes) was carried out by Morpheus (https://software.broadinstitute.org/morpheus/).

GO and KEGG pathway enrichment analysis of DEGs

The GO analysis is a common effective method for annotating genes and identifying characteristic biological attributes.[ KEGG (http://www.genome.jp/) is a knowledge database for a systematic analysis of gene functions.[ Comprehensively, mapping of the user's gene to the related biological annotation in the Database for Annotation, Visualization and Integrated Discovery (DAVID) database (https://david.ncifcrf.gov/) is an essential foundation for the success of any high-throughput gene functional analysis.[ To analyze the DEGs, GO enrichment and KEGG pathway analysis were respectively performed using the DAVID online tool. P < .05 was considered statistically of significance.

PPI network and module analysis

Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org) database is the online tool, which is designed to evaluate the PPI information.[ STRING (version 10.5) covers 9.6 million proteins from 2031 organisms. To evaluate the interactive relationships among the DEGs, we mapped all the DEGs to STRING, and only validated interactions with a combined score >0.4 were considered as significant. Then, the PPI network was constructed using the Cytoscape software (version 3.6.0, https://cytoscape.org/).[ The plug-in cytoHubba was used to select top 10 hub genes, while the plug-in Molecular Complex Detection (MCODE) was used to screen the modules of the PPI network in Cytoscape. The criteria were set as follows: MCODE scores ≥4 and number of nodes ≥4. Moreover, KEGG pathway analyses were performed for DEGs in these modules. P < .05 was also considered to be significant.

Expression analysis and survival analysis

The GEPIA is a new outstanding interactive web tool for analyzing the RNA-seq expression data of 9736 tumors and 8587 normal samples from The Cancer Genome Atlas and the Genotype-Tissue Expression projects.[ GEPIA can provide customizable functions such as tumor/normal differential expression analysis, patient survival analysis, profiling according to cancer types or pathological stages, correlation analysis, and so on.

Gene set enrichment analysis

The GSEA was also applied to identify the significant pathways in GSE88804 based on GO-BP and KEGG pathway. The coefficients of the Spearman correlation were defined as the weight of genes between genes and sample label.[ Statistical significance was assessed with the enrichment score of enrichment results, which generated from 1000 random permutations of the gene sets to obtain P values. The pathways with levels of False Discovery Rate (FDR) <25% and P < .01 were considered to be significant.

Construction of the miRNA-mRNA regulatory network and identification of miRNA-associated pathways

The miRNAs, a class of noncoding RNA with 20 to 22 nucleotides, can bind to the 3′Untranslated Regions of targeted mRNAs to induce translational repression or degradation of mRNAs.[ In our study, miRNAs interacting with hub mRNAs were predicted using an experimentally validated microRNA-target interactions database (miRTarBase).[ Cytoscape was used for the construction of the miRNA-mRNA regulatory network. We also performed KEGG pathway analysis for predicted miRNAs based on DIANA-miRPath v3.0 is a useful web tool which can provide experimentally supported miRNAs-mRNA interaction.[ The results of KEGG enrichment for predicted miRNAs were visualized using package ggplot2 in R.

Results

Supplementary Figure 1, illustrated the RLE among all of the samples after normalization. Based on data preprocessing and Student t test, a total of 20,329 genes were identified. Based on the criteria of P < .05 and Log|FC| > 2, we identified a total of 382 DEGs in ACC samples compared with NSG, which were shown in Supplementary Table 1. About 119 DEGs were found to be upregulated in ACC, while 263 genes were downregulated (Fig. 1). DEGs expression heat map (top 50 upregulated and downregulated genes) are shown in Figure 2.

Figure 1

Volcano plot of differentially expressed genes. Red dots = upregulated genes, blue dots = downregulated genes.

Figure 2

Heat map of the top 100 differentially expressed genes (50 upregulated genes and 50 downregulated genes). Red = upregulation; blue = downregulation.

GO term enrichment analysis

We uploaded all the DEGs to the online software DAVID to identify overrepresented GO categories and KEGG pathways. GO analysis results showed that upregulated DEGs were significantly enriched in BPs, including extracellular matrix organization, cell adhesion, mitotic nuclear division, cell division, and skeletal system development (Table 1 and Supplementary Figure 2); the downregulated DEGs were significantly enriched in BPs, including retina homeostasis, ethanol oxidation, detection of chemical stimulus involved in sensory perception of bitter taste, biomineral tissue development, and transmembrane transport (Table 1 and Supplementary Figure 2). For MF, the upregulated DEGs were enriched in calcium ion binding, glycosaminoglycan binding, extracellular matrix structural constituent, cyclin-dependent protein serine/threonine kinase activity, and chromatin binding, and the downregulated DEGs were enriched in extracellular exosome, extracellular space, extracellular region, microvillus, and endoplasmic reticulum (Table 1 and Supplementary Figure 2). In addition, GO cell component (CC) analysis also displayed that the upregulated DEGs were significantly enriched in extracellular matrix, proteinaceous extracellular matrix, plasma membrane, spindle microtubule, and spindle, and downregulated DEGs enriched transporter activity, aldehyde dehydrogenase (NAD) activity, oxidoreductase activity (acting on the aldehyde or oxo group of donors, NAD or NADP as acceptor), protein homodimerization activity and alcohol dehydrogenase activity, zinc dependent (Table 1 and Supplementary Figure 2,).

Table 1

Gene ontology analysis of differentially expressed genes associated with adenoid cystic carcinoma.

KEGG pathway analysis

Table 2 and Supplementary Figure 3, contain the most significantly enriched pathways of the upregulated DEGs and downregulated DEGs analyzed by KEGG analysis. The upregulated DEGs were enriched in p53 signaling pathway, glycosphingolipid biosynthesis-lacto and neolacto series, while the downregulated DEGs were enriched in salivary secretion, tyrosine metabolism, peroxisome proliferator-activated receptor (PPAR) signaling pathway, fatty acid degradation, regulation of lipolysis in adipocytes, glycolysis/gluconeogenesis, histidine metabolism, arginine and proline metabolism, drug metabolism-cytochrome P450, glycine, serine and threonine metabolism, phenylalanine metabolism, gastric acid secretion, ABC transporters, metabolic pathways and AMP-activated protein kinase signaling pathway.

Table 2

Kyoto Encyclopedia of Genes and Genomes pathway analysis of differentially expressed genes associated with adenoid cystic carcinoma.

Module screening from the PPI network

Based on the information in the STRING database, the top 10 hub nodes with higher degrees were screened using plug-ins CytoHubba through Cytoscape. These hub genes included TOP2A, CDK1, KIF11, BUB1B, CCNB2, DTL, KIF23, ANLN, CENPF, and NUSAP1. Among these genes, TOP2A showed the highest node degree, which was 39. Moreover, a total of 226 nodes and 519 edges were analyzed using plug-ins MCODE. The top 4 significant modules were selected, and the functional annotation of the genes involved in the modules was analyzed (Fig. 3). Enrichment analysis showed that the genes in modules 1 to 4 were mainly associated with cell cycle, p53 signaling pathway, PPAR signaling pathway, tyrosine metabolism, drug metabolism-cytochrome P450, histidine metabolism, renin secretion, and morphine addiction.

Figure 3

Top 4 modules from the protein–protein interaction network. (A) Module 1. (B) The enriched pathways of module 1. (C) Module 2. (D) The enriched pathways of module 2. (E) Module 3. (F) The enriched pathways of module 3. (G) Module 4. (H) The enriched pathways of module 4.

Expression analysis and survival analysis based on GEPIA

We applied GEPIA to validate gene expression level and survival rates of the TOP 10 hub genes between ACC tissues and normal tissues, and 7 genes significantly increased expression levels with obvious changes of survival analysis in ACC tissues. Then, box plots of expression and corresponding survival plots were conducted based on GEPIA (Fig. 4A–G).

Figure 4

Expression analysis and survival analysis based on Gene Expression Profiling Interactive Analysis (GEPIA). (A) ANLN. (B) CCNB2. (C) CDK1. (D) CENPF. (E) DTL. (F) KIF11. (G) TOP2A.

The DEGs between ACC samples and NSG samples were also analyzed by the GSEA method which used a database of several thousand predefined sets of genes. GSEA is able to detect small and significant expression changes in these connected genes that cannot be revealed by gene-by-gene comparisons. Then, the results of GSEA showed that 1363 gene sets are upregulated, in which 898 gene sets are significant at FDR < 0.25 and 354 gene sets are significantly enriched at nominal P-value <.01; and 2305 gene sets are downregulated, in which 768 gene sets are significant at FDR <0.25 and 763 gene sets are significantly enriched at nominal P-value <.01. The top 3 upregulated and downregulated GO and KEGG pathways are listed in Figure 5.

Figure 5

The top 3 representative upregulated and downregulated enrichment plots of GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways for adenoid cystic carcinoma (ACC) were analyzed by gene set enrichment analysis (GSEA). (A–C) Upregulated. (D–G) Downregulated.

MiRNA-target regulatory network

The miRNAs binding to DEGs in subnetworks were predicted using miRTarBase. The miRNA-mRNA regulatory network included 82 nodes and 111 edges. Of these, hsa-miR-192-5p and hsa-miR-215-5p can antagonize ANLN, CENPF, and DTL, while hsa-miR-193b-3p can antagonize CDK1, KIF11, and TOP2A (Fig. 6A). Also, we performed KEGG pathway enrichment analysis of these predicted miRNAs, which mainly enriched in proteoglycans in cancer, pathways in cancer, fatty acid metabolism, hippo signaling pathway, and Transforming growth factor-β (TGF-β) signaling pathway (Fig. 6B and Supplementary Table 3).

Figure 6

MiRNA-mRNA regulatory networks in adenoid cystic carcinoma (ACC). (A) MiRNA-mRNA regulatory networks. (B) Bubble graph for TOP 20 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of predicted miRNAs.

Discussion

The ACC is a product of somatic, cumulative genetic, epigenetic, and endocrine aberrations.[ The relative rarity and slow growing of ACC aggressive nature have complicated the molecular markers. The understanding of the molecular mechanism of ACC is of critical importance for its diagnosis and treatment. Because microarray and high-throughput sequencing provide expression levels of thousands of genes in the human genome, they have been widely applied to predict the potential therapeutic targets for ACC.[ Owing to the original paper (PMID: 28954282) for their contribution to the microarray data, we can conduct this study.[ Gao et al[ identified a unique ACC signature with parallel MYB-dependent and MYB-independent biomarkers and identified VCAN/HAPLN1 complexes as a potential target, which showed that forced MYB-NFIB expression in NSG cells alters cell adhesion and cell morphology in vitro and depletion of VCAN blocked tumor cell growth of ACC tumor. Rettig et al's research proved that NFIB was a vital role in ACC oncogenesis.[ Mitani et al[ conducted whole-genome sequencing in 21 salivary ACCs, defining novel molecular subclasses characterized by MYBL1 rearrangements and 5′-NFIB gene fusions. Brayer et al[ suggested that proteins of MYB and MYBL1 were oncogenic targets in ACC. Ho et al[ also observed MYB-NFIB translocations and somatic mutations in MYB-associated genes, suggesting these aberrations as critical events. Bell et al[ implied that EN1, DLX6, and OTX1 may be potential drivers of ACC. Andersson et al[ indicated that the MYB-NFIB fusion drives ACC cells’ proliferation, which is regulated through AKT-dependent signaling induced by IGF1R overexpression. In a recent research, Frerich et al[ performed detailed RNA-sequencing (RNA-seq) analysis on 68 ACC tumor samples, which resulted that MYB or MYBL1 would be direct targets of Myb proteins in ACC tumors.

In our study, we extracted the data from GSE88804 and identified 119 upregulated and 263 downregulated DEGs between ACC and NSG using bioinformatics analysis. Function annotation showed that these DEGs were mainly involved in extracellular exosome, extracellular region, extracellular space, and salivary secretion. Cumulative evidence has detailedly demonstrated that coexpression gene normally consists of a group of genes with similar expression profiles, which participate in the parallel BP as well. To better understand the interactions of the DEGs, we performed GO and KEGG pathway analysis.

The GO term analysis showed that upregulated DEGs were mainly involved in extracellular matrix organization, extracellular matrix, cell adhesion, mitotic nuclear division, and calcium ion binding, while downregulated DEGs were involved in extracellular exosome, extracellular space, and extracellular region. Furthermore, the enriched KEGG pathways of upregulated DEGs included p53 signaling pathway and glycosphingolipid biosynthesis-lacto, and neolacto series. A previous study has shown that related upregulated genes of the p53 signaling pathway in human development could predict the overall survival of patients with ACC.[ Recent evidence indicated that p53 signaling pathway might be associated with ACC metastasis and progression.[ Downregulated DEGs were related to salivary secretion, tyrosine metabolism, and PPAR signaling pathway. Porto-Figueira et al[ reported that butanoate metabolism and tyrosine metabolism might be highly activated in cancers, as well as tyrosine metabolism in a lesser extent. Antonosante et al[ highlighted the different roles of PPAR isotypes in various cancer cells.

By constructing the PPI, we identified 10 hub genes that can provide new ideas for the therapeutic studies in ACC. The top 10 hub genes were TOP2A, CDK1, KIF11, BUB1B, CCNB2, DTL, KIF23, ANLN, CENPF, and NUSAP1. TOP2A was identified as one of the hub genes that exhibit the highest degree of connectivity. TOP2A, as a protein-coding gene, might promote the development of the tumor, especially in proliferation and differentiation.[ Ren et al[ reported that high TOP2A expression suggested that the more significant relationship with worse prognosis of cancer. The 2nd hub gene CDK1, one of the Ser/Thr protein kinase family, is markedly related to transferase activity, transferring phosphorus-containing groups, and protein tyrosine kinase activity. Chu et al[ demonstrated the treatment of human lung high metastasis cell line of ACC cells with 5 to 20 μM sulforaphane resulted in G(2)/M cell cycle arrest. The third hub gene KIF11 is one of the kinesin-like protein family, and the function of this gene product includes centrosome separation, chromosome positioning, and establishing a bipolar spindle during cell mitosis. A recent study has proposed that KIF11 upregulation represented an independent prognostic indicator for the survival of patients with cancer and it might be a therapeutic target for cancers.[ BUB1B encodes a kinase, which is involved in spindle checkpoint function. Lee et al[ developed evidence that BUB1B might offer a predictive marker for aggressiveness and drug response. Another hub gene, CCNB2, is a member of the cyclin family. Qian et al[ reported that overexpression of CCNB2 protein may be associated with clinical progression, which may lead to poor prognosis. The other 5 hub genes are DTL, KIF23, ANLN, CENPF, and NUSAP1, respectively. Expression analysis and survival analysis based on GEPIA showed that the top 10 genes may be related to the poor overall survival rate, but only ANLN, CCNB2, CDK1, CENPF, DTL, KIF11, and TOP2A are highly expressed in ACC tumor samples. Chen et al[ indicated that DTL is a potential novel target gene for the treatment of cancers. Vikberg et al[ showed that the level of KIF23 could be elevated due to the additional copy of chromosome 15 demonstrated. Recently, Long et al[ demonstrated that ANLN may be involved in developmental processes through the regulation of nuclear division pathway. Moreover, CENPF encodes the protein that associates with the centromere–kinetochore complex. Aytes et al[ reported that coexpression of FOXM1 and CENPF would be a robust prognostic indicator of poor survival and metastasis. Overexpression of NUSAP1 might impact prostate cancer progression by increasing proliferation or invasion of cancer cells.[ In addition, some groups have also published their data about RNA-seq of ACC samples and normal control samples. Feng et al[ found that miR-155 may be involved in ACC metastasis through UBA2-related pathways, and UBA2 may act as a mediator of ACC metastasis. Han et al[ revealed that SCUBE3 may play an essential role in epithelial mesenchymal transition of ACC. Chen et al[ concluded that TGF-β1 might play a important role during lung metastasis of ACC and be considered as a candidate target for treatment of metastatic ACC. Kasamatsu et al[ showed that CACNA1C, CCND2, COL18A1, CTNNB1, HES1, ITGA9, MLL, NCAM1, PLXNB1, PTTG1, SCARB1, SEMA3F, TIAM1, and TP73L were cancer-related genes of ACC.

Module analysis of the PPI network revealed that the development of ACC was associated with renin secretion, PPAR signaling pathway, and salivary secretion. The aspartyl-protease renin is the key regulator of the renin-angiotensin-aldosterone system, which is critically involved in extracellular fluid volume and blood pressure homeostasis of the body. Renin secretion was still considered to be one of the main functional circuits of miR-210.[ PPARs are nuclear hormone receptors which are activated by fatty acids and their derivatives. Salivary secretion always occurs in response to stimulation by neurotransmitters.[

Many differential expression of miRNAs, such as miR-130a, miR-320a, and miR-21, were associated with ACC development and progression.[ In the present study, we identified 75 miRNAs based on DEGs from subnetworks and these miRNAs mainly enriched in proteoglycans in cancer, pathways in cancer, fatty acid metabolism, hippo signaling pathway, and TGF-β signaling pathway. Also, miRNA-mRNA regulatory network analysis indicated that miR-192-5p, miR-193b-3p, and miR-215-5p may act as potential key miRNAs to regulate corresponding mRNAs. Actually, previous studies did not verify these predicted miRNAs and related pathways, but some experts verified MYB/miR-130a activated the STAT3 and AKT pathways by downregulating NDRG2, while some experts indicated that downregulation of miR-125a-5p promotes ACC progression through p38 signal pathway and miR-125a-5p can be a potential therapeutic target of ACC.[ Therefore, these predicted miRNAs and signaling pathways can be served as potential diagnostic biomarkers and therapeutic targets for ACC, which may provide potential hallmarks for further experimental studies.

In this study, we are the first to observe that ACC was associated with farmer breast cancer cluster 2, Kauffmann DNA repair genes, Mori immature B lymphocyte double negative, glycolysis and gluconeogenesis, citrate cycle/TCA cycle, and pentose phosphate pathway by GSEA, but the association ACC with these pathways and BPs need to be validated.

However, there are also several limitations in this study. First, only 15 ACC samples and 7 NSG samples were included. Second, owing to the relative rarity of ACC and its slow growing of ACC tissue, we could not collect enough ACC tissues to validate the results based on the bioinformatic analysis. Third, our study was concluded from the bioinformatics analysis of the expression profile data, which was downloaded from the GEO database. So, further experiments will be needed to verify our results.

Conclusion

Our work provides a comprehensive bioinformatics analysis of the DEGs, which may be involved in the progress of ACC. The study provides a set of potential targets for future investigation. However, further molecular biological experiments would be required to confirm the function and pathways of the DEGs in ACC, with the goal of improving treatment response and patient outcome.

Ethics

Ethics approval and patient written informed consent were not required because all analyses in our study were performed based on data from the GEO database.

Acknowledgment

The authors are grateful to Dr Hui-Zi Li for his statistical support for the study.

Author contributions

Conceptualization: Hong-Bing Liu, Meng-Si Luo.

Data curation: Guan-Jiang Huang, Meng-Si Luo.

Formal analysis: Hong-Bing Liu, Meng-Si Luo.

Investigation: Guan-Jiang Huang, Meng-Si Luo.

Methodology: Hong-Bing Liu, Guan-Jiang Huang, Meng-Si Luo.

Project administration: Hong-Bing Liu.

Resources: Hong-Bing Liu, Guan-Jiang Huang.

Software: Guan-Jiang Huang, Meng-Si Luo.

Supervision: Hong-Bing Liu.

Validation: Guan-Jiang Huang.

Visualization: Guan-Jiang Huang, Meng-Si Luo.

Writing – original draft: Guan-Jiang Huang, Meng-Si Luo.

Writing – review & editing: Hong-Bing Liu, Meng-Si Luo.