Pan-cancer analysis of homeodomain-containing gene C10 and its carcinogenesis in lung adenocarcinoma

Patients and clinical specimens

The mRNA expression profile and clinical information of 33 types of cancers and normal tissues (eQTL: expression quantitative trait loci, https://www.gtexportal.org/home/) were extracted from The Cancer Genome Atlas data platform. We also downloaded the immune cells and immunophenotype scores from the ImmuCellAI. The clinicopathological features of the 46 patients with benign lung disease were summarized in Supplementary Table 1, which were obtained during surgery and the clinicopathological features of 408 LUAD specimens are summarized in Supplementary Table 2, which were obtained during surgery or needle biopsy. All tissues were obtained from January 2010 and December 2020 at Jiangmen Central Hospital (Guangdong, China). The diagnoses were based on clinical and pathological results. The use of these clinical materials for research purposes have obtained prior consent of patient and approval by the Institutional Research Ethics Committee of Jiangmen Central Hospital. Two independent professional pathologists assessed the ratio of tumor vs non-tumor of hematoxylin and eosin (H&E) stained tissue samples. More than 75% of the clinical lung cancer tissue samples analyzed in this study were used for further analysis.

Bioinformatic analysis

In this study, the following bioinformatic analysis methods were used: (1) Differentially expressed genes analysis was performed using R “limma” packages; (2) Using the follow-up data and HOXC10 expression, we analyzed the correlations of HOXC10 with four clinical outcomes: overall survival (OS: the period from date of diagnosis until death from any case), disease specific survival (DSS: the percent of people who died from a specific disease in a defined period of time), disease free interval (DFI: the measure of time after treatment during which no sign of cancer is found), and progression free interval (PFI: the length of time during and after the treatment of a disease) [15]. The Kaplan-Meier survival curves were plotted, and the forest plot of the results of univariate Cox regression with hazard ratio (HR) and their 90% confidence interval (CIs); (3) Gene set variation analysis was performed to compare the differentially expressed pathways of different levels of HOXC10 expression using “GSVA” package; (4) we explored the correlations of HOXC10 with tumor mutation burden (TMB) and microsatellite instability (MSI) using Pearson correlation analysis. (4) We explored the tumor microenvironment by calculating ESTIMATE score, stromal score, and immune score. (5) Pearson correlation analysis was used to present the correlations of HOXC10 expression with 22 immune score, 22 kinds of immune cells and function, immune activated genes, immune inhibited genes, chemokines genes and chemokines receptors genes.

Immunohistochemistry (IHC)

IHC analysis was conducted to examine HOXC10 expression in 46 benign lung disease and 408 LUAD samples, VEGFA and CD31 expression in 238 LUAD samples with pathological grade II. The IHC performance and expression scoring criterion were stated previously [16–18]. The formalin-fixed paraffin-embedded sections were boiled in TE (pH 9.0) buffer for 10 min to extract antigens, blocked with hydrogen peroxide and goat serum, then incubated overnight at 4° C in a humidified box with anti-HOXC10 antibody (1:100, Cat#PA5-31078, Invitrogen, Carlsbad, CA, USA), anti-VEGFA antibody (1:200, Cat#19003-1-AP, Proteintech, Wuhan, China), and anti-CD31 antibody (1:2000, Cat#11265-1-AP, Proteintech) in Antibody Diluent (Abcam, Cambridge, MA, USA). After incubation, slides were rinsed in TBS/0.05% Tween 20, incubated with biotin-conjugated secondary antibody (Proteintech) and horseradish peroxidase (HRP)-conjugated streptavidin (Proteintech) 30 min at room temperature, respectively, stained by 3,3′-Diaminobenzidine (DAB) Liquid Substrate Dropper System (Sigma-Aldrich, St Louis, MO, USA). Staining index (SI) evaluated by two experienced independent investigators was averaged for comparative evaluation for HOXC10 and VEGFA expression. SI was the product of the ratio of positive tumor cells and staining intensity score. Positive tumor cell ratio was scored as follows: 0 (no positive tumor cells); 1 (< 10% positive tumor cells); 2 (10-35% positive tumor cells); 3 (35-70% positive tumor cells) and 4 (> 70% positive tumor cells). Staining intensity score was according to the following criteria: 0 (no staining); 1 (weak staining, light yellow); 2 (moderate staining, yellow brown) and 3 (strong staining, brown). Based on this method, expressions of HOXC10 and VEGFA in LUAD samples were evaluated with SI scores of 0, 1, 2, 3, 4, 6, 8, 9 or 12. SI score 4 was the median of tumor tissues. Low and high expression of HOXC10 and VEGFA were stratified by the following standard: SI ≤ 4 was conducted tumor samples with low level of HOXC10 and VEGFA, and SI score of > 4 as tumors with high level of HOXC10 and VEGFA. M8 digital microscope (PreciPoint, Freising, Bavaria, Germany) was used to capture images of human lung tissue under 10× and 40× objective magnification.

Microvascular density (MVD)

CD31 stained slides were examined to determine MVD in clinical samples. For each slide, 3 hot spot fields (the largest number of stained vessels) were evaluated under a 50× light microscope. The pathologist subsequently counts the blood vessels at each "hot spot" under a 200 × (high power) light microscope. The calculation of MVD of each hot spot was based on the total number of vessels divided by the area. The MVD of the specimen was used as the average of the MVD of the three histological regions. For statistical analysis, the MVD was divided into low (below the median) and high (greater than the median) groups [19].

Subtraction enrichment (SE) of blood samples for CTCs

For blood samples, using Human Circulating Tumor Cell Subtraction Enrichment Kit (Cat#SHE-011, Cytelligen, SanDiego, CA, USA) enriched CTCs according to the manufactures’ direction, collecting 6 ml peripheral blood, and 200 ×g centrifugal for 15 minutes. 3.5ml hCTC buffer was used for softly mixing with the precipitated blood cells and subsequently loaded on nonhematopoietic cell separation matrix, then centrifuging at 450 × g for 5 min. The supernatant containing tumor cells was resuspended to 14 mL hCTC buffer. Samples were then rotated at room temperature (500 × g) for 4 min, then the supernatant is aspirated to 50 μl. Precipitated cells were softly resuspended and bleed with the cell fixative (Cytelligen, San Diego, CA, USA). Samples were painted on CTC slides, and dried overnight for subsequent immunostaining-fluorescence in situ hybridization analyses.

Identification of CTCs using immunostaining-fluorescence in situ hybridization (iFISH)

iFISH was conducted according to the manufacturer’s guidance (Cytelligen). Dry monolayer-coated slides were incubated with PBS at room temperature for 3 min, followed by hybridized Vysis chromosome 8 centromere probe (CEP8) SpectrumOrange (Abbott Laboratories, Chicago, IL, USA) for 4 hours utilizing an S500 StatSpin ThermoBrite Slide Hybridization/Denaturation System (Abbott Molecular, Abbott Park, IL, USA). Subsequently, samples were hatched with Alexa Fluor (AF) 488 (green)-conjugated monoclonal anti-CD45 antibodies (1:200, Cat#FSH-002, Cytelligen, San Diego, CA, USA) for 30 min in the dark. Cell nucleus staining were performed with DAPI (Beyotime, Shanghai, China). CTC image scanning and analysis were used to obtain the recognized tumor cell image automatically. CTCs were defined as DAPI+, CD45-, CEP8 >2 signals.

Statistical analysis

The chi-squared test was used to analyze the connection between HOXC10 expression and clinicopathological characteristics. For normal distributed variables, data are presented as the mean ± standard deviation (SD) and the statistical difference between the two groups was determined using the student’s t-test. The statistical differences between multiple groups were determined by one-way ANOVA. The median and interquartile range present continuous variables without normal distribution. Nonparametric test (Mann-Whitney U-test) was used to assess the significance of the differences. Survival curves for both HOXC10-high and HOXC10-low patients (median) were utilized by the Kaplan-Meier method, which were compared by the log-rank test. P < 0.05 indicated a statistically significant difference. All analyses were conducted by SPSS 23.0 software (IBM Corp, Chicago, IL, USA) and diagrams were produced by GraphPad Prism 8.0 software (GraphPad Inc, San Diego, CA, USA).

Differential expression of HOXC10 and its associations with prognosis in pan-cancer

We first evaluated the expression levels of HOXC10 between normal and tumor pan-cancer tissues. The result indicated that HOXC10 was significantly up-regulated in most cancers including BRCA, CESC, DLBC, ESCA, GBM, KIRC, KIRP, LUAD, LUSC, PAAD, SARC, STAD, and THYM (Figure 1A). However, the expression of HOXC10 was remarkably down-regulated in KICH, LAML, OV, PRAD, READ, SKCM, TGCT, THCA, UCEC (Figure 1A). Next, we analyzed the associations between HOXC10 expression and four clinical outcomes. Our results shown that elevated HOXC10 was associated with poor OS in LGG, ACC, COAD, GBM, OV, BRCA, LIHC, LUSC, THCA, MESO, and PAAD but was negatively associated with OS in KIRC, KIRP, SKCM, and STAD (Figure 1B). The Kaplan-Meier survival curve also shown same trends, and indicated high-expressed HOXC10 was associated with poor OS in PRAD and UCS (Supplementary Figure 1). Similarly, our results suggested that the high expression of HOXC10 was positively associated with poor DSS in LGG, COAD, ASS, LIHC, GBM, OV, THYM, STAD, MESO, CESC, UCS, and BLCA. But favorable DSS was found for patients with elevated HOXC10 expression in KIRC, KIRP and SKCM (Figure 1C). The Kaplan-Meier survival curve identified unfavorable and favorable roles in these cancers (Supplementary Figure 2). For DFI, high-expressed HOXC10 increased the risk of poor prognosis in LGG, LUAD, and ACC but decreased the risk of poor prognosis in KIRP, UCS, PCPG, and BLCA (Figure 1D). The Kaplan-Meier survival curve also suggested that elevated increased the risk of poor DSS in HNSC, MESO, and TGCT (Supplementary Figure 3). Results for PFI indicated that elevated-HOXC10 was associated with poor PFI in LGG, ACC, COAD, THYM, GBM, CESC, BLCA, STAD, OV, LUAD, MESO, LIHC, and HNSC but was favorable for DSS in KIRC, KIRP, SKCM, PCPG, and SARC (Figure 1E). Same trends were also found in Kaplan-Meier analyses (Supplementary Figure 4).

Pathways enrichment and molecular features of HOXC10 in pan-cancer

To explore the biological function roles of HOXC10 in cancers, we analyzed the correlations of the pathway’s enrichment differences with HOXC10 expression (Figure 2). We found TNFA signaling via NFKB, reactive oxygen species pathway, PI3K/AKT/MTOR signaling, P53 pathway, oxidative phosphorylation, MTORC1 signaling, glycolysis, fatty acid metabolism, hypoxia and apoptosis were totally positively enriched in HOXC10 for whole cancers. On the contrary, Notch signaling, Wnt-beta catenin signaling, myogenesis, and Hedgehog signaling were almost negatively enriched in HOXC10. Some pathways showed dual roles such as Kras signaling, interferon gram and alpha response, IL6/JAK/STAT3, IL2/STAT5 signaling were negatively enriched in HOXC10 for CESC, COAD, ESCA, HNSC but were positively enriched in HOXC10 for other cancers. Some immune-related pathways such as apical surface and junction, angiogenesis, and allograft rejection were also negatively enriched in CESC, ESCA, and HNSC.

We further explored the correlations of HOXC10 with TMB and MSI. We found that HOXC10 was positively associated with TMB in LIHC, LUAD and was negatively associated with TMB in GBM, KIRC, STAD, and BLCA (Figure 3A). While HOXC10 showed a positively connection with MSI in GBM, LGG, CESC, LUSC, LUAD, and OV but a negative correlation to MSI in SKCM, STAD, BRCA, LIHC, PAAD, TGCT, and UCEC (Figure 3B).

We estimated the correlations of HOXC10 with tumor microenvironment using ESTIMATE and stromal scores. For KIRP, LUAD, STAD, PCPG, CESC, PAAD, SARC, KIRC, ESCA, BRCA, LUSC, and BLCA, the expression of HOXC10 was negatively associated with ESTIMATE score while showed a positive connection only for COAD. For stromal score, negative associations between HOXC10 expression and stromal score were found in KIRP, PCPG, PAAD, STAD, LUAD, BRCA while HOXC10 presented positive connections with stromal score for TGCT, COAD, LGG, KIRC, and PRAD.

Association between HOXC10 and immune status in pan-cancer

To explore the immune regulation role of HOXC10 in cancers, we analyzed the associations between HOXC10 expression and immune score and immune-related genes. Our results showed that HOXC10 was negatively associated with immune score in GBM, LUAD, LUSC, TGCT, CESC, ESCA, KIRP, PAAD, KIRC, SARC, STAD, and BRCA (Figure 4A–4L). However, a positive connection was found between HOXC10 expression and immune score in COAD.

We next evaluated the correlations of HOXC10 with immune infiltration in cancers. Our results indicated that HOXC10 was almost positively correlated with all immune cells infiltration levels in KICH, LGG, COAD, PRAD, THCA, TGCT and THYM. On the contrary, these immune cells negatively infiltrated with HOXC10 in ACC, MESO, KIRP, CESC, STAD, GBM, LUSC, SARC, and SKCM. Only, neutrophil and CD56 dim natural killer cell were positively with HOXC10 in KIRC. In LUAD and HNSC, immune cells show similar trends (Figure 4N).

Immune-activated genes were positively associated with HOXC10 expression in THYM, LGG, KICH, LIHC, COAD, THCA, DLBC, CHOL, and PRAD. However, HOXC10 seems to inhibit these immune genes expression in ACC, KIRP, STAD, SKCM, CESC, SARC, MESO, GBM, LUAD, LUSC, KIRCT, TGCT, and UCS. For immune-inhibited genes, HOXC10 was positively associated with these genes such as CD274 in PCPG, LGG, COAD, PRAD but showed negative with CD274 in CESC, STAD, SKCM, GBM, KIRP, ESCA, and LUAD (Figure 5A, 5B).

The directed migration of immune cells is a necessary condition for the occurrence and completion of immune response. Chemokines are a class of cytokines that control cell directed migration. Its function is mediated by chemokine receptors. The interaction between chemokines and their receptors controls the directional migration of various immune cells in the circulatory system and tissues and organs, so that they can go to the site of abnormal proliferation, destroy abnormal proliferation cells, and maintain the function of tissue cell balance. Our results suggested that chemokines and their receptors were positively associated with HOXC10 expression in LGG, THYM, PRAAD, CHOL, LIHC, READ, KICH, DLBC, and TGCT but were remarkably negatively associated with HOXC10 expression in ACC, SKCM, CESC, STAD, KIRP (Figure 6A, 6B).

HOXC10 is upregulated in LUAD, especially in distant metastasis tissues, by bioinformatics analysis

Since a large amount of literature supports the standpoint that HOX family genes are abnormally expressed in several types of human cancer and play vital roles in cancer progression and distant metastasis, the significance of HOX genes in LUAD seized our great momentum. By analyzing the mRNA expression of 39 HOX genes in the LUAD dataset from TCGA, we found that 18 HOX genes were significantly higher in LUAD tissues with distant metastasis in 3 years (n=249) compared with that without distant metastasis in 3 years (n=68), and HOXC10 was identified as the most upregulated HOX gene (Figure 7A). Therefore, we focused on the HOXC10 for further study. Moreover, the mRNA level of HOXC10 was significantly upregulated in the LUAD tissues (n=517) from TCGA compared with that in adjacent normal tissue (n=57) (Figure 7B). Similarly, HOXC10 was upregulated in LUAD tissue from TCGA (n=57) compared with matched adjacent normal tissue (Figure 7C). Importantly, the significance of HOXC10 expression on DMFS was analyzed in the LUAD dataset from TCGA, and LUAD patients with high expression of HOXC10 exhibited shorter DMFS compared with LUAD patients with the low level of HOXC10 (Figure 7D). Therefore, these findings from the LUAD dataset of TCGA indicated that HOXC10 is upregulated in LUAD tissues, especially with distant metastases, which predicted worse DMFS.

Overexpression of HOXC10 in human LUAD samples indicates advanced progression and poor prognosis

To investigate the clinical significance of HOXC10 in LUAD, we further examined the expression of HOXC10 in a larger number of clinical samples, including 46 benign lung lesions and 408 LUAD tissues, by immunohistochemistry (IHC). As shown in Figure 7E–7G, HOXC10 was mainly expressed in the nucleus of cancer cells, and HOXC10 expression was dramatically higher in LUAD tissues than in benign lung lesions. HOXC10 upregulation was observed in 48.28% (197/408) LUAD tissues (Figure 7F and Table 1). Notably, the expression of HOXC10 was gradually increased with the clinical stage (Figure 7G).

It is well known that the grade and stage of advanced tumors are significantly associated with the accelerated development of LUAD. Our statistical data indicated that elevated HOXC10 expression was positively correlated with pathological grade, N classification, M classification and clinical stage (Table 1). Furthermore, to analyze the effect of tumor stage and grade in LUAD patients with overexpressed HOXC10 on poor DMFS, LUAD patients were first stratified into different groups of Stage I to IV and Grade I to III. As shown in Figure 8A–8C, based on Kaplan-Meier survival analysis, high expression of HOXC10 exhibited poorer DMFS of LUAD patients with all Stage, Stage I-II, and Stage III-IV. Consistently, LUAD patient with high level of HOXC10 predicted shorter DMFS in the same pathological grade compared with low levels of HOXC10 (Figure 8D–8F). Collectively, these results demonstrated that the HOXC10 overexpression was significantly positively correlated with advanced progression and poor prognosis in LUAD patients.

HOXC10 is an independent prognostic marker for DMFS of LUAD

To investigate whether HOXC10 is an independent prognostic value of DMFS in LUAD patients, univariate analysis of clinicopathological features and HOXC10 expression associated with DMFS in 408 LUAD patients was performed by Cox-regression analysis. LUAD patients with high levels of HOXC10 had a significantly higher risk of distant metastasis than patients with low levels of HOXC10 (hazard ratio [HR] = 2.426, 95% confidence interval [CI] = 1.417-4.153, P-value = 0.001). Similar results were shown for T classification (HR=1.887, 95%CI: 1.048-3.398, P=0.034), M classification (HR=31.258, 95%CI: 7.143-136.786, P<0.001), and clinical stage (HR=2.132, 95%CI: 1.264-3.598, P=0.005), while no significant correlation was found between N classification, age, gender, and pathological grade (Table 2). Multivariate analysis indicated that HOXC10 was an independent predictor for DMFS (HR=2.195, 955CI: 1.272-3.788, P=0.005, Table 2). Taken together, these studies suggest that HOXC10 is an independent prognostic marker for DMFS of LUAD.

Correlation of HOXC10 expression with the expression of angiogenic markers in LUAD

The growth of new capillaries in tumor tissue contributes to poor progression and metastasis, and it has been reported that HOXC10 is highly expressed in glioma and promotes angiogenesis [20]. To determine the possible association between HOXC10 and angiogenesis in LUAD patients, VEGFA and CD31 in LUAD patients were examined by IHC. Representative staining images immunohistochemical of HOXC10, VEGFA and CD31 in each group are shown in Figure 8G. HOXC10 expression was positively correlated with VEGFA expression (r = 0.372, P < 0.001) (Figure 8H, 8I). Moreover, MVD was efficiently enhanced in HOXC10-high tissues versus HOXC10-low tissues, as shown by CD31 staining (r = 0.387, P < 0.001) (Figure 8J, 8K). These data indicated that HOXC10 expression is positively correlated with angiogenic marker in LUAD.

Correlation of HOXC10 expression with the number of CTC clusters in LUAD patients

The invasion of circulating tumor cells (CTCs) into peripheral blood represents a critical process to tumor metastasis, progression, and angiogenesis. Therefore, the clinical correlation of HOXC10 expression with the number of CTC clusters in LUAD patients was further analyzed. Immunofluorescence staining of CD45 and DAPI and iFISH with centromeres of chromosome 8 (CEP8) were used to identify CTCs in patients with LUAD. CD45-/CEP8 >2 cells were defined as CTCs, while CD45+/CEP8 = 2 cells were defined as white blood cell (WBC) (Figure 9A). We found that the number of CTC clusters in LUAD patients with high levels of HOXC10 (n= 34) was more than that in those with low levels of HOXC10 (n= 36) (Figure 9B). While the number of CTC-WBC clusters between them was no significant difference (Figure 9C). Our results indicated that HOXC10 expression levels is positively correlated with the number of CTC clusters in LUAD patients.