Skip to main content

ORIGINAL RESEARCH article

Front. Mol. Biosci., 15 October 2021
Sec. Molecular Diagnostics and Therapeutics
Volume 8 - 2021 | https://doi.org/10.3389/fmolb.2021.763652

Transcriptomic Analysis Identifies Complement Component 3 as a Potential Predictive Biomarker for Chemotherapy Resistance in Colorectal Cancer

www.frontiersin.orgXiao-Shun He1 www.frontiersin.orgSheng-Yi Zou2 www.frontiersin.orgJia-Lu Yao3 www.frontiersin.orgWangjianfei Yu4 www.frontiersin.orgZhi-Yong Deng5 www.frontiersin.orgJing-Ru Wang1 www.frontiersin.orgWen-Juan Gan6 www.frontiersin.orgShan Wan1* www.frontiersin.orgXiao-Qin Yang4* www.frontiersin.orgHua Wu1,6*
  • 1Department of Pathology, Medical College of Soochow University and The First Affiliated Hospital of Soochow University, Soochow University, Suzhou, China
  • 2Department of Endocrinology and Metabolism, The First Affiliated Hospital of Soochow University, Soochow University, Suzhou, China
  • 3Department of Cardiology, the First Affiliated Hospital of Soochow University, Soochow University, Suzhou, China
  • 4Department of Bioinformatics, Medical College of Soochow University, Soochow University, Suzhou, China
  • 5Department of Pathology, The First People’s Hospital of Kunshan, Kunshan, China
  • 6Department of Pathology, Dushu Lake Hospital Affiliated of Soochow University, Soochow University, Suzhou, China

Objective: 5-fluorouracil- and oxaliplatin-based FOLFOX regimens are mainstay chemotherapeutics for colorectal cancer (CRC) but drug resistance represents a major therapeutic challenge. To improve patient survival, there is a need to identify resistance genes to better understand the mechanisms underlying chemotherapy resistance.

Methods: Transcriptomic datasets were retrieved from the Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases and combined with our own microarray data. Weighted gene co-expression network analysis (WGCNA) was used to dissect the functional networks and hub genes associated with FOLFOX resistance and cancer recurrence. We then conducted analysis of prognosis, profiling of tumor infiltrating immune cells, and pathway overrepresentation analysis to comprehensively elucidate the biological impact of the identified hub gene in CRC.

Results: WGCNA analysis identified the complement component 3 (C3) gene as the only hub gene associated with both FOLFOX chemotherapy resistance and CRC recurrence after FOLFOX chemotherapy. Subsequent survival analysis confirmed that high C3 expression confers poor progression-free survival, disease-free survival, and recurrence-free survival. Further correlational analysis revealed significant negative association of C3 expression with sensitivity to oxaliplatin, but not 5-fluorouracil. Moreover, in silico analysis of tumor immune cell infiltration suggested the change of C3 expression could affect tumor microenvironment. Finally, gene set enrichment analysis (GSEA) revealed a hyperactivation of pathways contributing to invasion, metastasis, lymph node spread, and oxaliplatin resistance in CRC samples with C3 overexpression.

Conclusion: Our results suggest that high C3 expression is a debilitating factor for FOLFOX chemotherapy, especially for oxaliplatin sensitivity, and C3 may represent a novel biomarker for treatment decision of CRC.

Introduction

Colorectal cancer (CRC) is a critical global medical issue. While CRC cases and deaths have decreased over the last decades due to improved screening efforts and novel therapies, it has been estimated that CRC still ranked third in overall morbidity and mortality in the United States in 2021. The prognosis of CRC patient is strongly dependent on the clinical stage at diagnosis: the 5-years survival rate for stage I is 90%, compared with only 10% for stage IV (Virk et al., 2019). Surgery is the main treatment option during the early stages of colorectal cancer, but unfortunately, in many cases the disease is only detected in advanced stages or when distant metastasis has already occurred. Therefore, the discovery and establishment of new strategies for early diagnosis of colorectal cancer will be of great clinical significance.

CRC is a heterogeneous disease with a complex etiology involving an interplay of defective DNA damage repair, somatic mutations, chromosomal instability, microsatellite instability, DNA methylation, and epigenetic alterations (Colussi et al., 2013; Turano et al., 2019; Xu et al., 2021), all of which are thought to contribute to the malignant transformation of intestinal epithelial cells. These factors may also influence treatment response and prognosis of CRC patients. Transition from adenoma to adenocarcinoma is the classic model of CRC progression. Moreover, certain genetic alterations have been shown to be involved in the progression of CRC (Vogelstein et al., 2013). For instance, APC mutations are thought to occur first, followed by KRAS mutation and subsequent TP53 inactivation, resulting in abnormal activation of the Wnt/β-catenin, PI3K/Akt, NF-κB, JAK/STAT, and TGF-β/BMPs pathways (Prasad et al., 2010; Fearon, 2011; Vogelstein et al., 2013; Groner and von Manstein, 2017). Many genes such as APC, TP53, PIK3CA, SMAD4, FBXW7, BRAF, KRAS, and HER-2, amongst others, have been found to be commonly mutated in CRC (Muzny et al., 2012). Interestingly, FBXW7, HER-2 and TP53 have been found to be more frequently mutated during the early stages than in advanced stages (Testa et al., 2018; Yaeger et al., 2018). Moreover, TP53 and CTNNB1 mutations are thought to more frequently occur in younger patients while mutations in APC, KRAS, BRAF, and FAM123B are more likely to occur in older patients (Lieu et al., 2019). The liver is the most common organ for CRC metastasis, with at least 25% CRC patients exhibiting liver metastases following tumor progression (Martin et al., 2020). The main genetic mutations associated with CRC metastasis are mutations in KRAS and TP53, although some studies have found that SMAD4 and BRAF may also play an important role in CRC metastasis. We recently reported that TRAF6 can inhibit epithelial-mesenchymal transformation and CRC metastasis by driving selective β-catenin degradation via autophagy (Wu et al., 2019), while truncated RXRα (tRXRα) promotes colitis-associated colorectal tumorigenesis via activation of NF-κB-IL-6-STAT3 signaling in mouse (Ye et al., 2019). We recently also reported that the long non-coding RNA NONHSAT062994 acts as a tumor suppressor and inhibits CRC development via inactivation of Akt signaling (He et al., 2017).

