Carboxypeptidase E is a prognostic biomarker co-expressed with osteoblastic genes in osteosarcoma

Patient samples

This study was conducted on osteosarcoma tissues and adjacent normal tissues collected from the Department of Orthopaedic Oncology Surgery of the Beijing JiShuiTan Hospital. All of the patients had received four cycles of standard preoperative adjuvant chemotherapy (high-dose methotrexate, ifosfamide, cisplatin, and doxorubicin, MAP + IFO) prior to tissue collection (Suehara et al., 2019). The tissue samples were stored in liquid nitrogen following surgery (Table 1). Sample was collected with written informed consent from patient’s family members, according to the ethics committee of Beijing Jishuitan Hospital approved protocols (Grant No. K2022-117-00).

Filtering of DEG from microarray

Agilent Human lncRNA Microarray 2019 (4*180k, Design ID: 086188) was used. Data analysis of 12 samples was completed by OE Biotechnology Co., Ltd., (Shanghai, China). Quantification of total RNA was carried out using a spectrophotometer, NanoDrop ND-2000 (Thermo Fisher Scientific, Waltham, MA, USA), while the integrity analysis of RNA samples was performed using an RNA analyzer, Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). The subsequent steps, including sample labeling, microarray hybridization, and washing, were executed following established protocols recommended by the manufacturer. In brief, the total RNA was subjected to reverse transcription to synthesize double-stranded cDNA, which was then converted into complementary RNA (cRNA) and labeled using Cyanine-3-CTP. The labeled cRNAs were subsequently hybridized onto a microarray, and following a thorough washing step, the microarrays were scanned using the Agilent Scanner G2505C (Agilent Technologies, Santa Clara, CA, USA).

