Identification of functional genes regulating gastric cancer progression using integrated bioinformatics analysis

Abstract

BACKGROUND

Gastric cancer (GC) is one of the most common cancers and has a poor prognosis. Treatment of GC has remained unchanged over the past few years.

AIM

To investigate the potential therapeutic targets and related regulatory biomarkers of GC.

METHODS

We obtained the public GC transcriptome sequencing dataset from the Gene Expression Omnibus database. The datasets contained 348 GC tissues and 141 healthy tissues. In total, 251 differentially expressed genes (DEGs) were identified, including 187 down-regulated genes and 64 up-regulated genes. The DEGs’ enriched functions and pathways include Progesterone-mediated oocyte maturation, cell cycle, and oocyte meiosis, Hepatitis B, and the Hippo signaling pathway. Survival analysis showed that BUB1, MAD2L1, CCNA2, CCNB1, and BIRC5 may be associated with regulation of the cell cycle phase mitotic spindle checkpoint pathway. We selected 26 regulated genes with the aid of the protein-protein interaction network analyzed by Molecular Complex Detection.

RESULTS

We focused on three critical genes, which were highly expressed in GC, but negatively related to patient survival. Furthermore, we found that knockdown of BIRC5, TRIP13 or UBE2C significantly inhibited cell proliferation and induced cell apoptosis. In addition, knockdown of BIRC5, TRIP13 or UBE2C increased cellular sensitivity to cisplatin.

CONCLUSION

Our study identified significantly upregulated genes in GC with a poor prognosis using integrated bioinformatics methods.

Key Words: Gastric cancer, Bioinformatics analysis, Differentially expressed gene, Protein-protein interaction network, Cisplatin resistance

Core Tip: Gastric cancer (GC) is one of the most common malignancies of the digestive system with few genetic markers for its early detection and prevention. 348 GC tissues and 141 normal tissues were analyzed in this study; 251 differentially expressed genes (DEGs) were identified including 187 down-regulated genes and 64 up-regulated genes. Significantly upregulated genes with a poor prognosis in GC were detected using integrated bioinformatics methods. Furthermore, the molecular mechanism of GC via bioinformatics analysis showed three DEGs (BIRC5, TRIP13, or UBE2C) which play key roles in the progression of CG.

INTRODUCTION

Gastric cancer (GC) remains the most common cause of cancer death globally, although the mortality of GC has markedly declined in the past year[1]. Although tremendous efforts have been made to analyze the pathogenesis, there is limited knowledge on the exact molecular mechanisms involved in GC. It was reported that several prognostic potential biomarkers have been examined, but the precise pathogenic mechanism of GC remains to be illustrated[2,3]. It is essential to investigate the variations in global gene expression to enhance the efficacy of treatment and to assist us in understanding the pathogenesis of GC. As a mature technology, gene chip is widely used in modern medicine[4]. It can detect differentially expressed genes (DEGs) quickly, which are stored in public databases. Therefore, many related studies on GC have been published in recent years[5-8]. Genetic factors were revealed to play an important role in susceptibility to GC.