The identification of circulating biomarkers, such as circulating tumor cells, exosomes, and miRNAs in blood (Turano et al., 2019), for early detection as well as monitoring treatment efficacy, recurrence, and metastasis of CRC has been a focus of research in recent years (Diaz et al., 2012; Tie et al., 2015). A limited number of targeted therapeutics, such as anti-EGFR drugs, have received approval for the treatment of colorectal cancer by the US FDA. However, as KRAS and NRAS act downstream of the EGFR pathway, CRC patients with KRAS mutations do not benefit from anti-EGFR therapy (Lieu et al., 2019). The mutation rate of BRAF in CRC is 8–13%, with the most common mutation site being V600E. BRAF mutation, especially in MSI-L (Microsatellite Instability-Low) and MSS (Microsatellite Stability), can enhance tumor resistance to anti-EGFR drugs, and result in a poor prognosis (Lech et al., 2016). Excitingly, new immunotherapies may elicit remission in more patients. The monoclonal antibodies Pembrolizumab and Nivolumab exert antitumor effects by blocking programmed cell death protein 1 (PD-1). They are suitable for patients with MSI-H (Microsatellite Instability-High) or dMMR (Dificient Mismatch Repair) and can also be used for patients with liver metastases, although patients with MSI-L or pMMR do not benefit from the anti-PD-1 treatment (Turano et al., 2019).

Despite the abovementioned efforts in cancer screening, treatment, and monitoring, the 5-years survival rate of CRC patients is still unsatisfactory. Therefore, the discovery of novel therapeutic strategies and identification of regulatory molecules is urgently needed. In this study, we combined transcriptomic data from public databases with our own work and found that complement component 3 (C3) gene was associated with oxaliplatin resistance via reprogramming the tumor immune microenvironment in CRC. These findings suggest C3 represents a potential predictive biomarker for the prognosis of oxaliplatin-based therapy in CRC patients.

Materials and Methods

Data Collection

RNA-Seq transcriptomic data of colon adenocarcinoma patients were retrieved from TCGA-COAD project via the R TCGAbiolinks package (version: 1.12.0) (Colaprico et al., 2015). The corresponding clinical information of the TCGA-COAD samples was obtained from the UCSC Xena database (Goldman et al., 2020). Microarray datasets (GSE39582, GSE81653, GSE87211, GSE28702 and GSE17536) for colorectal cancer samples and matching clinical characteristics were downloaded from the NCBI Gene Expression Omnibus (GEO) database. Microarray datasets for 5-Fluorouracil and oxaliplatin resistant cellular models (GSE76489 and GSE83131) were also obtained from GEO. Gene expression profiles of colorectal cancer cell lines and corresponding half-maximal inhibitory concentration (IC50) of 5-Fluorouracil and oxaliplatin were obtained from the Genomics of Drug Sensitivity in Cancer (GDSC) database (Yang et al., 2013).

Clinical Specimens

Colorectal cancer (CRC) samples (four with lymph node metastasis and three without lymph node metastasis) were obtained from patients who underwent curative surgery at the First Affiliated Hospital of Soochow University (Suzhou, Jiangsu, P.R. China). All patients included in this study provided written informed consent for participation. The study was approved by the Biomedical Research Ethics Committee of Soochow University (Suzhou 215123, Jiangsu Province, China).

Microarray Experiment

Total RNA of each sample was extracted using TRIzol reagent. The Agilent Human lncRNA-mRNA Microarray V2.0 4 × 180 K (Agilent Technologies, Inc.) was used to compare transcriptomic profiles between subgroups with and without lymph node metastasis. RNA labeling, microarray hybridization, and data acquisition were performed by Shanghai OE Biotech. Co., Ltd. The raw microarray data reported in this paper has been deposited in the OMIX, China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences, under the accession number OMIX505 and can be accessed at https://ngdc.cncb.ac.cn/omix.

High Throughput Data Preprocessing and Assessment of Differentially Expressed Genes

For TCGA-COAD RNA-Seq data, the raw reads count matrix was converted to the counts per million mapped reads (CPM) format. Trimmed mean of M values (TMM) normalization was carried out to correct the variation of gene expression abundance across different samples. Both log2-transformed or raw formats were generated for further bioinformatics analyses. Samples were categorized according to the median expression level of target genes. Differentially expressed genes were filtered using the Benjamini-Hochberg false discovery rate approach (adjusted p-value < 0.05, and absolute value of log2 fold-change > 1). All preprocessing and differential gene expression analysis was conducted using an R pipeline based on “limma” (version: 3.48.1) (Ritchie et al., 2015) and “edgeR” (version: 3.34.0) (Robinson et al., 2010) R packages.

For Affymetrix microarray data, data was normalized using the Robust Multi-array Average (RMA) method. For Agilent microarray data, data were normalized using the “quantile” method. Gene level expression estimates were log2-transformed and summarized into a single matrix for subsequent analysis. Raw p-values < 0.01 were used as the cutoff to identify differentially expressed genes. All microarray analyses were performed using the “limma” R package (version: 3.48.1).

Weighted Gene Co-expression Network Analysis

Genes were ranked according to the median absolute division (MAD) values across samples in the GSE81653 and GSE28702 dataset, and the top 5,000 genes were subjected to the subsequent weighted gene co-expression network analysis (WGCNA). Briefly, a minimum power with R2 > 0.85 was set to optimize the power (β) for automatic network and further topological overlap matrix construction. After that, minModuleSize = 30 was utilized to initially classify co-expression modules. In one data set with less than 100 samples (GSE28702), the cutHeight was set to 0.2. Conversely, when the sample size exceeded 100 (GSE81653), this parameter was set to 0. A p-value less than 0.05 indicated statistical significance for the association between co-expression modules and clinical traits. Module membership (MM) values >0.7 served as the cutoff for key gene screening. The intersection between the results of GSE81653 and GSE28702 was taken for further analysis. All computational analyses were performed using the “WGCNA” R package (version: 1.70-3) (Langfelder and Horvath, 2008).

Interaction Network Analysis

The STRING database (version 11.0) (Szklarczyk et al., 2019) was used to identify a protein-protein interaction network linked by the target gene connecting the significant gene modules yielded by WGCNA and the up-regulated genes in the colorectal cancer samples with lymph node metastasis (filtered by p-value < 0.05 in the “limma” workflow). The “clusterProfiler” R package (version: 4.0.2) (Wu et al., 2021) was used to annotate the biological function of the genes in the network. A raw p-value less than 0.05 indicated statistically significant enrichment.

Gene Set Enrichment Analysis

Gene set enrichment analysis (GSEA) was performed to determine whether certain biological pathways or a priori defined gene sets showed statistically significant differences between two subgroups with high and low expression of the target gene. The “clusterProfiler” R package (version: 4.0.2) was used for enrichment analyses (Wu et al., 2021). Gene sets including Chemical and Genetic Perturbations, KEGG pathway, and Gene Ontology Biological Process, were accessed using the Molecular Signatures Database (MSigDB). Gene markers for lymph node metastasis were defined as the top 300 significantly up-regulated genes (p-value less than 0.01) ranked by log2-transformed fold change obtained from our own microarray data. A raw p-value less than 0.05 indicated statistically significant enrichment.

GSEA was also used to determine whether the gene signatures that changed upon the expression of the target gene showed statistically significant, concordant differences in 5-Fluorouracil or Oxaliplatin resistant and naïve cell lines. Significant differentially expressed genes between the TCGA-COAD subgroups with high and low expression of the target gene were ranked according to the log2-transformed fold change. The top 300 up- and down-regulated genes were defined as the gene sets. GSEA was then applied to examine the enrichment of these gene sets in the GSE76489 (5-Fluorouracil resistance model) and GSE83131 (Oxaliplatin resistance model) datasets, respectively. A raw p-value less than 0.05 indicated statistically significant enrichment.

In Silico Assessment of Chemotherapy Reagent Sensitivity

The Pearson correlation between log2-transformed levels of the target gene expression and ln-transformed IC50 from the GDSC database were calculated. A p-value less than 0.05 indicated statistical significance.

Survival Analysis

We used the GSE39582 dataset for evaluation of recurrence risk. Samples were categorized according to the median expression level of the target gene. Odds ratios (ORs) and corresponding 95% confidence intervals (95% CIs) were calculated to assess the risk of colorectal cancer recurrence. Stratified analyses according to TNM stage, size of primary tumor (T), regional lymph node involvement (N), presence of distant metastatic spread (M), and mutation status of TP53, KRAS, and BRAF were performed to assess whether the effect of the target gene is influenced by the above-mentioned confounding factors.

TCGA-COAD, GSE17536, GSE87211 and GSE39582 datasets were used for prognostic analyses to assess the impact of the target gene on colorectal cancer survival. The optimal cutoff values for each dataset were determined by the algorithms embedded in the “survminer” (version: 0.4.9) R package. The logrank test was performed to compare the survival distributions between the assigned subgroups. A p-value less than 0.05 indicated statistical significance.

Analysis of Tumor Infiltrating Immune Cells

CPM reads of the TCGA-COAD dataset were assessed for immune cell infiltration analysis. Samples were categorized according to the median expression level of the target gene. The cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT) method (Newman et al., 2015) was used to profile tumor infiltrating immune cells and compare the differences between the two subgroups. A p-value less than 0.05 indicated statistical significance.

Statistical Analysis

All bioinformatics analyses were carried out in the R programming environment (version: 4.1.0).

Results

WGCNA Analyses Identifies C3 as the Common Gene Associated With FOLFOX Chemotherapy Resistance and Cancer Recurrence Risk

To identify key genes associated with specific clinical traits of colorectal cancer, we employed the WGCNA method to identify gene module-trait correlations. We categorized the top 5,000 genes with the highest MAD in the GSE28702 dataset into 22 co-expression modules (Figure 1A). The MEgrey60 and Melightgreen modules, containing 195 genes, were positively associated with FOLFOX chemotherapy resistance (Figure 1B) and we selected 44 key genes with a MM value greater than 0.7. Next, we performed WGCNA analysis for GSE81653 using the same methods and parameters as previously to explore the relevant factors of recurrence risk after FOLFOX therapy which resulted in clustering of 7 co-expression modules (Figure 1C). The Mebrown and Meblue modules consisting of a total of 216 genes were positively associated with recurrence risk (Figure 1D) and we again selected 114 hub genes with a MM value higher than 0.7. We then intersected the key gene members from both WGCNA analyses corresponding to FOLFOX resistance and tumor recurrence which identified the complement component 3 (C3) gene as the only common element (Figure 1E).

FIGURE 1
www.frontiersin.org

FIGURE 1. Weighted gene co-expression network analysis (WGCNA) predicts that C3 gene aggravates FOLFOX chemotherapy resistance and recurrence of colorectal cancer. (A) Clustering dendrogram of the top 5,000 most variant genes in the GSE28702 dataset, with dissimilarities based on topological overlap labelled with assigned module colors (B) Heatmap for the correlation of WGCNA modules with clinical traits (GSE28702). (C) Clustering dendrogram of top 5,000 most variant genes in the GSE81653 dataset, with dissimilarities based on topological overlap labelled with assigned module colors (D) Heatmap for the correlation of WGCNA modules with clinical traits (GSE81653). (E) Venn plot indicated that the C3 gene was the only common element between the chemotherapy respondence-associated module (WGCNA of GSE28702) and the colorectal cancer recurrence-associated module (WGCNA of GSE81653). Genes with module membership values >0.7 were used for this intersection analysis. A p-value less than 0.05 indicated statistical significance.

C3 Links the Subnetworks Contributing to Poor Prognosis of Colorectal Cancer

