1.
Introduction
The tumor microenvironment (TME) is a complex ecosystem comprising numerous cells, including immune cells, stromal cells and immunosuppressive cells [1]. Tumour endothelial cells (TECs) are a prominent component of TME and TECs, which have critical functions on tumor growth, progression and metastasis [2]. As tumor angiogenesis is essential for tumor growth and metastasis, inducing tumor-associated angiogenesis is a promising tactic in cancer progression [3]. Angiogenesis is now accepted as an important hallmark of cancer [4] and endothelial cells (which form tight adhesions to ensure vessel integrity) are essentials in inducing angiogenesis in TME [3]. In TME, TECs interact with tumor cells via juxtacrine and paracrine signaling during tumor intravasation and metastasis [5].
Colorectal cancer (CRC) is the fourth most common cancer and a leading cause of cancer mortality worldwide [6]. Recently stated that distinct stromal transcriptional and interactions are altered in colon cancer development [7] and stromal cells contribute to the CRC transcriptome [8]. It was also stated that stromal gene expression defines poor prognosis in colorectal cancer [9] and induces poor-prognosis subtypes in colorectal cancer [10]. Altogether, stromal cells have substantial influence and regulatory roles in colon cancer and endothelial cells (ECs) may cause colon carcinogenesis as a part of the stromal component.
So, analyzing the transcriptomes of colon TECs (CTECs) in TME reveals new targets and avenues and explores more uncontrolled factors in colon cancer research. We present a comprehensive analysis through bioinformatic tools and identified molecular and genetic alterations in TECs. We identified DEGs from microarray datasets of gene expression omnibus (GEO), hub genes from the interactions of deregulated genes and find out significant KEGG pathways. We also identified significant hub genes that are linked with poor survival of patients. Moreover, we found that hub genes are associated with immune cells infiltration and positively correlated with immune suppressive markers. This integrated analysis provides vital molecular insights into CTECs characterization, which may directly affect treatment recommendations for colon cancer patients. It provides opportunities for genome-guided clinical trials as well as drug development for the treating of colon cancer patients.
2.
Materials and methods
2.1. Datasets
We searched the NCBI gene expression omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) using the keywords "colon cancer endothelial cells, " "endothelial cells, " and "normal endothelial cells, " and identified one CTECs gene expression dataset GSE89287 (n = 24) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE89287) [11]. This dataset included different cell types (endothelial cells, macrophages and epithelial cells) which were isolated from non-paired primary normal colon tissues and colorectal carcinomas and subsequently, RNA was isolated. In this study, we included only endothelial cells and excluded the other cell types from our study. Endothelial cells were isolated using the facs-sorting method. The dataset included 24 samples of endothelial cells, including 10 normal colon samples and 14 colon tumor samples. The platform of this data is GPL4133 (Agilent-014850 whole human genome microarray 4x44K G4112F) with a feature number version. Gene expression profiling interactive analysis (GEPIA) (http://gepia2.cancer-pku.cn/#index) dataset (n = 316) was used for the study of prognostic and expression levels of hub genes [12]. In addition, the correlation between the expression of hub genes and tumor-infiltrating immune cells was explored via gene modules in the TIMER dataset (n = 457) (https://cistrome.shinyapps.io/timer/) [13]. For identifying the gene-gene correlation, we used the TCGA-COAD dataset (n = 287) (https://gdc-portal.nci.nih.gov/) which was normalized into base log2.
2.2. Identification of differentially expressed genes (DEGs) between CTECs and CNECs
We used the web tool Network Analyst [14]to identify the significant DEGs between CTECs and CNECs. Dataset was normalized by quantile normalization and the R package "limma" was utilized to identify the DEGs between CTECs and CNECs. We utilized the Network Analyst tool to identify the average expression level of single gene having multiple probes in this gene expression study. We identified the DEGs with a threshold of |LogFC| > 0.585 and adjusted P value < 0.05.
2.3. Gene-set enrichment analysis
We performed gene-set enrichment analysis of the DEGs by GSEA [15]. We identified the KEGG [16] pathways that were significantly associated with the up-regulated and the down-regulated DEGs, respectively. We identified the significant pathways with a threshold of false discovery rate (FDR) < 0.05.
2.4. Construction of protein-protein interaction (PPI) network and modular analysis of PPI
To better know the relationship among these screened DEGs, the PPI network was established using the STRING database [17]. The hub genes in the PPI network were identified according to a degree using the node explorer module of Network Analyst software [14]. A hub gene was defined as a gene that was connected to no less than 10 other DEGs. A cytoscape plug-in molecular complex detection (MCODE) was employed to detect the modules from the PPI network [18]. We identified the significant modules based on the MCODE score and node number. The threshold of the MCODE was Node Score Cut-off: 0.2, Haircut: true, K-Core: 2 and maximum depth from Seed: 100.
2.5. Survival analysis of hub genes
We compared the overall survival (OS) and the disease-free survival (DFS) of colon cancer patients classified based on gene expression levels (expression levels > median versus expression levels < median). Kaplan-Meier survival curves were used to show the survival differences and the log-rank test was utilized to evaluate the significance of survival differences. The prognostic roles of screened hub genes in COAD were analyzed using expression profiling interactive analysis (GEPIA) databases [12]. Cox P < 0.05 was considered as significant between two groups of patients.
2.6. Validation of prognosis-associated hub gene expression levels in TCGA COAD dataset
The expression levels of top hub genes were further validated using the gene expression profiling interactive analysis (GEPIA), a newly developed interactive web server for analyzing the RNA sequencing data [12]. We used TCGA COAD tumor data with matched normal. Hub genes with |log2FC| > 0.585 and P < 0.05 were considered statistically significant between the two groups.
2.7. Evaluation of immune scores and stromal scores in colon cancer
We utilized the "ESTIMATE" R package to quantify the immune scores (predict the immune cell quantity) and stromal scores (predict the stromal cells quantity) for each of the tumor samples [19]. Then we calculated Spearman's correlations between the expression levels of prognostic hub genes and immune and stromal scores. The threshold value of correlation is R > 0.30 and P-value is less than 0.001 (Spearman's correlation test).
2.8. Analysis of immune cell infiltration and correlation of top hub genes with immune-inhibitory markers
We analyzed the significant correlations of hub genes with the abundance of immune infiltrates, including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages and dendritic cells, via gene modules in TIMER [13]. In addition, Spearman's correlations between the expression level of hub genes and the immune inhibitory marker genes of tumor-infiltrating immune cells were explored via correlation modules in TIMER [13]. Moreover, we identified pearson's correlations between the expression level of hub genes and the marker genes of inhibitory immune cells in TCGA-COAD datasets using R software. The R package "ggplot2" was employed for drawing a correlation graph [20]. The gene markers of tumor-infiltrating immune cells included markers of monocytes, TAMs, M2 macrophages, T-helper 1 (Th1) cells, Tregs and exhausted T cells. These gene markers are illustrated in prior studies [12,21,22,23,24].
2.9. Statistical analysis
We used R programming software for the statistical analysis of this study. For the differential expression analysis of gene expression data, only genes with |logFC| > 0.585 and adjusted P-value < 0.05 were considered statistically significant between the two groups. In the survival analysis, Cox P < 0.05 was considered as statistically significant. Spearman's correlation test between the ssGSEA scores and the expression level of specific prognostic genes was performed because these data were not normally distributed (P < 0.05) [25]. For analyzing the correlations between the expression levels of hub genes with the expression levels of other marker genes, we employed Pearson's correlation test because these data were normally distributed [25].
3.
Results
3.1. Identification of Differentially Expressed Genes
Based on the log2FC and adjusted P-value, we identified 362 differentially expressed genes (DEGs) in CTECs relative to the CNECs, which included 117 up-regulated genes (Supplementary Table S1) and 245 down-regulated (Supplementary Table S2) in the CTECs samples. Table 1 demonstrates the top ten up-regulated genes (SPARC, IGFBP7, COL4A2, SELE, CXCL1, VWF, CCDC3, FSTL1, VWA1 and IFITM3) with significant statistical description. The most up-regulated gene SPARC is up-expressed in the stromal portion of CRC tissues [26]. The extracellular matrix-associated gene, COL4A2, is linked with cancer stemness [27]. Another up-regulated gene, CXCL1, is associated with increasing the metastatic ability of colon cancer by enhancing cell migration, matrix metalloproteinase-7 expression and EMT [28]. VWF, another up-regulated gene, is a plasma marker for the early detection of adenoma and colon cancer [29]. Table 2 illustrated the top ten down-regulated genes (IGKV1-5, IGLV1-44, IGKC, IGHV3-30, KRT81, IGHG1, IGHM, POU2AF1, IGKV2-24 and IGLJ3) with all significant statistical descriptions. Interestingly, most of the top down-regulated genes are associated with immunological activities. Some other down-regulated genes of CTECs have been demonstrated to be under-expressed in CRC. For example, the expression of POU2AF1 is down-regulated in the rat colon and in turn, reveals proximal-distal differences in the process of histone modifications and proto-oncogene expression [30]. Altogether, abnormally expressed numerous genes in CTECs compared to NTECs identified by the bioinformatic analysis have been associated with CRC pathology and carcinogenesis.
3.2. KEGG pathway enrichment analysis of significantly identified DEGs
The functional enrichment analysis identifying the enriched up-regulated and down-regulated biological pathways that are associated with significant DEGs (Figure 1). GSEA pathway analysis revealed nine up-regulated pathways (focal adhesion, ECM-receptor interaction, pathways in cancer, small cell lung cancer, adherens junction, arrhythmogenic right ventricular cardiomyopathy (ARVC), pathogenic escherichia coli infection, cytokine-cytokine receptor interaction and proximal tubule bicarbonate reclamation) (Figure 1A) in CTECs. In addition, we also identified ten down-regulated pathways (protein export, cell adhesion molecules (CAMs), arrhythmogenic right ventricular cardiomyopathy (ARVC), N-Glycan biosynthesis, the intestinal immune network for IgA production, primary immunodeficiency, hypertrophic cardiomyopathy (HCM), ECM-receptor interaction, dilated cardiomyopathy and vasopressin-regulated water reabsorption) in CTECs (Figure 1B). Recently, it was stated that focal adhesion, ECM-receptor interaction, pathways in cancer, small cell lung cancer, adherens junction and cytokine-cytokine receptor interaction pathways are dysregulated in colon tumors stroma [7].
3.3. Construction of PPI by using significant DEGs and identification of hub genes and functional clusters from the PPI
Using all 362 DEGs, PPI networks were constructed using STRING software and hub genes were identified using the node explorer of network analyst software. Many significant DEGs are involved in the networks (Supplementary Table S3). The identification of hub genes explored that FN1, GAPDH, COL1A1, COL1A2, CTNNB1, LAMB1, LAMC1, SPARC, IGFBP3 and COL4A1 are top up-regulated hub genes and HSP90B1, PDIA6, ITGA4, DDOST, SDC1, CKAP4, CD38, MANF, ITGB7 and CCR2 are top down-regulated hub genes in CTECs respectively (Figure 2 and Table 3). Interestingly, we found that top-up-regulated (Table 1) SPARC also acts as a hub gene in the PPI network. FN1 suppressed apoptosis and promoted viability, invasion and migration of tumor cells in CRC [31]. It was showed that COL1A1 and COL1A2 are associated with colon cancer [32,33]. The high expression levels of CTNNB1 positively correlated with metastasis of colon cancer [34]. Deregulation of SDC1 contributes to cancer progression by promoting cell proliferation, metastasis, invasion and angiogenesis [35]. Altogether, it may be stated that CTECs-derived these hub genes may contribute to colon cancer pathogenesis.
Moreover, MCODE identified 15 clusters from the original PPI networks. The description of MCODE derived clusters is illustrated in Table 4. The top significant cluster 1 contained 15 nodes and 105 edges (Figure 2 and Table 4). We identified the functional enrichment of KEGG pathways for all clusters by using the GSEA. Interestingly, we found that eight of the clusters out of 15 are associated with the enrichment of KEGG pathways. Gene set of Cluster 1 is related to the enrichment of 4 pathways: ECM-receptor interaction, small cell lung cancer, pathways in cancer and focal adhesion (Figure 2A, B). Cluster 3 (Figure 2C, D), Cluster 4 (Figure 2E, F) and Cluster 9 (Cell adhesion molecules and natural killer cell-mediated cytotoxicity) are associated with the enrichment of immune, stromal and cancer-associated pathways. Some of the enriched pathways in Cluster 3 included ECM-receptor interaction, focal adhesion, small cell lung cancer, cell adhesion molecules, the intestinal immune network for IgA production, pathways in cancer and hematopoietic cell lineage. In addition, some of the pathways enriched in Cluster 4 are the chemokine signaling pathway, adherens junction, cytokine-cytokine receptor interaction, focal adhesion, NOD-like receptor signaling pathway, pathways in cancer, leukocyte trans-endothelial migration and tight junction. Protein export is enriched in Cluster 2 and the B cell receptor signaling pathway is enriched in Cluster 5. Cluster 13 (Fatty acid metabolism and PPAR signaling pathway) and cluster 15 (Glycosphingolipid biosynthesis-Lacto and neolacto series and N-Glycan biosynthesis) are mainly associated with the enrichment of metabolic pathways. Uddin et al. also identified some of these pathways, including pathways in cancer, small cell lung cancer, ECM-receptor interaction, focal adhesion, natural killer cell-mediated cytotoxicity and cell adhesion molecules, are enriched in colon tumor stroma [7,9]. Altogether, it can be stated that this hub PPI from CTECs may contribute to colon carcinogenesis.
3.4. Hub genes are negatively associated with prognosis
We considered the top ten up-regulated and top ten down-regulated genes from the original PPI network for survival analysis. We selected the TCGA-COAD database for screening OS and DFS. We got COL1A1, COL1A2, IGFBP3, SPARC and DDOST hub genes are significantly involved with patient survival time (Figure 3). We found that the higher expression level of up-regulated COL1A1, COL1A2, IGFBP3 and SPARC are worse for patient's DFS and COL1A2 is also associated with poor OS. In addition, the low expression level of down-regulated DDOST is more inferior for patient's OS and DFS. It showed that higher expressions of COL1A1 were related to poor survival in colon cancer patients [33]. COL1A2 gene is also associated with the prognosis of IIA stage colon cancer [32]. Stromal SPARC had a pro-metastatic impact in vitro and was a characteristic of aggressive tumors with a poor prognosis in CRC patients [26]. Collectively, hub genes of CTECs are prognostic markers in COAD and may be contributed to colon carcinogenesis.
3.5. Comparisons of the expression levels of prognostic hub genes between colon cancer and normal tissue
Interestingly, we found that all four up-regulated hub genes (COL1A1, COL1A2, IGFBP3 and SPARC) were consistently up-regulated in TCGA COAD samples versus normal samples (P < 0.05) (Figure 4). It indicates that CTECs may have a contribution to the COAD-derived transcriptomes in colon cancer. However, one of the down-regulated hub genes, DDOST, showed no significant expression differences between TCGA-COAD samples versus normal samples (P < 0.05), suggesting that this hub gene is expressed explicitly in CTECs.
3.6. Prognostic hub genes are correlated with immune infiltrations and associated with immune-inhibitory markers in colon TME
Tumor-infiltrating lymphocytes are independent prognostic factors for better survival and sentinel lymph node status in cancers [36]. Genes highly expressed in the microenvironment are expected to negatively associate with tumor purity, while the opposite is expected for genes highly expressed in the tumor cells [13]. So, we assessed whether the expression of hub genes was correlated with immune infiltration levels in COAD. We investigated the correlations of prognostic hub genes (COL1A1, COL1A2, IGFBP3, SPARC and DDOST) with immune scores and stromal scores in COAD. Interestingly, we found that the expression levels of COL1A1, COL1A2, IGFBP3 and SPARC positively correlated with the immune scores (Figure 5A) and stromal scores (Figure 5B) (Spearman's correlation test, R > 0.30, P < 2.2e−16). This result indicated that the prognostic hub genes are associated with the regulation of the tumor microenvironment.
In addition, we identified the correlations between the expression of hub genes and the immune infiltration levels in COAD by using the TIMER tool. The results discovered that up-regulated COL1A1, COL1A2, IGFBP3 and SPARC expression have significant (P ≤ 0.05) correlations with tumor purity in COAD but down-regulated DDOST not correlated with tumor purity (Figure 6). In addition, COL1A1 expression has significant correlations with infiltrating levels of CD4+ T cells, macrophage, neutrophil and dendritic cells. We also found that a significant positive correlation between COL1A2, IGFBP3 and SPARC expression and infiltration of all five immune cells (CD4+ T cells, CD8+ T cells, Neutrophils, Macrophages and Dendritic cells) in COAD. Besides, we found significant negative correlations between DDOST expression and infiltration of all five immune cells (CD4+ T cells, CD8+ T cells, Neutrophils, Macrophages and Dendritic cells) (Figure 6). Recently it was stated that macrophages and neutrophils are associated with immunosuppressive inflammation to modulate anti-tumor immunity [37]. DCs can promote tumor metastasis by increasing Treg cells and reducing CD8+ T cell cytotoxicity [38]. These results suggest that the expression of COL1A1, COL1A2, IGFBP3, SPARC and DDOST has potential roles in tumor purity and immune cell infiltrations in colon cancer.
To elucidate the relationship between up-regulated hub genes COL1A1, COL1A2, IGFBP3 and SPARC and the diverse immune infiltrating cells, we find out the correlations between these genes and immune marker sets of various immune cells including monocytes, TAMs, M2 macrophages, Th1, Tregs and T cell exhaustion (Table 5). Surprisingly, we got that the hub genes COL1A1, COL1A2, IGFBP3 and SPARC are positively correlated (Pearson correlation test, P ≤ 0.001) with the immune markers of monocytes, TAMs, M2 macrophages, Th1, Tregs and T cell exhaustion. These findings suggest that the high expression level of hub genes plays an essential role in the infiltration of monocytes, TAMs, M2 macrophages, Th1, Tregs and T cell exhaustion in COAD. It was recently reported that LAYN acts as a prognostic biomarker for determining prognosis and immune infiltration in gastric and colon cancers [21]. Potent immunosuppressive T-regulatory cells (Tregs) are found in a vast array of tumor types and tumor-infiltrating Tregs are often associated with poor clinical outcomes. Tregs also promote cancer progression through their ability to limit anti-tumor immunity and promote angiogenesis [40]. Another marker, FOXP3, plays an important role in Treg cells which leads to the suppression of cytotoxic T cells attacking tumor cells [40]. TIM-3 is a crucial surface protein on exhausted T cells, a critical gene that regulating T cell exhaustion [41].
Since the elevated expression of prognostic hub genes were positively associated with the immune inhibitory markers PD-1 and TGFB1 (Table 5), we expect the expression levels of PD-L1 (CD274), PD-L2 (PDCD1LG2) and TGFBR1 could be positively associated with the prognosis-associated hub genes. Interestingly, we found that PDL1 (CD274), PDL2 (PDCD1LG2) and TGFBR1 are moderate to highly correlated with all four-prognosis related highly expressed hub genes: COL1A1, COL1A2, SPARC and IGFBP3 (Figure 7). T cell immunity recovered by PD-1/PD-L1 immune checkpoint blockade has been demonstrated to be a promising cancer therapeutic strategy [42]. The expression of PD-L2 was observed in tumors, stroma and endothelial cells and PD-L2 expression is relevant to anti-PD-1 therapy in cancer [43]. TGFBR1 is one of the receptors for TGF-β ligands. The TGF-β signaling pathway is crucially associated with stimulation, induction and maintenance of EMT in cancers [44]. TGFBR1 is up-regulated in multiple malignancies and dysregulation of TGFBR1 leads to tumorigenesis by controlling cellular signaling, which in turn is associated with colorectal cancer risk [45]. Altogether, it suggests that the prognostic associated hub genes may regulate the immune-suppressive activities of TME through the interactions with the immune marker in colon cancer.
4.
Discussion
TECs are the substantial component of tumor stroma in TME and these cells are associated with tumor malignancies, progression and metastasis [5]. This present study conducted a differential expression analysis using gene expression profiling data from the GEO database. We identified 362 DEGs in CTECs, including 117 up-regulated genes (SPARC, COL4A2, CXCL1 and VWF) and 245 down-regulated genes. Interestingly, we found that many immunological genes are down-regulated in CTECs (Table 2). Subsequent pathway enrichment analysis revealed that up-regulated DEGs were significantly enriched with cancer, cellular development and immune regulation. Down-regulated DEGs were enriched considerably with mostly immunodeficiency and metabolism-related pathways (Figure 1). These results revealed the abnormal cellular growth, immune regulatory, cancerous and metabolism-associated pathways in CTECs.
Next, we employed DEGs to construct a PPI network and extracted significant clusters from the original PPI network (Figure 2). Finally, we identified hub genes (degree > 10) that are dysregulated in CTECs. We identified CTECs-derived five hub genes (COL1A1, COL1A2, SPARC, IGFBP3 and DDOST) whose expression was significantly associated with the poor prognosis in TCGA-COAD data (Figure 3). These hub genes are mainly involved in protein digestion and absorption (COL1A1 and COL1A2), cellular signaling (SPARC and IGFBP3) and enzymatic action of metabolism (DDOST), suggesting that the deregulation of these cellular functions in CTECs may contribute to the altered prognosis in colon cancer. In addition, the expression analysis of these prognostic hub genes in colon cancer was further evaluated using the GEPIA database [12]. We found that four hub genes (COL1A1, COL1A2, SPARC and IGFBP3) are also significantly up-regulated in TCGA-COAD samples. Still, another down-regulated hub gene (DDOST) is not considerably down-regulated in TCGA-COAD samples (Figure 4). It indicates that CTECs-associated gene signature substantially contributed to colon carcinogenesis.
Furthermore, our study showed that the immune infiltration levels and diverse immune marker sets are correlated with the expression levels of COL1A1, COL1A2, SPARC and IGFBP3 hub genes. Our analysis demonstrated that a significant positive correlation between the COL1A1, COL1A2, SPARC and IGFBP3 expression levels and infiltration level of CD4+ T cells, macrophages, neutrophils and DCs in COAD (Figure 6). These results are suggesting that the potential immunological roles of hub genes in colon TME. In addition, our study indicated that COL1A1, COL1A2, SPARC and IGFBP3 could activate Tregs and induce T cell exhaustion. The increase in up-regulated four hub genes expression positively correlates with Treg and T cell exhaustion markers (Table 5). Other markers of monocytes, TAM, M2 macrophage, Th1, are also positively correlated with the expression of these hub genes (Table 5). Some other crucial immune inhibitory markers, PD-L1, PD-L2 and TGFBR1, are significantly correlated with the expression level of COL1A1, COL1A2, SPARC and IGFBP3 (Figure 7). Altogether, the prognostic hub genes are associated with immunosuppressive colon TME.
Since our analysis identified numerous CTECs-derived transcriptomes that are critically associated with colon cancer, we speculated that the endothelial cells might contribute to the other vascular pathologies, including the central nervous system (CNS) and other diseases. It was stated that cerebral endothelial cells play an active role in the pathogenesis of CNS inflammatory diseases [46]. Endothelial cells are associated with the leukocyte migration in CNS and regulating the leukocyte-endothelial cell interactions and the crosstalk between endothelial cells and glial cells or platelets in CNS [46]. The major vascular anomalies of the nervous system included cerebral cavernous malformations (CCMs) and arteriovenous malformations (AVMs) [47,48]. These CNS-associated vascular diseases are characterized by genetic alterations, such as gene polymorphisms and mutations [47,48,49,50,51]. For example, germline mutation enrichment in pathways controlling the endothelial cell homeostasis in patients with brain arteriovenous malformation [48]. Altogether, our findings may provide insights to find the roles of endothelial cells in diseases, including cancer, CNS-associated disease and other vascular pathologies.
This study identified numerous deregulated transcriptomes in the CTECs that could be used as biomarkers for the diagnosis and prognosis of colon cancer and may provide therapeutic targets for colon cancer. However, to translate these findings into clinical application, further experimental and clinical validation would be necessary.
5.
Conclusions
In summary, we identified potential CTECs-derived significantly deregulated transcriptomes that are involving with the pathogenesis, prognosis and immune inhibition within the tumor microenvironment of colon cancer. This study reveals a potential regulatory mechanism of CTECs in the tumor microenvironment of colon cancer and may contribute to revealing CTECs-colon cancer cellular cross-talk.
Conflicts of interest
The authors declare no conflict of interest.
Acknowledgments
We thank all laboratory staff and technicians for their participation in this study. The present study was supported by the project from Xinjiang Health Poverty Alleviation Work Appropriate Technology Promotion (Grant no. SYTG-Y202035).
Author contributions
Jie Wang and Md. Nazim Uddin contributed equally to this work. Jie Wang and Md. Nazim Uddin conceived the research, designed analysis strategies. Rehana Akter wrote the manuscript and Yun Wu conceived the study and wrote the manuscript. All the authors read and approved the final manuscript.
DownLoad:
