The GSE41657 (analyzed on Agilent-014850 Whole Human Genome Microarray 4^x44K G4112F platform; conducted by Peking Union Medical College and Chinese Academy of Medical Sciences in China) and GSE41258 (analyzed on Affymetrix Human Genome U133A Array platform; conducted by Weizmann Institute of Science in Israel) gene expression profiles in human colorectal cancer were downloaded from the GEO database of the National Center for Biotechnology Information (www.ncbi.nlm.nih.gov/geo). The patients selected for the cancer and normal groups were all aged <50 years.

Firstly, GEO2R was used to identify the differentially expressed genes (DEGs) between youth colorectal cancer and normal tissues. GEO2R (www.ncbi.nlm.nih.gov/geo/geo2r/) is a convenient web tool for DEG screening by comparing two groups of samples in a GEO dataset (11). An adjusted P<0.05 and |log₂ fold change (FC)|≥1 were set as the cut-off criteria.

Secondly, these DEGs were classified according to Gene Ontology (GO; http://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG; www.genome.jp/kegg) pathways using Database for Annotation, Visualization and Integration Discovery (DAVID; http://david.ncifcrf.gov/) software. Based on the enriched GO terms and significant KEGG pathways, the disease-associated pathways and genes were screened.

Thirdly, DEGs were imported into STRING software (version 3.6.0; http://string-db.org/) to produce protein-protein interaction networks, and then Cytoscape software (version 3.6.0; http://cytoscape.org/) was used to construct a visible network diagram. In this network, cytoHubba plugin was used to screen the top 10 hub genes using degree algorithm. Combined with the analysis of signaling pathways, disease-associated genes were screened out.

TCGA (http://cancergenome.nih.gov/) dataset was used to validate the disease-associated genes screened by GEO datasets. TCGA dataset was analyzed using Gene Expression Profiling Interactive Analysis (http://gepia.cancer-pku.cn/), which is a commonly used interactive website that plots expression profiles of selected genes. A total of 362 colorectal cancer and 51 normal tissues from TCGA database were selected for analysis. Survival analysis was further performed on these genes, based on gene expression levels, as previously described (12). Key genes were obtained by overall survival (OS) and disease-free survival (DFS) analysis, and the association between the expression of these genes and tumor staging was explored. Patients with high expression levels (50%) of key genes were assigned to high expression group. The expression levels of key genes were also validated using the HPA database (http://www.proteinatlas.org/) (13).

Six pairs of youth colorectal cancer and adjacent normal tissues from patients aged <50 years (age range, 35–48; 3 men and 3 women) were collected from the First Department of Gastrointestinal Surgery of the Affiliated Hospital of Putian University (Putian, China). All colorectal cancer and normal tissues were stored in a liquid nitrogen tank within 30 min from their removal from the patient. For RNA extraction, the sample was ground into pieces in a mortar filled with liquid nitrogen. Total RNA was then extracted using RNAiso Plus (Takara Bio, Inc.) and reverse transcribed into cDNA using PrimeScript™ RT reagent kit with gDNA Eraser (Takara Bio, Inc.), according to the manufacturer's instructions. RT-qPCR was performed using the SYBR Premix Ex Taq™ II (Takara Bio, Inc.), according to the manufacturer's instructions. The PCR primers were designed by the National Center for Biotechnology Information (www.ncbi.nlm.nih.gov/) and provided by SangonBiotech Co., Ltd. (Table SI). The gene GAPDH was selected as an internal reference, to compare gene expression in different samples. The RT-qPCR process included the following steps: Pre-denaturation at 95°C for 30 sec, followed by 40 cycles of denaturation at 95°C for 5 sec and extension at 60°C for 35 sec (7500 Real Time PCR System; Thermo Fisher Scientific, Inc.). Each sample was tested in triplicate. Gene expression values in each sample were calculated using the 2^−∆∆Cq method (14). The present study was approved by the Ethics Committee of the Affiliated Hospital of Putian University (approval no. 202016). Written informed consent was obtained from each patient for sample collection and analysis.

SPSS version 17.0 (SPSS Inc.) was used to analyze the data. Continuous variables are presented as the mean ± standard deviation. One-way ANOVA with least significant difference post hoc test was used to analyze the data of more than two groups. Student's t-test was used for comparisons between two groups. A paired t-test was used when comparing the six pairs of youth colorectal cancer and adjacent normal tissues, and an unpaired t-test was used for other comparisons. If the variance was not equal between two groups, Mann-Whitney U test was used for statistical analysis. The statistical significance of survival time was determined by the log-rank test. P<0.05 was considered to indicate a statistically significant difference.

The gene expression levels of genes were downloaded from the GEO database. Based on GEO2R analysis, 6,322 and 12,585 DEGs in GSE41657 and GSE41258 were screened in youth colorectal cancer compared with normal colorectal tissues. Screened with the cut-off criteria of adjusted P<0.05 and |log₂ FC|>1, 1,256 upregulated and 1,782 downregulated genes were identified in the GSE41657 dataset (Fig. S1A), and 311 upregulated and 568 downregulated genes were identified in the GSE41258 dataset (Fig. S1B). In the integrated analysis of the two datasets, 443 overlapping genes, including 131 upregulated and 312 downregulated, were obtained (Fig. S1C and D).

The functional enrichment of 443 candidate genes was analyzed using DAVID. Three GO category results are presented in Table SII, including biological processes, cellular components and molecular functions. The biological process results revealed that the selected genes were mainly enriched in ‘cellular response to zinc ion’, ‘negative regulation of growth’ and ‘cell proliferation’. The cellular component results revealed that selected genes were mainly enriched in ‘extracellular exosome’, ‘extracellular space’ and ‘apical plasma membrane’. The molecular function analysis revealed that the selected genes were mainly enriched in ‘protein binding’, ‘carbonate dehydratase activity’ and ‘CXCR chemokine receptor binding’. These results demonstrated that the majority of the candidate genes were significantly enriched in ‘binding’ and ‘cell proliferation’.

The signaling pathway enrichment of 443 candidate genes was analyzed using KEGG pathway online databases. A total of 23 genes were significantly enriched in ‘Pathways in cancer’ (Table SIII). Based on previous reports (15,16), ‘Pathways in cancer’ is an important signaling pathway associated with the occurrence and development of cancer.

A total of 443 candidate genes were analyzed using the STRING 11.0 database, and then the protein-protein interaction (PPI) was imported into Cytoscape 3.6.0 software to build a visible network diagram. A total of 390 nodes and 1,462 edges are presented in the network (Fig. S2).

Within the PPI network, the cytoHubba plugin was used to screen for hub genes. Based on the degree algorithm, the top 10 hub genes were VEGFA, CXCL8, MYC, CD44, CXCL12, CCND1, IGF1, CXCL1, KIT and SOX9 (Table I). Combining signaling pathway analysis, CXCL8, IGF1, KIT, CXCL12, CCND1, VEGFA and MYC in the ten hub genes were enriched in the ‘pathways in cancer’ (Fig. S3). Among these genes, CXCL8, CCND1, VEGFA and MYC were upregulated, and the others downregulated.

Differential expression of CXCL8, IGF1, KIT, CXCL12, CCND1, VEGFA and MYC in colorectal cancer compared with normal colorectal tissues, was screened in TCGA database. The expression levels of these eight genes were significantly higher in both colon and rectal cancer tissues, as compared with normal colorectal tissues (P<0.05; Fig. 1).

These genes were further subjected to survival analysis to screen out genes associated with colorectal cancer prognosis. The expression level of CXCL8 was found to impact overall survival (OS; P<0.05; Fig. 2A), and that of VEGFA to impact disease-free survival (DFS; P<0.05; Fig. 2D). The expression level of MYC and CD44 was found to have no impact on OS or DFS (P>0.05; Fig. 2E-H). Thus, CXCL8 and VEGFA were identified as key genes.

Compared with the expression at different Tumor-Node-Metastasis (TNM) stages of colorectal cancer, it was found that the expression of the CXCL8 gene in high TNM-stage colorectal cancer was increased, but there was no statistical significance (P>0.05; Fig. 2I), while the gene expression of VEGFA in high TNM stages was significantly higher (P=0.000661; Fig. 2J).

The comparison of CXCL8 and VEGFA expression between youth and elderly colorectal cancer, revealed no significant differences (Fig. 3A and B). Further analysis indicated that high CXCL8 and VEGFA expression was associated with tumor stage. The proportion of stages III–IV was 56.25% in youth colorectal cancer, 42.95% in elderly colorectal cancer with high CXCL8 expression, and 65.63% in youth colorectal cancer, 47.65% in elderly colorectal cancer with high VEGFA expression (Fig. 3F). The HPA database was also used to verify the protein expression levels of key genes. The expression of CXCL8 and VEGFA in colorectal cancer from the HPA database was also consistent with the aforementioned bioinformatics analysis (Fig. 4C and D).

To confirm the reliability of the findings from data mining, the two genes, CXCL8 and VEGFA, were selected for validation by RT-qPCR in six pairs of youth colorectal cancer and adjacent normal tissues. According to the experimental results, the difference in expression of these two genes is consistent with the results of data mining. In youth colorectal cancer, the expression level of CXCL8 and VEGFA were significantly upregulated (P<0.05; Fig. 4A and B).