We then constructed a protein-protein interaction network by integrating the interaction relationships of the significant modules from the WGCNA analyses above and the differentially expressed genes between samples with and without lymph node metastasis to delineate the core genes connecting subnetworks responsible for FOLFOX chemotherapy resistance, tumor recurrence, and spread to lymph nodes. As highlighted in the interaction network (Figure 2A), the C3 gene represents one of the bridge hubs linking the three subnetworks. Further functional annotation analysis identified a significant enrichment of gene sets related to cancer stemness, invasiveness, migration, proliferation, angiogenesis, and inflammation response. In addition, cellular responses to chemical stimulus, a pathway associated with chemotherapy resistance, was also significantly enriched (Figure 2B). Based on these results, we propose that the C3 gene may act as a potential regulator of the prognostic CRC network and serves as a shared hub gene, rather than as a node with scattered functional edges participating in a single signaling pathway.

FIGURE 2
www.frontiersin.org

FIGURE 2. C3 interacts with genes associated with poor prognosis of colorectal cancer. (A) The interaction network of C3. Genes from significant modules related FOLFOX chemotherapy resistance are shown in the blue circle, genes from the significant module related to colorectal cancer recurrence are shown in the in red circle, and genes upregulated in colorectal cancer samples with lymph node metastasis are shown in the green circle (B) Enrichment analysis for genes in the C3 gene interaction network. A raw p-value less than 0.05 indicated statistical significance.

High C3 Gene Expression Leads to Increased Colorectal Cancer Recurrence Risk and Poor Survival

We next compared the effect of high and low expression of C3 on tumor recurrence in patients with CRC in the GSE39582 dataset. Patients with high C3 expression exhibited an increased risk of recurrence (OR = 1.52, 95% CI: 1.05–2.22; p = 0.023) (Figure 3). Further subgroup analysis showed that the levels of C3 expression had significantly positive effect on colorectal cancer recurrence in patients during early TNM stages (p = 0.002) and those without spread to lymph nodes (p = 0.009). In addition, we found a statistically significant susceptibility in subgroups without mutations in TP53 (p = 0.010), KRAS (p = 0.004), and BRAF (p = 0.019).

FIGURE 3
www.frontiersin.org

FIGURE 3. Forest plot of the association between C3 gene expression and the recurrence risk of colorectal cancer. Samples in the GSE39582 dataset were categorized according to the median expression C3 gene. Odds ratio (OR) and 95% confidence interval (95% CI) were calculated to assess the risk of recurrence. M: mutation; WT: wild type. *p < 0.05, **p < 0.01.

Subsequent survival analysis was performed to validate the prognostic value of C3 in four independent transcriptome datasets with a total of 1,358 samples. Kaplan-Meier plots revealed that high C3 expression was associated with a statistically significant unfavorable progression−free survival (Figure 4A, p = 0.023), disease−free survival (Figures 4B,C; both p-values < 0.05), and recurrence−free survival (Figure 4D, p = 0.022). Therefore, we propose that C3 represents a candidate biomarker for poor prognosis in colorectal cancer patients.

FIGURE 4
www.frontiersin.org

FIGURE 4. C3 gene expression is a significant prognostic factor of colorectal cancer. Kaplan-Meier curves assessing progression-free survival in the TCGA-COAD dataset (A), disease-free survival in the GSE17536 dataset (B) and GSE87211 dataset (C), and recurrence-free survival in the GSE39582 dataset (D) were plotted. Optimal separation based on C3 gene expression for samples in each dataset was identified to achieve best statistical significance. A p-value less than 0.05 indicated statistical significance.

C3 Overexpression Promotes Oxaliplatin Resistance

To explore the impact of C3 gene expression on resistance to FOLFOX chemotherapy, we assessed the correlation between C3 gene expression in different colorectal cancer cell lines and the corresponding IC50 of 5-Fluorouracil and Oxaliplatin. C3 gene expression was positively associated with Oxaliplatin IC50 (Figure 5 A, p = 0.019), but not with 5-Fluorouracil (Figure 5B, p = 0.120). Next, we compared the transcriptome concordance between chemotherapy reagent resistance and C3 overexpression. The top 300 up-regulated genes were obtained by comparing the samples with a high and low abundance of C3 expression (no significantly down-regulated gene were identified because of insufficient log2-transformed fold changes). The GSEA results found these genes were significantly enriched in the comparison between Oxaliplatin-resistant cells and naïve cells (Figure 5C, p = 0.003). However, the enrichment was not statistically significant in the cell model of 5-Fluorouracil resistance (Figure 5D, p = 0.832). This was consistent with the results of in silico IC50 tests. Based on the observations above, it could be reasonably speculated that a high C3 gene expression in colorectal cancer might promote resistance to Oxaliplatin.

FIGURE 5
www.frontiersin.org

FIGURE 5. High expression of C3 induces resistance to oxaliplatin but not 5-Fluorouracil. Pearson correlation test revealed a statistically significant association of C3 gene expression with the corresponding IC50 of oxaliplatin (A), but not 5-Fluorouracil (B), in colorectal cancer cell lines. GSEA demonstrated that the top 300 upregulated genes in colorectal cancers with high C3 expression were significantly enriched in Oxaliplatin resistant cells (C), but not in 5-Fluorouracil resistant cells (D). A raw p-value less than 0.05 indicated statistical significance.

High C3 Gene Expression Shapes the Immune Microenvironment of Colorectal Cancer

To explore the impact of C3 gene expression on the tumor immune microenvironment, we employed the CIBERSORT algorithm to compare the abundance of 22 human infiltrating immune cell types between colorectal cancers with C3 gene expression above and below median level (Figure 6). This revealed a large heterogeneity of infiltrating immune cells between the two different C3 gene expression groups. There were significantly more infiltrating dendritic cells activated, mast cells activated, monocytes, and T cells CD4 memory resting in colorectal cancers with C3 lower expression compared to those with higher C3 expression (all p-values < 0.001). Conversely, tumors with high C3 expression showed significantly higher fractions of M0 macrophages, M1 macrophages, and mast cells resting (all p-values < 0.01). No significant difference of other immune cells fraction was found (all p-values > 0.05).

