Forkhead box D subfamily genes in colorectal cancer: potential biomarkers and therapeutic targets

Data download and pre-processing

Using the TCGAbiolinks package and R software (version 4.0.4) (Colaprico et al., 2016; R Core Team, 2021), we downloaded the following data from The Cancer Genome Atlas (TCGA) database: HTSeq-counts, clinical, and single-nucleotide polymorphisms data for TCGA-Colon Adenocarcinoma Collection (COAD) and TCGA-Rectum Adenocarcinoma (READ). We excluded samples with null values and no survival information, and 643 samples (the baseline data of patients was detailed in Table S1) were included in the study. All the samples in the dataset were derived from Homo sapiens using the Illumina platform. The count data were converted into transcripts per million reads format and log2 conversion was performed. Additionally, we downloaded the gene expression matrix and annotation files of normal tissues from the Genotype-Tissue Expression (GTEx) database. In the bioinformatics analysis, we used strict parameter control. We adopted the Benjamini–Hochberg method to adjust the p-value to control the false discovery rate (FDR) and set adj.p < 0.05 and |log2 fold change| > 1 as the screening threshold.

FOXD subfamily identifications and characteristic analysis

The combined data of TCGA-COAD and TCGA-READ Counts were analyzed by the DESeq2 package (Love, Huber & Anders, 2014) to obtain differentially expressed genes (DEGs). Then, we performed batch survival analysis using the Survival and Survminer packages in R to obtain genes with different expression levels. The FOXD subfamily (FOXD1, FOXD2, FOXD3, and FOXD4) was defined as the primary variable based on the intersection of the above two gene sets. The FOXD subfamily receiver operating characteristic curve (ROC) was analyzed in the pROC package and plotted in the ggplot2 package. In addition, lollipop plots of amino acid changes (gene mutations) in the FOXD subfamily were drawn the using maftools package (Mayakonda et al., 2018).

Functional enrichment analysis of FOXD subfamily

DESeq2 package was used to screen FOXD subfamily-related DEGs with high or low expression. The ggplot2 package was used to draw the volcano maps and the pheatmap package (https://CRAN.R-project.org/package=pheatmap) was used to draw a heat map showing the overall expression pattern of FOXD subfamily-related DEGs. The selection criteria were as follows: p < 0.05 and |log2FC|>1.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the clusterProfiler package (Yu et al., 2012), and the significance criterion was set at p < 0.05. For Gene Set Enrichment Analysis (GSEA), we choose “c2. Cp. Kegg. V7.0. Entrez. GMT” as a reference gene set, and the significant concentration threshold values were as follows: FDR < 0.25 and p < 0.05.

Protein–protein interaction analysis and hub gene screening

Protein–protein interaction (PPI) network analysis was performed using the STRING database and Cytoscape software (version 3.8.2). Hub genes were identified using the Molecular Complex Detection (MCODE) plugin (Bader & Hogue, 2003).

Immune cell infiltration analysis

Based on the principle of linear support vector regression, CIBERSORT was used to deconvolve the transcriptome expression matrix and estimate the composition and abundance of 22 immune cells in a mixed cell sample, as previously described (Newman et al., 2015). Samples with p < 0.05 were filtered to generate an immune cell infiltration matrix. Correlation heat maps were drawn using R-package corrplot (Yang et al., 2021) to visualize the interrelationship between the 22 immune cells, and lollipop maps were drawn using the GSVA package (Hänzelmann, Castelo & Guinney, 2013) to show the correlation of 22 immune cells with the FOXD gene family.

Survival analysis and clinical relevance

Prognostic survival analysis was performed using the Survminer and Survival packages. To determine the thresholds for the high and low expression groups, raw gene expression data were collected and normalized using log2 transformation to reduce data discretization and median centralization. Next, we determined the tangency point of the greatest difference by calculating the correlation between patient survival data and gene expression levels using the ‘surv_cutpoint’ function (R package ‘survminer’), which divides the expression data based on minimum-tangency maximization statistics. To ensure the reliability of results, the robustness of the selected tangent points was verified using various methods such as resampling and cross-validation. The Kaplan–Meier method was used to draw the survival curve, and the log-rank test was used to compare the survival differences between the high and low expression groups. In addition, multivariate Cox proportional hazards model analyses included age, sex, and other clinical characteristics as covariates to adjust for confounders (The detailed parameters are shown in Table S2). To plot a concordance curve, we fitted the nomogram and performed calibration analysis.

Immunohistochemistry analysis

CRC tumor tissues and adjacent normal colon tissues (n = 28) were collected from the Second Hospital of Dalian Medical University. All participants were informed of the reasons to conduct the study, anonymity was ensured, and how the collected data were stored. We received written informed consent from participants of our study. Eight cases with missing information and unqualified section quality were excluded, and 20 samples were retained. The protocol was approved by the Second Hospital of Dalian Medical University ethics committee (approval number: 2021-079). FOXD subfamily gene expression in the collected samples was measured using IHC. Human CRC and adjacent normal colon tissues were fixed overnight in 10% formalin and embedded in paraffin. The tissue blocks were cut into 3 μm tissue sections. The tissue sections were dewaxed, rehydrated, and subjected to antigen retrieval in either citrate or EDTA buffer, depending on the antibody. Following 3% hydrogen peroxide treatment, the sections were incubated overnight at 4 °C with the following primary antibodies: anti-FOXD1 (1:100 dilution, ab129324; Abcam, Cambridge, UK), anti-FOXD2 (1:50; TA322739, OriGene, Rockville, Maryland, US), anti-FOXD3 (1:100 dilution, abs133773; Absin, Shanghai, China), and anti-FOXD4 (1:100, orb156929; Biorbyt, Cambridge, UK). The bound antibodies were detected using horseradish peroxide-conjugated secondary antibodies and developed with nickel-diaminobenzidine using a two-step test kit (PV-9000; Beijing Zhongshan Jinqiao Biological, Beijing, China), according to the manufacturer’s instructions. IHC scores were calculated through multiplying the staining score by the percentage score. The percentage of positive cells were judged as 0, ≤5% positive; 1, 6% to 25% positive; 2, 26% to 50% positive; 3, 51% to 75% positive; and 4, 76% to 100% positive. Staining intensity was judged as 0 for no staining, 1 for light yellow or yellow, 2 for brown, and 3 for dark brown (Gao et al., 2020; Khan et al., 2018; Mo et al., 2021). The IHC staining results were assessed by one senior pathologist. The final staining score ranged from 0 to 12, less than six was considered as low expression, while staining score of six or more was considered as high expression.

Cell culture

CCD841CON cells were obtained from Meisen CTCC (Zhejiang, China) and cultured in 90% Dulbecco’s Modified Eagle Medium (DMEM) supplemented with 1% penicillin/streptomycin (PS) and 10% fetal bovine serum (FBS). All colon cancer cell lines (SW480, SW620, HT29 and HCT 116) were purchased from ProCell company (Wuhan, China). Among of them, SW480 and SW620 cells were cultured in 90% DMEM supplemented with 10% FBS and 1% PS. HT29 and HCT116 cells were cultured in 90% McCoy’s 5A medium supplemented with 10% FBS with 1% PS. The cells were incubated at 37 °C in a humidified 5% CO₂ incubator.

Western blot assay

Total proteins were extracted from cells using radioimmunoprecipitation assay. The concentrations of the protein samples were determined using a bicinchoninic acid protein concentration assay. The prepared 10% sodium dodecyl sulfate-polyacrylamide gel was poured onto a glass plate and allowed to solidify. Next, 50 µg of the marker and each protein sample was loaded into the gel wells. Subsequently, the samples were electrophoresed under reducing conditions and transferred onto polyvinylidene fluoride membranes. The membranes were blocked for 2 h. Next, the samples were incubated with primary antibodies used were against FOXD1 (dilution 1:1000, abs102294; Absin, Shanghai, China), FOXD2 (dilution 1:1000, TA322739; Origene, Rockville, MD, US), FOXD3 (dilution 1:1000, abs133773; Absin, Shanghai, China), and FOXD4 (dilution 1:1000, orb156929; Biorbyt, Cambridge, UK) overnight at 4 °C. Following incubation, the membrane was washed thrice with Tris-buffered saline (TBS) 0.1% Tween 20 (TBST) for 5 min each. Next, the membranes were placed in a diluted secondary antibody solution and incubated at room temperature for 2 h. After incubation with the secondary antibody, the membrane was washed with TBST three times for 10 min each and then washed once with TBS for 10 min. All the bands were detected using the ECL detection kit (SW 134; Seven), according to the manufacturer’s protocol. The results were analyzed using the ImageJ software. All experiments were performed in triplicate.

RNA isolation and quantitative real-time polymerase chain reaction (qRT-PCR)

We used TransZol Up (ET111-01-V2; TransGen Biotech, Beijing, China) to extract total RNA from the cell lines (normal colon epithelial cells and four colorectal cancer cell lines), five freshly resected human colon cancer tissues, and seven adjacent normal colon tissues. The A260/A280 ratio of the extracted RNA ranged from 1.8–2.0. We used Cytation three image reader to detect the total RNA concentration and added an appropriate amount of RNase-free water to obtain a concentration of 1 ug/ul. The experiment was performed according to the manufacturer’s instructions (Chen et al., 2022). The primer sequences used are listed in Table S6. cDNA was synthesized from the extracted RNA samples (1 µl) using a PerfectStart ® Uni RT&qPCR kit (AUQ-01-V2; TransGen Biotech, Beijing, China), according to the manufacturer’s instructions and qPCR was performed using a Applied Biosystems StepOne Real-Time PCR kit (AQ311-01; TransGen Biotech, Beijing, China). We used qPCR to detect differential expression of the FOXD subfamily in cells and tissues. The relative gene expression of the FOXD subfamily was calculated using the 2^−ΔΔCT method. Expression was normalized to that of GADPH. The qPCR was performed in triplicate.

Statistical analysis

All data analyses were performed using R (version 4.0.4, http://r-project.org/). To ensure reliability of results, t-test, chi-square test, Kaplan–Meier survival analysis, and Cox proportional risk model were performed. The t-test was used to compare continuous variable (such as gene expression values) differences between the two groups, and the chi-square test was used to analyze the associations between categorical variables, such as the relationship between gene mutation status and clinical characteristics. Kaplan-Meier survival analysis was used to assess the difference in survival between different groups (e.g., high and low gene expression groups), and significance was determined by log-rank test. The Cox proportional risk model was used in multivariate survival analysis to assess the combined impact of multiple variables (gene expression level and clinical features) on survival outcomes. Statistical analyses of qPCR and western blotting were performed using GraphPad Prism software (version 5.0). The significance threshold for all statistical analyses was set at p < 0.05, unless otherwise specified.

Characteristics of the FOXD subfamily

We determined the diagnostic potential, and mutational characteristics of FOXD1, FOXD2, FOXD3, and FOXD4 in CRC. The area under the curve (AUC) of the receiver operating characteristic (ROC) was calculated. The AUC values for FOXD1, FOXD2, FOXD3, and FOXD4 were 0.732, 0.784, 0.949, and 0.750, respectively (Figs. 1A–1D). With an AUC value of 0.949, FOXD3 was the most significant diagnostic biomarker for CRC. Compared to known CRC markers such as CEA and CA19-9, FOXD3 demonstrated greater diagnostic accuracy.

Mutation analysis of the FOXD subfamily were visualized in lollipop plots (Figs. 1E–1H). FOXD2 was mutated in both COAD and READ, FOXD3 and FOXD4 were mutated only in COAD, and FOXD1 mutations were not detected. Green dots represent missense mutations and blue dots represent gene deletions. These results provide biological clues regarding the characteristics and functions of FOXD in CRC.

Identification and visualization of FOXD subfamily-related DEGs

We divided the FOXD subfamilies into two groups based on the median DEG expression values. These two groups were the high- and low-expression groups. We performed differential gene expression analysis to identify the FOXD subfamily related DEGs. The results are represented in volcano plots (Figs. 2A–2D). We identified 590, 1,388, 737, and 1,122 FOXD1-4-related DEGs, respectively. Up- and down-regulated genes are represented in red and blue, respectively. The top 50 up- and down-regulated genes are depicted in heat maps (Fig. 2E–2H). Blue represents low FOXD expression and red represents high FOXD expression.

Functional enrichment analysis of FOXD subfamily-related DEGs

Next, we performed GO/KEGG analysis of the FOXD subfamily related DEGs (Fig. 3). GO analysis (Table S3) revealed correlations between FOXD1-related DEGs and metabolic development. FOXD2-related DEGs were associated with humoral antimicrobial responses and cornification. FOXD3-related DEGs were associated with axonogenesis and the modulation of chemical synaptic transmission. FOXD4-related DEGs were associated with the formation of quadruple SL/U4/U5/U6 snRNPs and mRNA transsplicing. KEGG pathway analysis (Table S4) revealed that FOXD2-related DEGs were linked to cytokine receptor pathways. Moreover, FOXD3-related DEGs were associated with the calcium signaling pathway. FOXD4-related DEGs correlated with salivary secretion and GABAergic synaptic pathways. No pathway enrichment was observed for FOXD1-related DEGs.

Gene set enrichment analysis

We subsequently performed GSEA of the FOXD subfamily related genes (Fig. 4, Table S5). GSEA showed that the FOXD subfamily-related genes were mainly linked to KEGG pathways, including cytokine-cytokine receptor interactions and ECM-receptor interactions.

Protein-protein interaction analysis and Hub genes identification

The FOXD subfamily-related DEGs were uploaded to the STRING database for PPI network analysis. The results were then imported into Cytoscape. DEGs were shown in different colors based on logFC expression values (Figs. 5A–5D, top panel). Red and blue represented the up- and down-regulated genes, respectively. The deeper the color, the larger the |logFC| value. Among the FOXD subfamily, FOXD2 and FOXD3 had more hub genes. Fig. 5A showed that the hub genes of FOXD1 were KRT5, KRT6A, KRT6B, KRT6C, KRT7, KRT14, KRT16, and KRT24. There were 21 hub genes of FOXD2, all of which were down-regulated (Fig. 5B). FOXD3 had 13 hub genes, all of which were also down-regulated, such as CXCL5, NYP, ADCY2, and so on (Fig. 5C). SPRR1A, SPRR1B, SPRR2E, SPRR2F, and SPRR3 were the hub genes of FOXD4 (Fig. 5D).

Analysis of immune cell infiltration

We further estimated the changes in immune cell infiltration in CRC using the CIBERSORT algorithm (Fig. 6A). The relationship between the levels of immune cell infiltration is shown in red and blue, with blue indicating a positive correlation and red indicating a negative correlation (Fig. 6B). The correlations between immune genes and FOXD1, FOXD2, FOXD3 and FOXD4 are shown as lollipop plots. The deeper the color, the smaller the p-values, and the stronger the correlations (Fig. 6C–6F).

Survival analysis and clinical correlation analysis of FOXD subfamily genes

We performed a prognostic survival analysis according to the high- and low-expression groups of FOXD subfamily genes (Fig. 7A–7D). The relationship between the levels of immune cell infiltration is shown, and positive and negative correlations are represented as blue and red, respectively. The clinical correlation parameters included tumor T-stage, N-stage, M-stage, sex, age, pathologic stage, history of colonic polyps, lymph node invasion, and FOXD subfamily genes and are displayed in a nomogram (Figs. 7E–7H). Using these nomograms, we predicted 1-year, 3-year, and 5-year patient survival prognoses. Finally, the clinical prediction models were assessed using calibration curves (Figs. 7I–7L). The closer the curves of 1-year, 3-year, and 5-year to the grey ideal line, the better the predictions.

Analysis of FOXD subfamily expression by IHC staining of patient samples

In human colon cancer tissues, protein expression levels of the FOXD subfamily were determined using IHC. FOXD1 expression was not detected in both tumor and matched normal tissues. The moderate staining was observed in FOXD2 and FOXD4. FOXD3 was strong staining. Overall, FOXD2, FOXD3 and FOXD4 were expressed at higher levels in tumors, compared with the adjacent normal tissue (Fig. 8).

Identification of FOXD subfamily gene expression in vitro

Western blotting results showed that FOXD1 expression was higher in all CRC cell lines than normal colon cells (CCD841CON). Among them, the difference between SW620 cell line and CCD841CON cell line was the most significant (Fig. 9A). Compared to the CCD841CON cell line, FOXD2 was expressed at a higher level in SW480, SW620, HT29, and HCT116 cell lines. However, the difference was statistically significant only in the SW480, SW620, and HCT116 cell lines (Fig. 9B). FOXD3 expression seemed to be higher in all CRC cell lines, but the difference was statistically significant only in the HT29 cell line (Fig. 9C). In addition, the FOXD4 expression level seemed to be higher in all CRC cell lines, but the difference was statistically significant only in the SW620 cell line (Fig. 9D). Next, we examined the expression of FOXD1 and FOXD4 in tissues and cell lines using qPCR (Fig. 10). The results showed that the expression levels of FOXD1 and FOXD4 were higher in tumors compared with the paracancerous tissues (Figs. 10A, 10C). However, FOXD1 expression seemed to be higher in SW480, HT29, and HCT116 cell lines than in normal colon cells (CCD841CON) but the difference was not statistically significant (Fig. 10B). The expression level of FOXD4 was higher in all CRC cell lines than in normal colon cells, but the difference was statistically significant only in the SW620 and HT29 cell lines (Fig. 10D). The list of primer sequences was shown in Table S6.