Transcriptomic data from 638 colorectal cancer (CRC) tumour tissue (colon adenocarcinoma, COAD and rectum adenocarcinoma, READ) and 51 normal tissue samples of TCGA-CRC was acquired from the TCGA (http://cancergenome.nih.gov/) database. Normal tissue samples from 4,010 healthy donors were obtained from the GTEx database (https://gtexportal.org). All data were standardised using the FPKM method to eliminate the influence of gene length on expression levels and to correct for sequencing depth-induced biases across samples.

Differentially expressed genes (DEGs) were obtained based on TCGA-CRC and differentially expressed miRNAs (DEMs) were obtained based on TCGA-CRC miRNA by differential analysis (p.adj < 0.05, |log2FC|>1) between CRC and normal groups using the DESeq2 package (19).

To investigate the expression of TBX5 in several cancers, TBX5 expression was analysed in 33 tumour and normal tissues using the TCGA and Genotype Tissue Expression (GTEx) databases. Protein expression of TBX5 was analysed between CRC and normal tissues using the Human Protein Abstracts Database (HPA, https://www.proteinatlas.org/).

Based on TBX5 expression, CRC patients were categorised into high and low expression groups. To analyse the relationship between TBX5 and the prognosis of CRC patients, Kaplan-meier (K-M) survival analysis was executed. To analyse the relationship between TBX5-related genes and the prognosis of CRC patients, spearman analysis was executed between TBX5 and DEGs. Based on the expression of the TBX5-related genes, patients were classified into high and low expression groups using the survminer package to calculate the optimal cutoff value, and K-M survival analysis was executed.

TBX5 and clinical characteristics (age, gender, T, N, and M stages) were selected to build a nomogram for predicting CRC patient survival. Firstly, based on the samples with survival information in TCGA-CRC, univariate Cox regression analysis was performed using the coxph function from the “survival” package (v 3.5-3) (20) [with p<0.05, hazard ratio (HR)≠1]. Then, to ensure the validity and reliability of the univariate results, proportional hazards (PH) assumption test was conducted using the cox.zph function from the “survival” package (v 3.5-3) (with p>0.05). Subsequently, the “survival” package (v 3.5-3) was further used to perform multivariate Cox regression models on the selected variables (with p<0.05, HR≠1), and PH assumption test (with p>0.05) were conducted to identify the independent prognostic factors. The “forestplot” package (v 3.1.1) (21) was used to visually present the results of univariate and multivariate Cox analyses by creating forest plots. Then, to explore the predictive power of the independent prognostic factors on CRC patient survival rates, a nomogram for 1-, 3-, and 5-year survival rates of CRC patients was constructed using the “rms” package (v 6.5-1) (22). Finally, to verify the accuracy and reliability of the nomogram, calibration curves and ROC curves for 1-, 3-, and 5-year survival rates, as well as 1-, 3-, and 5-year overall survival periods, were respectively drawn using the “rms” package (v 6.5-1) and the “survivalROC” package (v 1.0.3) (23).

To validate the relationship of TBX5 with cell stemness, EMT pathway and angiogenesis related genes, single sample gene set enrichment analysis (ssGSEA) scores were calculated for each gene set using the ssGSEA algorithm of the GSVA package (24). Spearman analysis was executed between TBX5 gene expression and per scores.

GSEA was executed to explore the function of TBX5. Firstly, the c5.go.bp.v7.4.symbols.gmt, c5.go.cc.v7.4.symbols.gmt, c5.go.mf.v7.4.symbols.gmt and c2.cp.kegg.v7.4.symbols.gmt genes sets were obtained from the molecular signatures database (MsigDB) (http://www.broadinstitute.org/gsea/msigdb/) as a reference gene set. Then, Based on the TCGA-CRC dataset, differential analysis between all samples in two groups were executed and ranked from largest to smallest by multiples of variance. Finally, GO and KEGG enrichment analyses were executed using ClusterProfiler package (p.adjust<0.05) (25).

To study the immune cell infiltration in the two groups, the ssGSEA algorithm was used to calculate the infiltration abundance of 28 immune cells between the two groups (26). The Wilcoxon signed-rank test was employed to compare differences in infiltration levels between the high- and low-expression groups of TBX5 within the training set. Concurrently, the Benjamini-Hochberg (BH) method was applied to correct the p-values, with both corrected and uncorrected results reported. Furthermore, to study the relationship between differential immune cells and CRC prognosis, K-M survival analysis was performed and significance of differences analysed using the log-rank test.

To explore the function of TBX5, TBX5 was included in the GeneMANIA (http://www.genemania.org) database for analysis. The miRNAs associated with TBX5 were forecasted by targetScan databases (v 7.2, published March 2018). Taking intersections of miRNAs targeting TBX5 and DEMs to obtain key miRNAs. lncRNA associated with key miRNAs were predicted via STARBASE databases (v 2.0, published January 2014). Moreover, cytoscape software was applied to optimize the results of mRNA-miRNA-lncRNA network.

The 50% inhibitory concentration (IC50) of each of the 128 drugs in TCGA-CRC was calculated for each sample using pRRophetic package (27) (p.adjust<0.05). To investigate the regulatory relationship between Gegen Qinlian Decoction (GQD) and TBX5. The corresponding target protein for each active ingredient in GQD was predicted by TCMSP database. The five drug-like principles (Lipinski rule): relative molecular weight (MW) ≤ 500, hydrogen-bonded donor (Hdon) ≤ 5, hydrogen-bonded acceptor (Hacc) ≤ 10, lipid-water partition coefficient (AlogP) ≤ 5, the value of oral bioavailability (OB) ≥ 30% and vaule of drug-likeness (DL) ≥0.18 were used as the screening conditions. Target proteins were abstracted and converted to gene names using UniProtKB (http://www.uniprot.org). Spearman’s correlation analysis was employed to assess the relationship between TBX5 and drug target genes. For each drug, the single target gene exhibiting the strongest correlation with TBX5 was selected and subjected to molecular docking with TBX5. First, the three-dimensional structure of TBX5 was obtained from the Uniprot database (https://www.uniprot.org/), with the corresponding crystal structure PDB ID being 2X6U (https://doi.org/10.2210/pdb2X6U/pdb). The 3D structure of the GQD active ingredient was retrieved from the TCM Database@Taiwan (https://tcm.cmu.edu.tw/) and HERB 2.0 (http://herb.ac.cn). The CB-Dock online platform (http://cao.labshare.cn/cb-dock/) was employed to predict protein-compound binding affinity and binding sites. Subsequently, AutoDock Vina was utilised to perform precise molecular docking validation and visualisation of the predicted results (docking binding energy < -5 kcal/mol).

The frozen cells were quickly thawed and dissolved. Subsequently, they were transferred to a centrifuge tube containing complete medium and centrifuged at 1000r/min to eliminate the supernatant. The cell suspension was then supplemented with DMEM complete medium, transferred to T-25 culture flasks, and incubated at 37°C with 5% CO₂. Upon reaching 80% confluence, cell passaging was carried out. The cells underwent PBS rinsing, digestion with 0.25% trypsin for single-cell suspension generation, followed by centrifugation. Finally, the cells were combined with a cryopreservation solution and placed into cryovials. After freezing at -80°C, the cells were transferred to liquid nitrogen for prolonged storage.

The frozen Human embryonic kidney (293T) cells were immediately rewarmed at 37°C, put in complete medium containing 10 ml, and centrifuged at 1000 r/min for 5 min at room temperature. After removing the supernatant, the cells were put to DMEM complete media and transferred to T-25 culture flasks, where they were cultivated at 37°C with 5% CO₂. After putting 293T cells in 4-5×10⁶/10cm dishes, incubate them for 24 hours at 37°C with 5% CO2. Plasmid No. 1 and the transfection reagent dilution No. 2 should be made in the following order Supplementary Table S1. After thoroughly combining No. 2 and No. 1 reagents, incubate the mixture at room temperature for 15 min. Drop by drop, 1 ml of the transfection mixture was added to the 293T cells when the cell density reached 80%-90%. Followed by, replace 10 ml of fresh 293T medium and virus medium 4–6 hours and 24 hours after transfection, respectively. Finally, cell culture supernatants were collected at 48, 72 hours post-transfection and centrifuged at 500 g for 10 min, after which the supernatants were collected for lentivirus infection.

The cells were divided into sh-TBX5 group and sh-NC group for infection. First, wall-adherent cells were evenly distributed into 6-well plates at a density of 3 × 10^5 cells per well for 18–24 hours prior to lentiviral infection to ensure that there were approximately 6 × 10^5 cells in each well. After waiting for the cell-adherent fusion to reach 70%, the original medium was replaced with 2 ml of fresh medium containing 8 μg/ml polybrene and the appropriate amount of virus suspension was added. The incubation was continued at 37°C for 8 hours, and then the virus-containing medium was replaced with fresh medium. If the lentivirus contains fluorescent proteins, fluorescent expression can generally be observed 48 hours after transfection, and the fluorescent expression will be more obvious after 72 hours. If the results are not satisfactory, the infection operation can be repeated and continuous screening can be carried out by adding 2ug/ml puromycin to the medium until one month’s time.

The total RNA was isolated with the help of TRIzol method. The SwsScript All-in-One First-strand-cDNA-synthesis SuperMIx for qPCR kit was used to synthesize cDNA. The SYBR Green qPCR Mix were used to perform the RT-qPCR reaction. Ct values were acquired, and the comparative gene expression was determined using the 2-ΔΔCt formula, with GAPDH serving as the internal control. Finally, GraphPad Prism 6 was used for data visualization and analysis. A list of primers was shown in Supplementary Table S2.

This study was approved by Institutional Animal Care and Use Committee of Zhejiang Center of Laboratory Animals and complied with the guidelines for animal care and use established by the Chinese government.15 female Balb/c nude mice were housed in individually ventilated cages, and the supplier ensured the high quality of the mouse strains through strict quality control measures. The study consisted of 3 groups of 5 mice each: the CRC group, the GQD group, and the TBX5 lentivirus transfection group.BALB/C-nu nude mice were inoculated with tumor cells in the right axilla.The nude mice were disinfected with an alcohol cotton ball wipe, fixed, and inoculated with 120 µL of cells in the right axillary subcutaneous area.The cotton swabs were pressed in the axilla after the completion of the inoculation to ensure that there was no leakage of cell suspensions.The growth status of the mice was observed on a daily basis, and the tumor was recorded once every 2 days. Tumor volume was recorded every 2 days.

Tumor tissue is taken and fixed in 4% paraformaldehyde solution, then dehydrated, embedded and sectioned. Sections are first antigenically repaired and blocked, then the TBX5 antibody is diluted 1:100 according to the antibody manufacturer’s instructions and the sections are incubated overnight in a 4°C refrigerator. The next day, sections are rewarmed in a 37°C warmer for 30 minutes, then goat anti-mouse/rabbit IgG secondary antibody is added and incubated at 37°C for 20 minutes. The final steps include DAB color development, hematoxylin re-staining, dehydration, transparency and sealing.

Statistical analysis was executed by R software. Differences between the two groups were executed by Wilcoxon test. The clusterProfiler was used to identified enrichment pathway.

A total of 3,681 DEGs (1,834 up regulated genes and 1,847 down regulated genes) were obtained by difference analysis (Figures 1A, B; Supplementary Table S3). Pan-cancer analysis showed that TBX5 expression was down-regulated in LUAD, LUSC, BRCA, ESCA, GBM, UCS, GBMLGG and TGCT, and up-regulated in other cancers, especially CRC-COAD (Figures 1C, D). Immunohistochemistry results showed that TBX5 was expressed at a higher level in CRC-COAD tissue compared to normal tissue. Moreover, in normal tissue, TBX5 was mostly expressed in the cytoplasm, whereas in CRC-COAD tissue it was mostly expressed in the nucleus, suggesting that the high expression of TBX5 into the nucleus played a regulatory role in CRC-COAD (Figures 1E, F).

Patients were categorised into high and low expression groups based on TBX5 expression. K-M survival analysis showed that patients in the TBX5 high-expression group had a lower survival rate (Figure 2A). The absolute value of correlation coefficient top 10 genes including CACNA1A, SPTBN4, RP11-35N6.1, CHD5, GREB1, SLC4A8, KLHL32, MAP2, CACNA2D3, and ATP8A2 were obtained by spearman analysis. K-M results showed that PTBN4, RP11-35N6.1, CHD5, SLC4A8, KLHL32, MAP2, and ATP8A2 were obviously different between the two groups (Figure 2B).

The univariate Cox regression analysis forest plot (Figure 3A) and the PH assumption test (p>0.05) (Supplementary Table S4-1) showed that age, T, N, and M stages were significantly related to the survival of CRC patients (p<0.01). Furthermore, the multivariate Cox regression analysis forest plot (Figure 3B) and the PH assumption test (p>0.05) (Supplementary Table S4-2) indicated that age (HR = 2.311; 95%CI = 1.424-3.75), N (HR = 1.954; 95%CI = 1.192-3.204), and M stages (HR = 3.189; 95%CI = 1.995-5.097) remained significant independent prognostic factors (p<0.05). Additionally, a nomogram based on these independent prognostic factors was constructed as shown in Figure 3C. Moreover, the calibration curves demonstrated a high degree of overlap between the predicted probabilities of 1-, 3-, and 5-year survival from the nomogram and the reference line (Figure 3D). Also, the ROC curves of the nomogram model revealed that the AUC values at various time points were all above 0.7, indicating that the nomogram model had high a predictive value (Figure 3E). In summary, these findings are helpful for enhancing the predictive efficiency and therapeutic effect of CRC and improving patients’ prognosis.

Based on TBX5, GO and KEGG enrichment were executed. The results of GO displayed that TBX5 was enriched in biological processes (BP) such as cardioblast proliferation etc, cellular components (CC) such as anchored component of synaptic membrane etc, and molecular functions (MF) such as inhibitory extracellular ligand gated ion channel activity etc (Figure 4A; Supplementary Table S5). And TBX5 were enriched in ecm receptor interaction, etc KEGG pathway (Figure 4B; Supplementary Table S6). Correlation analysis showed significant correlation between TBX5 and cell stemness, EMT and angiogenesis but the correlation coefficient was low. The scores of cell stemness, EMT and angiogenesis-related genes were all obviously upregulated in the TBX5 high-expression group (Figure 4C).

18 differently expressed immune cells (CD4 memory central T cells, monocyte, macrophage and mast cell etc) were obtained between the two groups. After BH correction, differences in these 16 immune cell types retained statistical significance. (Figure 5A). K-M results showed CD4 central memory T cells, eosinophil, B immature cells, mast cells, MDSC, monocyte etc were significantly associated with the prognosis of CRC patients (p value<0.05) (Figure 5B).

A total of 21 genes functionally similar to TBX5 were obtained, including NPPA, TBX15 and NKX2–5 etc, and were mainly associated with the pathways of cell fate commitment, embryonic morpheogenesis, and heart development (Figure 6A). A total of 475 DEMs were obtained by difference analysis (Figures 6B, C; Supplementary Table S7). A total of 13 key miRNA were obtained by overlapping miRNAs targeting TBX5 and DEMs (Figure 6D).1 mRNA, 3 miRNA and 19 lncRNA were acquired and a network of mRNA-miRNA-lncRNA was constructed, with some relationships like TBX5-hsa-mir-1270-TMPO-AS1 etc in the network (Figure 6E). Finally, pRRophetic predictions indicate that high TBX5 expression correlates with lower predicted IC₅₀ values for drugs such as AG.014699 (p.adjust<0.05), suggesting tumours with elevated TBX5 expression may exhibit greater sensitivity to these therapeutic agents (Figure 6F). However, these computational predictions require further validation through subsequent experimental studies.

A total of four drugs were screened, including Pueraria lobata, Coptis chinensis, Glycyrrhiza uralensis, and Scutellaria baicalensis. Pueraria lobata obtained 14 active ingredients corresponding to 220 protein targets, converted into 176 corresponding genes. According to the results of spearman analysis in Pueraria lobata, the correlation between MAP2 and TBX5 was the highest and significant (cor = 0.288, pvalue = 1.214098e-13). After screening, formononetin was found to have a regulatory relationship with MAP2 (docking binding energy = -6.1 kcal/mol) (Figure 7A). TBX5 docked better with formononetin (docking binding energy = -6.5 kcal/mol) (Figure 7B). Coptis chinensis obtained 35 active ingredients corresponding to 227 protein targets, converted into 177 corresponding genes. The highest and significant correlation between SLC6A3 and TBX5 was found in Rhizoma Coptidis (cor = 0.271, pvalue = 3.284383e-12). 3,4-dimethoxy-(8CI), hydroxytyrosol, (R)-Canadine, Corydaldine and FER had regulatory relationship with SLC6A3 and TBX5 (docking binding energy < -5 kcal/mol) (Figures 7C, D). Glycyrrhiza uralensis obtained 165 active ingredients corresponding to 316 protein targets, converted into 246 corresponding genes, SLC6A3 showed the highest and significant correlation with TBX5 (cor = 0.271, pvalue = 3.284383e-12). And HMO, 7-Methoxy-2-methyl isoflavone, Karenzu etc had regulatory relationship with SLC6A3 and TBX5 (docking binding energy < -5 kcal/mol) (Figures 7E, F). Scutellaria baicalensis obtained 105 active ingredients corresponding to 233 protein targets, converted into 171 corresponding genes, MAP2 showed the highest and significant correlation with TBX5 (cor = 0.28, pvalue = 1.214098e-13). And beta-sitosterol had regulatory relationship with MAP2 and TBX5 (docking binding energy < -5 kcal/mol) (Figures 7G, H). Mesh parameters for the molecular docking process are detailed in the Supplementary Table S8. In silico molecular docking predicted that several active components from GQD, such as formononetin, might interact with TBX5. This provides a hypothetical molecular basis for the observed effects and warrants further investigation to validate these interactions.

HCT116 cells were cultured and observed under the microscope. The cells were well attached to the culture dish, indicating that the cell growth status was good and suitable for further experiments (Figure 8A). sh-TBX5 lentivirus was successfully packaged in 293T cells with 100% green fluorescence intensity (Figure 8B). In order to determine the most effective knockdown vector, we performed RT-qPCR to evaluate the gene silencing efficiency. Compared with the sh-NC group, the expression level of TBX5 gene was significantly reduced in the sh-TBX5 group (Figure 8C). The successful establishment of a stable TBX5 knockdown cell line was confirmed and can be used for further analysis. To screen for stable cell lines, further RT-qPCR validation was performed. The results showed that compared with the sh-NC group, TBX5 gene expression was significantly decreased in the sh-TBX5-1, sh-TBX5–2 and sh-TBX5–3 groups, with the most obvious decrease in the sh-TBX5–2 group (Figure 8D). Therefore, the sh-TBX5–2 vector was selected for subsequent experiments.

The results of animal experiments showed that compared with the control group, the body weight of GQD group and TBX5 lentivirus transfected group increased and changed with time, and the tumors of the model group increased gradually, and the tumors of the GQD group and TBX5 lentivirus transfected group grew slowly (Figure 9A-D). Immunohistochemistry results showed that the GQD group showed low expression compared with the CRC group, and the sh-TBX5 group showed low expression compared with the CRC group. It indicates that GQD could inhibit tumor growth and TBX5 affected the development of CRC (Figure 9E).