FIGURE 6
www.frontiersin.org

FIGURE 6. Comparison of the abundance of infiltrating immune cells between high-C3 and low-C3 groups in the TCGA-COAD dataset. Samples were categorized according to the median expression level of C3 gene. Kruskal-Wallis test was performed to statistically compare differences. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.

C3 Overexpression Drives Key Signaling Cascades Regulating Colorectal Cancer Progression

We then performed GSEA to elucidate the biological pathways regulated by C3 that are involved in the progression of CRC. This revealed that clinical CRC samples with higher C3 expression were enriched for the general pathway of colorectal cancer development and progression (Figure 7A). Furthermore, lymphatic metastasis signatures identified by our microarray data were also significantly overrepresented (Figure 7B) and likewise, genes associated with invasion of lymphatic vessels during metastasis also showed marked enrichment (Figure 7C). In addition, further gene set mining suggested that higher C3 expression could promote cancer proliferation (Figure 7D), invasion (Figure 7E), and adhesion (Figure 7F). Taken together, these observations corroborated the hypothesis that a higher C3 gene expression could drive poor prognosis of colorectal cancer. To further uncover the underlying mechanisms, we next investigated alterations of biological pathways upon C3 hyperexpression. Colorectal cancers with high C3 gene level exhibited a hyperactivation of several pathways including pathways in cancer (Figure 7G), TGF-BETA signaling (Figure 7H), JAK-STAT signaling (Figure 7I), MAPK signaling (Figure 7J), VEGF signaling (Figure 7K), and WNT signaling (Figure 7L).

FIGURE 7
www.frontiersin.org

FIGURE 7. Identification of key pathways and gene sets associated with colorectal cancer prognosis based on TCGA-COAD RNA-Seq dataset via gene set enrichment analysis (GSEA). Gene sets for colorectal cancer (A), lymphatic metastasis (B), lymphatic vessels during metastasis (C), proliferation (D), multicancer invasiveness signature (E), focal adhesion (F), pathways in cancer (G), TGF-BETA signaling (H), JAK-STAT signaling (I), MAPK signaling (J), VEGF signaling (K), and WNT signaling (L) were significantly enriched. Samples were categorized according to the median expression level of C3 gene. A raw p-value less than 0.05 indicated statistical significance.

To summarize, overexpression C3 gene is likely to activate several signaling pathways and plays a prominent role in initiating cancerous deterioration, driving poor prognosis.

Discussion

Surgery and radiation are curative options for patients during early stages of CRC. Conversely, CRC patients at advanced stages typically receive chemotherapy to reduce the rate of recurrence and improve the survival rate. Oxaliplatin, a third-generation platinum-based anticancer agents, is approved in combination with fluorouracil (5-FU) or leucovorin (LV) for advanced CRC treatment (Mandala et al., 2004; Yothers, 2012). Although certain CRC patients typically benefit from treatment with oxaliplatin, many acquire resistance and thus their prognosis can be poor (Culy et al., 2000). It is therefore of great importance to identify potential biomarkers for oxaliplatin resistance in order to facilitate treatment decisions. In this study, we identified complement component 3 (C3) as a critical gene associated with FOLFOX chemotherapy resistance by integration and analysis of public transcriptomic data and our own microarray data.

C3 is a key molecule of the enzymatic complement cascade reaction which acts as a powerful pro-inflammatory factor. It can activate a variety of signaling pathways, such as NF-κB, IRF, AP-1 pathway, resulting in inflammation (Tam et al., 2014). Studies have shown that C3 acts as a potential biomarker for inflammation and immunity-related disease. For examples, levels of plasma C3 are significantly elevated in patients with rheumatoid arthritis and axial spondyloarthritis (Arias de la Rosa et al., 2020). Likewise, C3 levels were significantly upregulated in the mucosae of patients with inflammatory bowel disease (Sugihara et al., 2010). Inflammation and tumor often go hand in hand, and chronic inflammation can alter the microenvironment in favor of tumor formation and growth (Balkwill and Mantovani, 2001). Moreover, the complement can increase the production of tumor growth factors, inhibit apoptosis, promote angiogenesis, and inhibit anti-tumor immunity to promote oncogenesis (Rutkowski et al., 2010). C3 has previously been shown to play critical roles in several cancer types, including glioblastoma, ovarian cancer, skin squamous cell carcinoma, and CRC. For instance, deposition of C3 has been found to be higher in glioblastoma than in controls (Bouwens et al., 2015). Ovarian cancer cells can secret C3, which promotes tumor growth and metastasis (Cho et al., 2014; Cho et al., 2016). C3 is upregulated in skin squamous cell carcinoma and promotes the growth of cutaneous squamous cell carcinoma (Riihila et al., 2017). Plasma C3 levels have previously been shown to be significantly increased in CRC (Dowling et al., 2012). In a mouse model of colorectal tumor transplantation, complement inhibition and complement depletion were used to reduce the plasma levels C3 and it was found that a reduction of C3 levels inhibited tumor growth, while tumor progression was resumed upon replenishment of C3 (Downs-Canner et al., 2016). The abovementioned studies have provided important evidence about the oncogenic function of C3. However, the potential role and clinical application of C3 in CRC had remained largely undefined. Here, we demonstrate that C3 is a key gene associated with FOLFOX chemotherapy resistance and cancer recurrence risk. The prognostic power of C3 for CRC is supported by the clinical data showing that high C3 expression correlates with an unfavourable prognosis as evidenced by poor progression-free survival and disease-free survival as well as an increased CRC recurrence risk. Thus, C3 expression levels could be used as an index for predicting the outcome of CRC chemotherapy in future clinical studies.