The Gene Expression Omnibus (GEO) datasets GSE13911, GSE66229, and GSE79973 were selected for this study. The GEO2R web tool and Venn diagram software were used to assess the DEGs in the three datasets mentioned above. To examine the screened DEGs with Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, we used the Database for Annotation, Visualization, and Integrated Discovery (DAVID, http://david.abcc.ncifcrf.gov). Then, we analyzed the DEGs and screened the protein-protein interaction (PPI) network modules in GC using the PPI network DEGs.

We applied Cytotype Molecular Complex Detection (MCODE) to screen the key gene modules and co-expression networks to identify their relative significance in the three datasets above. Among these DEGs, 26 genes showed a significant difference in the PPI network. Subsequently, the 26 core DEGs were generated by the Kaplan Meier plotter online database (http://www.kmplot.com) for crucial prognostic information (P < 0.05). Also, we determined the expression of the DEGs in Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn) (P < 0.05) between GC tissues and normal gastric tissues.

We found 22 DEGs based on these data. The 22 DEGs were then re-analyzed for KEGG pathway enrichment. Subsequently, based on survival analysis, three DEGs (BIRC5, TRIP13, and UBE2C) were important related genes and played key roles in cellular pathways. Thus, we identified the underlying gene pathways and potential molecular mechanisms. These results may contribute to the treatment, diagnosis, and prevention of GC. Cellular experiments also revealed that targeting these three genes inhibited proliferation of GC cell lines and overcame resistance to cisplatin, which may provide novel therapeutic targets for GC.

MATERIALS AND METHODS

Microarray data information

The expression profiles of GC and healthy gastric tissues in the GSE13911, GSE66229, and GSE79973 datasets were available from the NCBI-GEO (http://www.ncbi.nlm.nih.gov/geo/), which were all performed on the GPL570 platform, an open database of microarray/gene profiles. Microarray data of GSE13911, GSE66229, and GSE79973 were all based on Affymetrix HGU 133 Plus 2 microarrays (Affymetrix, Inc., Santa Clara, CA, United States).

In GSE13911, there were 69 samples from 38 GC patients and 31 normal individuals; GSE66229 contains microarray data from healthy controls and GC patients. There were 400 samples in GSE66229, including 300 GC samples and 100 control samples. GSE79973 contains data from 10 cognitively healthy controls and 10 GC patients.

Data processing of DEGs

DEGs in GC tissues and healthy tissues were identified via GEO2R online tools[9] with an adjusted P value < 0.05 and |log FC| > 0.5. In the three datasets above, the original results were selected by Venn diagram software online to determine the common DEGs. The DEGs were distinguished by log FC as follows: the genes with log FC (fold change) > 0 were up-regulated, while the genes with log FC < 0 were down-regulated.

Gene ontology and pathway enrichment analysis

Gene ontology (GO; http://www.geneontology.org) analysis is a commonly used approach for the unification of biology that offers information on gene product function using ontologies to represent biological knowledge[10]. The KEGG is a large-scale database[11] including drugs, genomes, diseases, chemical materials, biological pathways, and so on. The DAVID is a useful bioinformatics tool designed to check a large number of gene or protein functions online[12]. We used DAVID to identify GO categories and pathways (P < 0.05).

PPI network and module analysis

Based on the Search Tool for the Retrieval of Interacting Genes (STRING) information database[13], the PPI network was constructed. The STRING database was then visualized using Cytoscape software version, a popular open-source software tool from the network: www.Cytoscape.org[14], to identify the prognostic correlation among the DEGs above, we set the confidence score to ³0.4 (modest confidence). The MCODE program in Cytoscape was also used to explore the PPI network modules (we set: Node score cut off = 0.2, degree cut off = 2, k-core = 2, and max. depth from seed = 100).

Survival analysis and RNA sequencing expression of core genes

Kaplan Meier-plotter (KM plotter, http://kmplot.com/analysis/) was used to assess the effect of numerous genes on survival based on EGA, TCGA database, and GEO. The log-rank P value and hazard ratio with 95% confidence intervals (CI) were computed and are shown on the plot[15]. To validate these DEGs, we used the GEPIA website to analyze RNA sequencing expression data based on thousands of samples from the GTEx projects and TCGA[16].

Patient samples and RT-PCR detection

We performed RT-PCR analysis of 15 pairs of patient samples including tumor and adjacent normal tissues in accordance with the ethics permission approved by our hospital. RNA was extracted using Trizol methods as previously described. The cDNA was the synthesized with a reverse transcription kit (Takara, China). The expression levels of BIRC5, TRIP13, and UBE2C genes were detected with a SYBR real-time PCR kit from Takara. The primers used for detecting these three genes were purchased from Thermo Fisher Co.; GAPDH was used as a control.

Cells and transfection

Human GC cell lines MKN-28 were purchased from the American Type Culture Collection. All cells were maintained at 37^oC in a 5% CO₂ humidity chamber, in RPMI-1640 medium supplemented with 10% fetal bovine serum (FBS, Gibco) and streptomycin and penicillin (10 mg/mL).

CCK-8 assay and apoptosis assay

To determine cell proliferation, we used a CCK-8 assay and cells with specific gene knockdown. Briefly, 5000 cells were seeded in each well of 96 well plates. At 24 h, 48 h, and 72 h after plating, the cells were stained with CCK8 solution from a CCK-8 kit (Beyotime, China) and detected under OD570 wavelength. Cell apoptosis was performed using an Annexin V/PI double staining kit according to the manufacturer’s instructions.

Colony formation assay

For the colony formation assay, cells were seeded in 60 mm dishes at 400, 800, 1600, 3200 cells/dish. Twenty-four hours later, cells were treated with cisplatin at concentrations of 0 μM, 2 μM, 4 μM, and 8 μM. After treatment, the cells were returned to the incubator for an additional 10 d. Colonies were counted after staining with ultraviolet staining solution.

Statistics analysis

As shown in Figure 1, the flow diagram of this study demonstrates that the raw experimental data were processed using several software packages. The moderate t-test was applied to identify DEGs, and Fisher’s Exact test was used to analyze GO and KEGG annotation enrichments[17].

Figure 1 A flow diagram of the present study. GEO: Gene Expression Omnibus; DEG: Differentially expressed gene; DAVID: Database for Annotation, Visualization, and Integrated Discovery; KEGG: Kyoto Encyclopedia of Genes and Genomes; STRING: Search Tool for the Retrieval of Interacting Genes; GEPIA: Gene Expression Profiling Interactive Analysis.

RESULTS

Identification of DEGs in GC

In the current study we included 348 GC tissues and 141 healthy tissues. With the help of GEO2R online tools, we obtained 5599, 2755, and 2473 DEGs from the three databases: GSE13911, GSE66229, and GSE79973, respectively. Among the three datasets, the common DEGs were identified via Venn diagram software. The findings revealed that 251 common DEGs were identified, including 187 and 64 genes which controlled down- and up-related genes in GC tissues, respectively (Supplementary Table 1; Figure 2).

Figure 2 Authentication of 251 common differentially expressed genes in the three datasets (GSE13911, GSE66229, and GSE 79973) through Venn Diagrams software. Different color represents different datasets. A: Sixty-four differentially expressed genes (DEGs) were upregulated in the three datasets (logFC > 0); B: 187 DEGs were downregulated in the three datasets (logFC < 0). DEGs: Differentially expressed genes. Available online: http://bioinformatics.psb.ugent.be/webtools/Venn/.

Identification and function of DEGs in GC

Associated with GO terms, the up-regulated and down-regulated DEGs were enriched in cell division, mitotic nuclear division, and mast cell cytokine production. The markedly up-regulated and down-regulated DEGs in biological process (BP) were involved in cell division, mitotic nuclear division, and mast cell cytokine production. In addition to the enrichment in mitotic metaphase plate congression, the up-regulated DEGs were significantly enriched in the BPs of cell division. In contrast, the down-regulated DEGs were significantly enriched in the metabolic process.

Gene ontology and pathway enrichment analysis

We performed functional and pathway enrichment analyses to investigate the biological classification of DEGs via DAVID. There were three categories in GO analysis: BP, cellular component (CC), and molecular function (MF) GO. In the BP-associated category shown by GO pathway enrichment analysis, the genes were significantly involved in nuclear chromosome segregation, sister chromatid segregation, cell division, cell cycle process, mitotic cell cycle process and down-regulated DEGs in organic acid catabolic process, carboxylic acid catabolic process, small molecule catabolic process, lipid modification, fatty acid catabolic process, and fatty acid oxidation (Supplementary Table 2). MF of DEGs were significantly enriched in identical protein binding, microtubule binding, enzyme binding, tubulin binding, cytoskeletal protein binding, protein kinase binding, cofactor binding, coenzyme binding, potassium channel regulator activity, ligand-gated ion channel activity, and ligand-gated channel activity (Supplementary Table 2). Changes in the CCs of DEGs were mainly enriched in the condensed chromosome kinetochore, condensed chromosome, centromeric region kinetochore, condensed chromosome, chromosome, centromeric region, and midbody (Supplementary Table 2; Figure 3).

Figure 3 Gene ontology terms and Kyoto Encyclopedia of Genes and Genomes pathways of differentially expressed genes significantly enriched in gastric cancer. GO: Gene ontology.

PPI network and modular analysis

A total of 251 DEGs containing 193 nodes and 633 edges were imported into the DEGs PPI network complex, and 187 down-regulated and 64 up-regulated genes were included (Figure 4A). The 187 down-regulated DEGs were not contained in the DEGs PPI network (Figure 4A). We then used Cytotype MCODE to gain further results: The outcomes showed that among the 64 nodes there were 26 central nodes of up-regulated genes (Figure 4B).

Figure 4 Common differentially expressed genes protein-protein interaction network constructed by Search Tool for the Retrieval of Interacting Genes online database and Module analysis. A: A total of 251 differentially expressed genes (DEGs) were found in the DEGs protein-protein interaction network complex. The nodes represent proteins; the edges represent the interaction of proteins; green circles represent down-regulated DEGs and red circles represent up-regulated DEGs; B: Module analysis via Cytoscape software. PPI: Protein-protein interaction; STRING: Search Tool for the Retrieval of Interacting Genes.

KM plotter and GEPIA database analysis

We identified the survival data related to 26 core genes via Kaplan Meier plotter (http://kmplot.com/analysis/). Among these genes, 24 genes showed a significantly poor survival, while 2 showed no significance (P < 0.05, Supplementary Table 3; Figure 5; Supplementary Figure 1). Then, GEPIA was used to determine the expression level of 24 genes between cancerous and normal tissues. In contrast to normal stomach samples, the results revealed that 22 of 24 genes were highly enriched in GC samples (P < 0.05, Supplementary Table 4; Figure 5; Supplementary Figure 2).

Figure 5 Analysis of BIRC5, TRIP13, and UBE2C genes in patient tissues and the correlation with survival. A-C: Overall survival of gastric cancer patients with high expression and low expression of BIRC5 (A), TRIP13 (B), and UBE2C (C) genes; D-E: The expression of these genes was analyzed with GEPIA tools. ^aP < 0.05.

KEGG pathway enrichment of 22 selected re-analyzed genes

To better understand the possible pathways of these 22 selected DEGs, we re-analyzed KEGG pathway enrichment via DAVID (P < 0.05). We found that five core genes (BUB1, MAD2L1, CCNA2, CCNB1, and BIRC5) were markedly enriched in the progesterone-mediated oocyte maturation and cell cycle pathway (P = 1.9E-3, Supplementary Table 5). Among these genes, the overexpression of BIRC5, TRIP13, and UBE2C was negatively correlated with poor outcomes in GC patients, and were further investigated.

Upregulation of BIRC5, TRIP13 and UBE2C in GC tissues from clinical patients is related to cell proliferation

To validate the expression level of these potential genes in cancer tissues, we collected 15 pairs of GC tissues and adjacent normal tissues. Following RT-PCR assay, we found that BIRC5, TRIP13, and UBE2C were significantly upregulated in tumor tissues (Figure 6A-C), and these data showed that these genes may be used as biomarkers for GC. To investigate whether these genes exert key functions in cancer, we transfected specific shRNAs into GC cells. Using the CCK-8 assay, we observed that cell proliferation was significantly suppressed in BIRC5, TRIP13, or UBE2C knockdown cells, respectively (Figure 6D-F).

Figure 6 Increased level of BIRC5, TRIP13, and UBE2C genes are related to cell proliferation of gastric cancer cells. A-C: Expression of BIRC5, TRIP13, and UBE2C genes in 15 pairs of tumor and adjacent normal tissues; D-E: Knockdown of BIRC5, TRIP13, or UBE2C genes significantly suppressed cancer cell proliferation. ^aP < 0.05; ^bP < 0.01.

Knockdown of BIRC5, TRIP13, or UBE2C reduced cisplatin resistance

Cisplatin is one of the most frequent drugs currently used to treat GC. However, cisplatin resistance is an obstacle for successful therapy. We investigated whether these critical genes were related to cisplatin resistance. The cells with knockdown of each specific gene were treated with cisplatin. Cell survival fraction and apoptosis were then examined. Our data showed that knockdown of BIRC5, TRIP13 or UBE2C significantly inhibited the cell survival fraction after cisplatin treatment (Figure 7A-C). Knockdown of these genes also significantly increased cisplatin-induced cell apoptosis (Figure 7D-F). These data suggested that high regulation of BIRC5, TRIP13 or UBE2C genes is closely related to cell proliferation as well as cisplatin resistance, and could be used as a potential target for treating GC.

Figure 7 Increased level of BIRC5, TRIP13, and UBE2C genes confers cellular resistance to cisplatin in gastric cancer cells. A-C: Colony formation assay in cells with shRNA transfected targeting BIRC5, TRIP13, or UBE2C genes; D-F: Cell apoptosis was detected in BIRC5, TRIP13, and UBE2C knockdown cells after treatment with cisplatin. ^aP < 0.05; ^bP < 0.01.

DISCUSSION

In the present research, in order to identify more core prognostic biomarkers in GC, we used bioinformatics methods based on three profile datasets (GSE13911, GSE66229, and GSE79973). Three hundred and forty-eight GC specimens and 141 normal specimens were included in this study. We discovered 251 changed DEGs (settings: |logFC| > 0.5 and adjusted P value < 0.05) including 64 up-regulated (logFC > 0) and 187 down-regulated DEGs (logFC < 0) using GEO2R and the Venn software platform. The GO and Pathway Enrichment Analysis using DAVID methods showed that: (1) For BP, up-regulated DEGs were particularly enriched in the regulation of nuclear chromosome segregation, sister chromatid segregation, cell division, cell cycle process, mitotic cell cycle process and down-regulated DEGs in organic acid catabolic process, carboxylic acid catabolic process, small molecule catabolic process, lipid modification, fatty acid catabolic process, and fatty acid oxidation; (2) for MF, up-regulated DEGs were enriched in identical protein binding, microtubule binding, enzyme binding, tubulin binding, cytoskeletal protein binding, protein kinase binding and down-regulated DEGs in cofactor binding, coenzyme binding, potassium channel regulator activity, ligand-gated ion channel activity, ligand-gated channel activity; and (3) for changes in CC, up-regulated DEGs were enriched in condensed chromosome kinetochore, condensed chromosome, centromeric region kinetochore, condensed chromosome, chromosome, centromeric region, and midbody, and down-regulated DEGs in mitochondrial matrix, mitochondrial part, mitochondrion, intrinsic component of plasma membrane, integral component of plasma membrane, and cytoplasmic vesicle membrane.

BUB1, which is from BUB (budding uninhibited by Benzimidazole) gene families, can encode proteins to form a large kinetochore-associated multi-protein complex. It has been suggested that altering BUB1 gene expression levels could disrupt its function, leading to aberrant chromosomal segregation. Stahl et al[18] indicated that tumors with a low frequency of BUB1 expression were associated with larger tumor size, higher incidence of lymph node metastases, distant metastases and higher UICC stage. Moreover, BUB1 expression was inversely correlated with the residual tumor stage. It was reported that the expression of BUB1 was also associated with taxane sensitivity in other types of cancer such as breast and esophageal cancer[19,20]. Furthermore, Kawakubo et al[21] found that BUB1 insufficiency leads to the lack of eNOS expression both in vitro and in vivo. They also revealed that there was a lack of eNOS expression in GC tissues with low BUBR1 expression.

Mitotic arrest deficient 2 (MAD2), is a key component of the mitotic spindle checkpoint pathway. It plays a vital role in maintaining spindle checkpoint function by generating special signals for monitoring the localization and segregation of chromosomes. Wang et al[22] also reported that the reduced expression of MAD2 resulted in a deficient mitotic checkpoint in several human cancers. Du et al[23] found that MAD2 interference could inhibit anticancer drug-induced apoptosis by up-regulating Bcl-2 and interfering with the mitochondria apoptosis pathway. Bargiela-Iparraguirre et al[24] demonstrated that depletion of MAD2 can delay Mcl-1 degradation, suggesting that exit to mitosis is delayed.

CCNB1 and CCNA2, which regulate cell progression via the G2/M transition during the cell cycle and encoded by related cyclin key genes, play significant roles in the G2/M phase. High levels of CCNA2[25] over-expression were detected in various cancers, including GC, non-small cell lung cancer and esophageal squamous cell carcinoma, resulting in uncontrolled cell growth. Up-regulation of CCNB1 expression was confirmed to promote the cell cycle in the G2/M phase, which was in line with poor survival in GC patients[26]. Moreover, Roncalli et al[27] indicated that CCNB1 arises in the cytoplasm of cells in S-phase, and is then transported to the nucleus at the G2/M transition where it is broken down during anaphase via a ubiquitin-dependent pathway. A recent study[28] showed that up-regulation of CCNA2 could be an indicator of pituitary adenoma invasiveness, as it played an essential role in monitoring protein in regulation of the cell cycle. Müssnich et al[29] found that downregulation of CCNA2 protein and mRNA levels decreased cell proliferation.

BIRC5, Baculoviral inhibitor of apoptosis repeat-containing 5, which is a well-known encoder of survivin, is a member of the inhibitor-of-apoptosis proteins (IAPs) family. During G2/M phase in tumor proliferating cells, BIRC5, located at 17q25, plays a key role in cell division, regulating cell apoptosis and the mitosis process. Several reports revealed that BIRC5 not only interacts with caspases but also enables the formation of the chromosomal passenger complex, which is vital for promoting mitosis[30-32]. Nabilsi et al[33] indicated that BIRC5 is related to sex hormones, and that overexpression of BIRC5 is associated with high levels of female sex hormones. Leung et al[34] reported that low expression of survivin may lead to a decrease in cell cycle function, which may result in the loss of related cells.

Based on patient survival analysis, we focused on three critical genes, including BIRC5, TRIP13, and UBE2C. It was found that these three genes were proved to play important roles in the progression of various cancers; however, our search of databases on these genes related to GC found only a few reports which suggested this correlation[35,36]. We also performed cell experiments and found that knockdown of these genes inhibited cell proliferation and reduced cellular resistance to cisplatin. Our findings provide novel targets for the treatment of GC.

CONCLUSION

In conclusion, exploring the molecular mechanism of GC via bioinformatics analysis showed that three DEGs (BIRC5, TRIP13, or UBE2C) play key roles in the progression of GC. Hub genes and key pathways of GC were identified using bioinformatics leading to a molecular link between the pathway and gene expression potentially involved in GC. The current study identified related genes and cellular pathways involved in the emergence and development of GC. However, certain limitations still existed due to the lack of cell experiments. Hence, further experimental studies should be carried out to confirm the expression and function of the identified genes at the molecular level, which will help to elucidate GC pathogenesis and identify novel biomarkers or drug targets for improved diagnostics and therapeutics for GC.

ARTICLE HIGHLIGHTS

Research background

Gastric cancer (GC) is one of the most common cancers and has a poor prognosis. Treatment of GC has remained unchanged over the past few years. Genetic factors have been revealed to play an important role in susceptibility to GC. Significantly upregulated genes associated with poor prognosis were detected in GC using integrated bioinformatics methods.

Research motivation

In order to identify core prognostic biomarkers in GC, several databases were searched for GC-related genes as tumor markers, and cellular tests were performed to confirm the results.

Research objectives

Bioinformatics analysis of the molecular mechanism involved in GC revealed that three differentially expressed genes (DEGs) (BIRC5, TRIP13, or UBE2C) play critical roles in the progression of GC. Bioinformatics were used to identify hub genes and important pathways in GC, resulting in a biological relationship between the pathways and gene expression likely involved in GC.

Research methods

In the theoretical analysis, microarray data information, data processing of DEGs, Gene Ontology and pathway enrichment analysis, PPI network and module analysis, and survival analysis were used. In the cellular experiments, RNA sequencing expression of core genes, patient samples and RT-PCR detection, cells and transfection, CCK-8 assay, apoptosis assay, colony formation assay, and statistical analysis were used.

Research results

Three hundred and forty-eight GC tissues and 141 normal tissues were analyzed; 251 DEGs were identified including 187 down-regulated genes and 64 up-regulated genes. We found that knockdown of BIRC5, TRIP13, or UBE2C significantly inhibited cell proliferation and induced cell apoptosis. Knockdown of BIRC5, TRIP13, or UBE2C increased cellular sensitivity to cisplatin.

Research conclusions

The molecular mechanism of GC, via bioinformatics analysis, showed that three DEGs (BIRC5, TRIP13, or UBE2C) played key roles in the progression of GC. These findings will help to elucidate GC pathogenesis and identify novel biomarkers or drug targets for improved diagnostics and therapeutics for GC.

Research perspectives

Bioinformatics was used to identify hub genes and important pathways in GC, resulting in a biological relationship between the pathways and gene expression likely involved in the progression of GC. Bioinformatics analysis revealed the relevant genes and cellular pathways involved in the genesis and development of GC.

ACKNOWLEDGEMENTS

We are grateful to Dr. Qingqiang Xu from Naval Medical University for kind discussions.