The whole-gene expression datasets were subjected to principal component analysis (PCA) for dimensionality reduction and quality control. One subject was excluded due to poor clustering. We conducted differential gene analysis using the R software “limma” package (https://bioconductor.org/packages/limma/) according to the filter criteria (log2|fold change| > 2, FDR < 0.05). Subsequently, the gene expression profile and survival (survival state and survival time) of patients with OS were analyzed by Cox regression, and prognostic related genes were detected in the TARGET database. Screening high risk factors related to survival based on P < 0.05.

Protein-protein interaction network and hub gene

STRING (https://string-db.org/) was employed to predict interactions between functional proteins and to establish a PPI network. The network was conducted with Cytoscape software (version 3.9.1). The PPIs between the putative protein biomarkers were analyzed using the STRING web service database (https://string-db.org/). Hub gene analysis was conducted after considering the high degree of connectivity in the PPI networks analyzed by the Maximal Clique Centrality (MCC) plugin of Cytoscape (Chin et al., 2014). Here, the top 10 MCC values were considered hub genes and were presented in the form of a PPI network.

Evaluation of the expression level of CPE

Two published OS datasets (GSE126209, GSE33382) were obtained from the GEO database for the expression-level evaluation of CPE. We selected 11 paired samples from the GSE126209 series and 68 OS samples, including three pathological types of OS (“Fibroblastic” 8, “Osteoblastic” 51, “Chondroblastic” 9), from GSE33382 series.

Total RNA extraction and qRT−PCR

Twenty paired OS samples (Table 1) were subjected to quantitative reverse transcription-PCR (qRT-PCR) to determine the expression of CPE (NM_001873.4), TRSP1 (NM_001330599.2), PFKFB3 (XM_047425347.1), and FADS1 (XM_047426935.1). In accordance with the manufacturer’s instructions, TRIzol reagent (Invitrogen, Waltham, MA, USA) was used to extract total RNA (Invitrogen, Waltham, MA, USA). RT-PCR analysis was performed using the LightCycler® 480 Probes Master Mix (LightCycler® 480II; Roche, Basel, Switzerland). RT-qPCR was performed using a 2720 thermal cycler (Applied Biosystems, Waltham, MA, USA). The relative mRNA levels were normalized against the housekeeping gene GAPDH (NM_001357943.2). PCR primer sequences are summarized as follows: GAPDH(F)5′-TGAAGGTCGGAGTCAACGG-3′,(R): 5′-TCCTGGAAGATGGTGATGGG-3′; TRSP1(F): 5′-AGCCCCAGTAAGGGAGGAAA-3′,(R): 5′-GGGTGCAGGCCATATCTTGAG-3′; PFKFB3(F):5′-GATGCCCTTCAGGAAAGCCT-3′,(R): 5′-TTGAACACTTTTGTGGGGACGC-3′;FADS1(F): 5′-GTCTTCTTCCTGCTGTACCTG-3′,(R): 5′-GGTTCCACTTTGAGGTGCTGA-3′; CPE(F): 5′-CATCTCCTTCGAGTACCACCG-3′,(R): 5′-CCGTGTAAATCCTGCTGATGG-3′.

Gene function analysis and survival analysis

The top 50 hub genes for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was built using R. P < 0.05 was the cutoff criterion, and RUNX2, IBSP, SPP1, SP7 were selected as osteogenic marker genes (Liu et al., 2020). Pearson correlation coefficients were used to indicate the correlation between the marker genes and CPE. The pROC (v 1.17.0.1) analysis was used to compare the predictive accuracy of CPE. The TARGET RNA sequencing data were used to divide patients into low- and high-expression groups according to the median relative gene expression level (50%). Kaplan–Meier survival analysis were conducted with the TARGET cohort in the Xiantao platform.

Analysis of scRNA-Seq transcriptome datasets

The scRNA-Seq transcriptome datasets for GSE162454 were from the NCBI website (https://www.ncbi.nlm.nih.gov/geo). The data quality control process was analyzed using the Seurat package (version 4.1.1). The Merge function was used to merge six OS cases. Quality control was performed using the following criteria: the feature gene number range from 300 to 4,500, and gene counts <15,000 and those with a mitochondrial gene number of <10%. Subsequently the Harmony package (version 1.0) was used to eliminate the batch effect of the cellular data. In addition, the FindClusters function of the Seurat package was used for cell clustering analysis (resolution = 0.091) and visualization by performing a uniform manifold approximation and projection (UMAP) dimension reduction analysis. The cell populations were annotated based on their respective gene expression levels in a known set of genes reported by Liu et al. (2021) (Table S1). The average gene expression was visualized using the VlnPlot and FeaturePlot functions of the Seurat packageS.

Statistical analysis

All analysis methods and R packages were implemented with R version 4.1.3 (R Core Team, 2022). Wilcoxon matched-pairs signed-rank test was used for the comparison between tumor and the adjacent tissues, while Student’s t-test was used to compare the two groups. P-values < 0.05 were considered statistically significant.

Identification of DEGs and hub genes in OS samples

We performed microarray analysis to screen differential genes between OS samples and adjacent tissues. The volcano plot shows that 1,121 genes were upregulated and 983 genes were downregulated (Fig. 1A, File S1). PCA showed sample clustering by DEGs (Fig. 1B), indicating that the tumor tissue and adjacent normal tissue had different expression patterns. We next analyzed the DEGs using GO and enrichment in KEGG pathways using the top 50 hub genes (Figs. 1C and 1E). The most significantly enriched terms in upregulated genes were “Cell cycle” in KEGG (Fig. 1C), “microtubule binding” in GO molecular function, “condensed chromosome” in GO cellular component, and “mitotic nuclear division” in GO biological process. In the downregulated genes, the annotation items were related to skeletal muscle physiological activity. The top 10 upregulated hub genes (KIF11, KIF20A, KIF2C, KIF4A, NCAPG, NDC80, PBK, TOP2A, TPX2, TTK) are presented in the PPI network (Fig. 1D). The top 10 downregulated hub genes were TTN, MYH6, TNNI, MYH8, MYH3, MYL1, TNNT2, TNNI3, TPM2, and MYL3 (Fig. 1F).

We then performed Gene-Set Enrichment Analysis (GSEA) with all genes (Fig. 1G, File S2). These hallmark gene sets were significantly enriched in “reactome_muscle_contraction” (NES = 3.499; p.adjust = 0.032; FDR = 0.021), “reactome_activation_of_ matrix_metalloprotein” (NES = 2.178; p.adjust = 0.032; FDR = 0.021), “PID_PLK1_pathway” (NES = 2.66; p.adjust = 0.032; FDR = 0.021), and “WP_cell_cycle” (NES = 2.642; p.adjust = 0.032; FDR = 0.021), suggesting the progression of osteosarcoma is related to these processes.

Screening for prognostic genes

Next, we obtained 1,104 high-risk factors from the TARGET database by Cox regression analysis (Fig. 2A, File S3). We overlapped the DEGs with the high-risk factors and presented the relative expression levels of the 138 dysregulated prognostic genes in the heatmap (Fig. 2B). These genes were highly enriched for extracellular matrix-related items (Fig. 2C). We validated four abnormal genes in the tumor tissue from the prognostic genes lists using qRT-PCR analysis (Figs. 3A–3D). The results of CPE, TRSP1, and FADS1 were in accordance with the sequencing results (logFC: 2.6767, 1.9593, 2.5813 respectively, P < 0.05). However, the PFKFB3 was inconsistent with high-throughput data which was downregulated in tumor tissue (logFC: −2.369, P < 0.05). Moreover, OS patients with higher expression of CPE, TRSP1, PFKFB3, and FADS1 had worse prognoses than those with lower expression in the TARGET database (Figs. 3E–3H). These results indicated that the CPE, TRSP1, and FADS1 may be associated with the prognosis of osteosarcoma patients and worth further research.

CPE was highly expressed in tumor tissues

To further validate the microarray data, we used the public database to detect CPE mRNA expression (Figs. 4A and 4B). CPE was highly expressed in tumor samples, which was consistent with the results of the microarray (logFC = 2.6767, adj.P.Val = 0.0029). GSE126209 included the 11 paired mRNA profile in OS and adjacent normal tissues. The results showed a significantly higher expression in OS tissues than normal samples (P = 0.001, Fig. 4A). We also analyzed the relationship between CPE expression clinicopathological features in 68 OS patients from GSE33382. Both osteoblastic and chondroblastic OS had higher expression compared to fibroblastic pathological types (P < 0.05, Fig. 4B). We also found a coexpression pattern between the CPE and osteoblastic markers in the TARGET database (Figs. 4C and 4D). The correlation coefficients (P < 0.001) for RUNX2, IBSP, SPP1, and SP7 were 0.43, 0.535, 0.452, and 0.713, respectively. A correlation coefficient of SP7 (r > 0.7) was considered to indicate high correlation.

CPE was an independent prognostic factor for OS

Next, univariate and multivariate Cox regression analyses were performed for overall survival and relapse-free survival to identify independent prognostic factors (Tables S2 and S3). The results indicated that CPE was independent of other clinical factors in relapse-free survival (P = 0.028), while there was no statistical difference in overall survival (P = 0.071), indicating that CPE was correlated with the relapse-free survival of the patients with OS. The AUC value of overall survival was 0.700 (CI [0.591–0.809], Fig. 4E), and the AUC value of relapse-free survival was 0.75 (CI [0.653–0.848], Fig. 4F). The dysregulated genes of CPE had moderate predictive power.

CPE was strongly enriched in osteoblastic OS cells

The scRNA-Seq transcriptome datasets were collected from NCBI for six treatment-naive OS cases to investigate the effects of CPE genes. After normalization and PCA dimensionality reduction, the dimensionality reduction data were displayed in 2D space using Uniform Manifold Approximation and Projection (UMAP) analysis, and the cells were then overlapped into nine clusters (Figs. 5A and 5B). Subsequently, we checked the transcription levels of CPE between the different OS tumor cell types. The average expression values of CPE in all cell types were shown in Figs. 5C and 5D. Our results suggested that CPE was highly expressed in osteoblasts, osteoclast OS cells, and carcinoma-associated fibroblast clusters but expressed at low levels in other cell types.