Mechanisms of oxaliplatin resistance are complex and involve many molecular events including reduced drug uptake, decreased apoptotic response, and enhanced tumor immunosuppressive effect (Martinez-Balibrea et al., 2015; Ning et al., 2021). There is conflicting evidence regarding the mechanism of oxaliplatin resistance in cancer. On one hand, oxaliplatin treatment results in an increased proportion of CD4+/CD8+ T lymphocytes and decreased proportion of regulatory T cells (Tregs) (Stojanovska et al., 2019). Conversely, an increased proportion of Tregs has also previously been observed in peripheral blood of CRC patients that received combined oxaliplatin and 5-FU treatment (Vizio et al., 2012). In this study, we found a large heterogeneity of infiltrating immune cells. Patients with different C3 expression exhibited significantly differential infiltration of activated dendritic cells and mast cells, mast cells resting, M0 and M1 macrophages, monocytes, and resting CD4+ memory T cell. GSEA results further showed that gene signatures representing oncogenic signaling and metastasis, such as pathways involved in cancer, TGFβ signaling, JAK-STAT signaling, MAPK signaling, VEGF signaling, and WNT signaling, were all significantly enriched in the patients with higher C3 expression compared to those with lower C3 expression. Taken together, these data indicate that overexpression of C3 confers CRC cells resistance to oxaliplatin treatment, The resistance mechanism may likely involve a reprogramming of the tumor immune microenvironment and oncogenic signaling.

In summary, we identify C3 as an oxaliplatin resistance-related gene using both public datasets and our own microarray data. These findings highlight a potential clinical role of C3 as a predictive biomarker for oxaliplatin-based therapy response and provide a theoretical guidance and decision-making for oxaliplatin-based chemotherapy in CRC patients.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://ngdc.cncb.ac.cn/omix/release/OMIX505

Ethics Statement

The studies involving human participants were reviewed and approved by Soochow University for Biomedical Research Ethics Committee. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

XH collected the human CRC samples and clinical data. SZ, JY, WY and XY analyzed the transcriptomic data. ZD and WG provided technical support. XH, SW, XY and HW designed the experiments, analyzed the data and wrote the manuscript. All authors read and approved the final manuscript.

Funding

This work was supported by Natural Science Foundation of Jiangsu Province (BK20190042, BK20181434 and BK20190182), and National Natural Science Foundation of China (82022050, 81902400, 81972601 and 81772541). This work was also supported by the Fujian Provincial Key Laboratory of Innovative Drug Target Research, the Tang Scholar Funds, and the Priority Academic Program Development of Jiangsu Higher Education Institutions.

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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2021.763652/full#supplementary-material

References

Arias de la Rosa, I., Font, P., Escudero-Contreras, A., López-Montilla, M. D., Pérez-Sánchez, C., Ábalos-Aguilera, M. C., et al. (2020). Complement Component 3 as Biomarker of Disease Activity and Cardiometabolic Risk Factor in Rheumatoid Arthritis and Spondyloarthritis. Ther. Adv. Chronic Dis. 11, 2040622320965067. doi:10.1177/2040622320965067

PubMed Abstract | CrossRef Full Text | Google Scholar

Balkwill, F., and Mantovani, A. (2001). Inflammation and Cancer: Back to Virchow? Lancet 357, 539–545. doi:10.1016/s0140-6736(00)04046-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Bouwens, T. A. M., Trouw, L. A., Veerhuis, R., Dirven, C. M. F., Lamfers, M. L. M., and Al-Khawaja, H. (2015). Complement Activation in Glioblastoma Multiforme Pathophysiology: Evidence from Serum Levels and Presence of Complement Activation Products in Tumor Tissue. J. Neuroimmunol. 278, 271–276. doi:10.1016/j.jneuroim.2014.11.016

CrossRef Full Text | Google Scholar

Cho, M. S., Vasquez, H. G., Rupaimoole, R., Pradeep, S., Wu, S., Zand, B., et al. (2014). Autocrine Effects of Tumor-Derived Complement. Cel Rep. 6, 1085–1095. doi:10.1016/j.celrep.2014.02.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Cho, M. S., Rupaimoole, R., Choi, H.-J., Noh, K., Chen, J., Hu, Q., et al. (2016). Complement Component 3 Is Regulated by TWIST1 and Mediates Epithelial-Mesenchymal Transition. J.I. 196, 1412–1418. doi:10.4049/jimmunol.1501886

PubMed Abstract | CrossRef Full Text | Google Scholar

Colaprico, A., Silva, T. C., Olsen, C., Garofano, L., Cava, C., Garolini, D., et al. (2015). TCGAbiolinks: an R/Bioconductor Package for Integrative Analysis of TCGA Data. Nucleic Acids Res. 44, e71. doi:10.1093/nar/gkv1507

PubMed Abstract | CrossRef Full Text | Google Scholar

Colussi, D., Brandi, G., Bazzoli, F., and Ricciardiello, L. (2013). Molecular Pathways Involved in Colorectal Cancer: Implications for Disease Behavior and Prevention. Ijms 14, 16365–85. doi:10.3390/ijms140816365

PubMed Abstract | CrossRef Full Text | Google Scholar

Culy, C. R., Clemett, D., and Wiseman, L. R. (2000). Oxaliplatin. A Review of its Pharmacological Properties and Clinical Efficacy in Metastatic Colorectal Cancer and its Potential in other Malignancies. Drugs 60, 895–924. doi:10.2165/00003495-200060040-00005

PubMed Abstract | CrossRef Full Text | Google Scholar

Diaz Jr, L. A., Williams, R. T., Wu, J., Kinde, I., Hecht, J. R., Berlin, J., et al. (2012). The Molecular Evolution of Acquired Resistance to Targeted EGFR Blockade in Colorectal Cancers. Nature 486, 537–540. doi:10.1038/nature11219

PubMed Abstract | CrossRef Full Text | Google Scholar

Dowling, P., Clarke, C., Hennessy, K., Torralbo-Lopez, B., Ballot, J., Crown, J., et al. (2012). Analysis of Acute-phase Proteins, AHSG, C3, CLI, HP and SAA, Reveals Distinctive Expression Patterns Associated with Breast, Colorectal and Lung Cancer. Int. J. Cancer 131, 911–923. doi:10.1002/ijc.26462

PubMed Abstract | CrossRef Full Text | Google Scholar

Downs-Canner, S., Magge, D., Ravindranathan, R., O’Malley, M. E., Francis, L., Liu, Z., et al. (2016). Complement Inhibition: A Novel Form of Immunotherapy for Colon Cancer. Ann. Surg. Oncol. 23, 655–662. doi:10.1245/s10434-015-4778-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Fearon, E. R. (2011). Molecular Genetics of Colorectal Cancer. Annu. Rev. Pathol. Mech. Dis. 6, 479–507. doi:10.1146/annurev-pathol-011110-130235

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldman, M. J., Craft, B., Hastie, M., Repečka, K., Mcdade, F., Kamath, A., et al. (2020). Visualizing and Interpreting Cancer Genomics Data via the Xena Platform. Nat. Biotechnol. 38, 675–678. doi:10.1038/s41587-020-0546-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Groner, B., and von Manstein, V. (2017). Jak Stat Signaling and Cancer: Opportunities, Benefits and Side Effects of Targeted Inhibition. Mol. Cell Endocrinol. 451, 1–14. doi:10.1016/j.mce.2017.05.033

PubMed Abstract | CrossRef Full Text | Google Scholar

He, X.-S., Guo, L.-C., Du, M.-Z., Huang, S., Huang, R.-P., Zhan, S.-H., et al. (2017). The Long Non-coding RNA NONHSAT062994 Inhibits Colorectal Cancer by Inactivating Akt Signaling. Oncotarget 8, 68696–68706. doi:10.18632/oncotarget.19827

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: an R Package for Weighted Correlation Network Analysis. BMC Bioinf. 9, 559. doi:10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Lech, G., Slotwinski, R., Slodkowski, M., and Krasnodebski, I. W. (2016). Colorectal Cancer Tumour Markers and Biomarkers: Recent Therapeutic Advances. World J. Gastroenterol. 22, 1745–1755. doi:10.3748/wjg.v22.i5.1745

CrossRef Full Text | Google Scholar

Lieu, C. H., Golemis, E. A., Serebriiskii, I. G., Newberg, J., Hemmerich, A., Connelly, C., et al. (2019). Comprehensive Genomic Landscapes in Early and Later Onset Colorectal Cancer. Clin. Cancer Res. 25, 5852–5858. doi:10.1158/1078-0432.ccr-19-0899

PubMed Abstract | CrossRef Full Text | Google Scholar

Mandala, M., Ferretti, G., and Barni, S. (2004). Oxaliplatin in colon Cancer. N. Engl. J. Med. 351, 1691–1692. doi:10.1056/nejm200410143511623

CrossRef Full Text | Google Scholar

Martin, J., Petrillo, A., Smyth, E. C., Shaida, N., Khwaja, S., Cheow, H., et al. (2020). Colorectal Liver Metastases: Current Management and Future Perspectives. World J. Clin. Oncol. 11, 761–808. doi:10.5306/wjco.v11.i10.761

CrossRef Full Text | Google Scholar

Martinez-Balibrea, E., Martínez-Cardús, A., Ginés, A., Ruiz De Porras, V., Moutinho, C., Layos, L., et al. (2015). Tumor-Related Molecular Mechanisms of Oxaliplatin Resistance. Mol. Cancer Ther. 14, 1767–1776. doi:10.1158/1535-7163.mct-14-0636

PubMed Abstract | CrossRef Full Text | Google Scholar

Muzny, D. M., Bainbridge, M. N., Chang, K., Dinh, H. H., Drummond, J. A., Fowler, G., et al. (2012). Comprehensive Molecular Characterization of Human colon and Rectal Cancer. Nature 487, 330–337. doi:10.1038/nature11252

PubMed Abstract | CrossRef Full Text | Google Scholar

Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust Enumeration of Cell Subsets from Tissue Expression Profiles. Nat. Methods 12, 453–457. doi:10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

Ning, T., Li, J., He, Y., Zhang, H., Wang, X., Deng, T., et al. (2021). Exosomal miR-208b Related with Oxaliplatin Resistance Promotes Treg Expansion in Colorectal Cancer. Mol. Ther. 29 (9), 2723–2736. doi:10.1016/j.ymthe.2021.04.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Prasad, S., Ravindran, J., and Aggarwal, B. B. (2010). NF-κB and Cancer: How Intimate Is This Relationship. Mol. Cel Biochem. 336, 25–37. doi:10.1007/s11010-009-0267-2

CrossRef Full Text | Google Scholar

Riihilä, P., Nissinen, L., Farshchian, M., Kallajoki, M., Kivisaari, A., Meri, S., et al. (2017). 118 Complement Components C3 and Complement Factor B Promote Growth of Cutaneous Squamous Cell Carcinoma. J. Invest. Dermatol. 137, S20. doi:10.1016/j.jid.2017.02.132

CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43, e47. doi:10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, M. D., Mccarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics 26, 139–140. doi:10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

Rutkowski, M. J., Sughrue, M. E., Kane, A. J., Mills, S. A., and Parsa, A. T. (2010). Cancer and the Complement Cascade. Mol. Cancer Res. 8, 1453–1465. doi:10.1158/1541-7786.mcr-10-0225

PubMed Abstract | CrossRef Full Text | Google Scholar

Stojanovska, V., Prakash, M., Mcquade, R., Fraser, S., Apostolopoulos, V., Sakkal, S., et al. (2019). Oxaliplatin Treatment Alters Systemic Immune Responses. Biomed. Res. Int. 2019, 4650695. doi:10.1155/2019/4650695

PubMed Abstract | CrossRef Full Text | Google Scholar

Sugihara, T., Kobori, A., Imaeda, H., Tsujikawa, T., Amagase, K., Takeuchi, K., et al. (2010). The Increased Mucosal mRNA Expressions of Complement C3 and Interleukin-17 in Inflammatory Bowel Disease. Clin. Exp. Immunol. 160, 386–393. doi:10.1111/j.1365-2249.2010.04093.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Gable, A. L., Lyon, D., Junge, A., Wyder, S., Huerta-Cepas, J., et al. (2019). STRING V11: Protein-Protein Association Networks with Increased Coverage, Supporting Functional Discovery in Genome-wide Experimental Datasets. Nucleic Acids Res. 47, D607–d613. doi:10.1093/nar/gky1131

PubMed Abstract | CrossRef Full Text | Google Scholar

Tam, J. C. H., Bidgood, S. R., Mcewan, W. A., and James, L. C. (2014). Intracellular Sensing of Complement C3 Activates Cell Autonomous Immunity. Science 345, 1256070. doi:10.1126/science.1256070

PubMed Abstract | CrossRef Full Text | Google Scholar

Testa, U., Pelosi, E., and Castelli, G. (2018). Colorectal Cancer: Genetic Abnormalities, Tumor Progression, Tumor Heterogeneity, Clonal Evolution and Tumor-Initiating Cells. Med. Sci. (Basel) 6, 31. doi:10.3390/medsci6020031

CrossRef Full Text | Google Scholar

Tie, J., Kinde, I., Wang, Y., Wong, H. L., Roebert, J., Christie, M., et al. (2015). Circulating Tumor DNA as an Early Marker of Therapeutic Response in Patients with Metastatic Colorectal Cancer. Ann. Oncol. 26, 1715–1722. doi:10.1093/annonc/mdv177

PubMed Abstract | CrossRef Full Text | Google Scholar

Turano, M., Delrio, P., Rega, D., Cammarota, F., Polverino, A., Duraturo, F., et al. (2019). Promising Colorectal Cancer Biomarkers for Precision Prevention and Therapy. Cancers (Basel) 11, 1932. doi:10.3390/cancers11121932

PubMed Abstract | CrossRef Full Text | Google Scholar

Virk, G. S., Jafri, M., Mehdi, S., and Ashley, C. (2019). Staging and Survival of Colorectal Cancer (CRC) in Octogenarians: Nationwide Study of US Veterans. J. Gastrointest. Oncol. 10, 12–18. doi:10.21037/jgo.2018.09.01

PubMed Abstract | CrossRef Full Text | Google Scholar

Vizio, B., Novarino, A., Giacobino, A., Cristiano, C., Prati, A., Ciuffreda, L., et al. (2012). Potential Plasticity of T Regulatory Cells in Pancreatic Carcinoma in Relation to Disease Progression and Outcome. Exp. Ther. Med. 4, 70–78. doi:10.3892/etm.2012.553

PubMed Abstract | CrossRef Full Text | Google Scholar

Vogelstein, B., Papadopoulos, N., Velculescu, V. E., Zhou, S., Diaz, L. A., and Kinzler, K. W. (2013). Cancer Genome Landscapes. Science 339, 1546–1558. doi:10.1126/science.1235122

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, H., Lu, X.-X., Wang, J.-R., Yang, T.-Y., Li, X.-M., He, X.-S., et al. (2019). TRAF6 Inhibits Colorectal Cancer Metastasis through Regulating Selective Autophagic CTNNB1/β-Catenin Degradation and Is Targeted for GSK3B/GSK3β-Mediated Phosphorylation and Degradation. Autophagy 15, 1506–1522. doi:10.1080/15548627.2019.1586250

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, T., Hu, E., Xu, S., Chen, M., Guo, P., Dai, Z., et al. (2021). ClusterProfiler 4.0: A Universal Enrichment Tool for Interpreting Omics Data. The Innovation 2, 100141. doi:10.1016/j.xinn.2021.100141

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, H., Liu, L., Li, W., Zou, D., Yu, J., Wang, L., et al. (2021). Transcription Factors in Colorectal Cancer: Molecular Mechanism and Therapeutic Implications. Oncogene 40, 1555–1569. doi:10.1038/s41388-020-01587-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Yaeger, R., Chatila, W. K., Lipsyc, M. D., Hechtman, J. F., Cercek, A., Sanchez-Vega, F., et al. (2018). Clinical Sequencing Defines the Genomic Landscape of Metastatic Colorectal Cancer. Cancer Cell 33, 125–136. doi:10.1016/j.ccell.2017.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, W., Soares, J., Greninger, P., Edelman, E. J., Lightfoot, H., Forbes, S., et al. (2013). Genomics of Drug Sensitivity in Cancer (GDSC): a Resource for Therapeutic Biomarker Discovery in Cancer Cells. Nucleic Acids Res. 41, D955–D961. doi:10.1093/nar/gks1111

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, X., Wu, H., Sheng, L., Liu, Y.-x., Ye, F., Wang, M., et al. (2019). Oncogenic Potential of Truncated RXRα during Colitis-Associated Colorectal Tumorigenesis by Promoting IL-6-STAT3 Signaling. Nat. Commun. 10, 1463. doi:10.1038/s41467-019-09375-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Yothers, G. (2012). Oxaliplatin in the Adjuvant Treatment of Colon Cancer. Ca-a Cancer J. Clin. 62, 3–4. doi:10.3322/caac.v62:3

CrossRef Full Text | Google Scholar

Keywords: colorectal cancer, C3, oxaliplatin resistance, prognosis, transcriptomic analysis

Citation: He X-S, Zou S-Y, Yao J-L, Yu W, Deng Z-Y, Wang J-R, Gan W-J, Wan S, Yang X-Q and Wu H (2021) Transcriptomic Analysis Identifies Complement Component 3 as a Potential Predictive Biomarker for Chemotherapy Resistance in Colorectal Cancer. Front. Mol. Biosci. 8:763652. doi: 10.3389/fmolb.2021.763652

Received: 24 August 2021; Accepted: 27 September 2021;
Published: 15 October 2021.

Edited by:

Deshui Jia, Shanghai General Hospital, China

Reviewed by:

Qian Zhao, Tianjin Medical University, China
Jiwei Zhang, Shanghai University of Traditional Chinese Medicine, China
Liqun Chen, Fuzhou University, China

Copyright © 2021 He, Zou, Yao, Yu, Deng, Wang, Gan, Wan, Yang and Wu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Hua Wu, wuhua@suda.edu.cn; Shan Wan, shanwan5@suda.edu.cn; Xiao-Qin Yang, yangxiaoqin@suda.edu.cn

These authors have contributed equally to this work

Download