A total of 160 paired paraffin-embedded CRC specimens and corresponding adjacent non-tumor tissues were collected to design a tissue array chip from the Shanghai Tenth People’s Hospital, Tongji University School of Medicine. The study was approved by the Research Ethics Committee of Shanghai Tenth People’s Hospital and carried out in accordance with the ethical standards formulated in the Helsinki Declaration. The related ethical approval code is 2020-KN155-01. Tissue microarray was constructed by 1.5-mm cores.

Immunohistochemistry (IHC) for PLOD3 was carried out on CRC tissue microarray slides. The slides were first incubated at 60°C for 4 h, deparaffinized in xylene, and then rehydrated in alcohol. After heating in citrate buffer for 23 min, we used 0.3% of hydrogen peroxide (H₂O₂) to block endogenous peroxidase activity. Slides were blocked with 3% bovine serum albumin for 30 min and incubated in the anti-PLOD3 antibody (diluted 1:200; ab128698; Abcam) overnight at 4°C. The next day, after 3 washes with PBS, slides were incubated with secondary antibody for 1 h, then we used the 3,3-diaminobenzidine (DAB) kit for visualization, and hematoxylin was used to stain nuclei. After the experiments, the slides were observed by microscope. All stainings were scored based on the staining intensity and extensity of positive cells, the intensity (0=genitive, 1=weak, 2=moderate, and 3=strong) and extensity (0 = 5% or less of cells stained positive; 1 = 5%-25%; 2 = 26% to 50%; 3 = 51% to 75%; and 4 = 75% or more) of tumor staining were evaluated. The positive cell density of each core was counted by two independent investigators blind to clinical outcome and knowledge of the clinicopathological data. The final IHC score was calculated by multiplying the strongest intensity score and the total extensity score (maximum value of 12).

The Cancer Genome Atlas (TCGA) colon adenocarcinoma (COAD) and rectum adenocarcinoma (READ) data, including gene expression quantified by fragments per kilobase million and clinical information of 51 normal tissues and 638 tumor tissues, were obtained from the UCSC Xena project (http://xena.ucsc.edu/). Four independent validation cohorts (GSE17536, GSE39582, GSE74602, and GSE113513) were obtained from the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo). Gene level mutations (Mutect2) of the COAD and READ cohorts were acquired from TCGA data portal (https://portal.gdc.cancer.gov/).

The TIMER (12) database was used to analyze differences in PLOD3 expression between tumors and normal controls from TCGA data set.

The pan-cancer immune cell infiltration scores for TCGA were obtained from a previously published study (19). The results were based on CIBERSORT (20) and were used for further analysis. TCGA CRC cancer samples were divided into two groups according to median PLOD3 expression (high versus low), and immune cell infiltration was compared between groups. Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (21) was used to estimate a tide score and the predicted response to immune checkpoint blockade. Immunescore and Stromalscore were calculated via “estimate” (22) package.

The top 1000 genes showing the highest correlation with PLOD3 were extracted from the LinkedOmics database (http://www. linkedomics.org) (23). Function annotations were performed to identify potentially involved biological processes and signaling pathways using DAVID (24) 6.8 (https://david.ncifcrf.gov/).

Tumor mutation analysis was performed using the “maftools” package (25). High and low PLOD3 expression oncoplots were generated via “oncoplot” function, and tumor mutation burden (TMB) was calculated using the “maftools” package. Differentially mutated genes between PLOD3-high and low groups were evaluated via Fisher’s exact test.

Four independent validation cohorts (GSE17536 (26), GSE39582 (27), GSE74602, and GSE113513) were obtained from the GEO database and used as validation cohorts to determine the diagnostic and prognostic value of PLOD3. GSE17536 (177 patients) and GSE39582 (585 patients) with relevant survival information were applied for PLOD3’s prognostic value validation; GSE74602 (30 pairs) and GSE113513 (14 pairs) including the CRC tissues and matched adjacent tissues were employed for validating the expression difference of PLOD3 between tumor and normal tissues. In this analysis, patients were divided into low- and high-PLOD3 groups according to an optimal PLOD3 cutoff, which was generated using the association between PLOD3 and survival data with the survminer package. GSE91061 (28) including 105 immunotherapy-treated samples was used for the validation of PLOD3’s therapeutic value. After immunotherapy treatment, samples were classified into the following categories according to the patient’s response: complete response (CR), partial response (PR), stable disease (SD) and progressive disease (PD). Among them, CR and PR are recognized as patients who respond to immunotherapy. SD and PD are recognized as patients who do not respond to immunotherapy. Moreover, the protein levels of PLOD3 in colorectal tumors and normal tissues were assessed using a proteomic dataset (29). Another proteomic dataset (30) was used to validate the prognostic value of PLOD3 according to protein expression levels.

Differences in variables between groups were tested using the Wilcoxon test or chi-squared test, as appropriate. Kaplan-Meier curves were drawn to estimate the overall survival distribution. The log-rank test was used to analyze the statistical difference in survival curves between two groups. Kruskal-Wallis tests were used for comparing PLOD3’s expression difference among more than two comparison groups and Wilcoxon test was used for comparison between two groups. The ROC curve was plotted via “pROC” package. All figures and statistical analyses were performed using R software (version 4.0.2; http://www.R-project.org). A value of p < 0.05 was considered statistically significant. All statistical tests were two-sided.

Analysis of PLOD3 using the TIMER2 database showed that PLOD3 expression was higher in 19 TCGA tumors than in the corresponding normal tissues, including bladder urothelial carcinoma, breast invasive carcinoma, cholangiocarcinoma, COAD, esophageal carcinoma, Glioblastoma multiforme, head and neck squamous cell carcinoma, Kidney Chromophobe, kidney renal clear cell carcinoma, kidney renal papillary cell carcinoma, Liver hepatocellular carcinoma, lung adenocarcinoma, lung squamous cell carcinoma, Pheochromocytoma and Paraganglioma, Prostate adenocarcinoma, READ, Skin Cutaneous Melanoma, Stomach adenocarcinoma, and Uterine Corpus Endometrial Carcinoma ( Figure 1A ). Particularly, higher PLOD3 expression was observed in TCGA COAD and READ cohorts, separately and collectively, compared with the adjacent normal tissues ( Figures 1B, C ), suggesting that PLOD3 plays a role in the pathogenesis of CRC.

The role of PLOD3 in CRC remains unclear. Investigating the correlation between PLOD3 expression and clinical features may clarify the function of PLOD3 in the progression of CRC. In this study, we examined the relationship between PLOD3 expression and the clinical parameters of CRC using TCGA cohort. The results showed that PLOD3 expression differed significantly according to tumor N stage, M stage, clinical stage, and microsatellite instability (MSI) status. Increased N, M, and clinical stages were associated with increased PLOD3 expression (all p < 0.05), suggesting that PLOD3 may be a poor prognostic factor ( Figures 2A–C ). The microsatellite stable (MSS) group also showed higher PLOD3 expression ( Figure 2D ). Survival analysis using PLOD3 median expression as the cut-off value showed that CRCs with higher expression of PLOD3 had a worse prognosis than those with lower expression (log rank p < 0.01) ( Figure 2E ).

After determining the prognostic value of PLOD3 in CRC, we next explored the biological functions associated with PLOD3. GO and KEGG enrichment analyses were performed, and the top GO terms and signaling pathways are shown in Figures 3A, B . PLOD3 gene expression was associated with many biological processes, such as vesicle-mediated transport, ephrin receptor signaling pathway, and autophagy. PLOD3 was also associated with Notch signaling, neurotrophin signaling, and glycosaminoglycan biosynthesis.

To examine the relationship between PLOD3 and mutation profiles in CRC, tumor mutations were compared between PLOD3 high and low groups. We chose the top 20 mutation genes in the whole CRC cohort and compared the mutation frequency difference between PLOD3 high ( Figure 4A ) and low ( Figure 4B ) groups using the Fisher’s test. Furthermore, we visualized these mutation genes in a forest ( Figure 4C ). TP53, APC, and KRAS, were significantly mutated in the PLOD3-high group, whereas PIK3CA, FAT4, and OBSCN were specifically mutated in the low-expression group. In addition, a significant (p < 0.001) negative correlation was observed between PLOD3 and TMB ( Figure 4D ).

Next, we analyzed the immune cell infiltration difference between PLOD3 high and low groups. The infiltration scores of B cell plasma, T cell CD8+, T cell CD4 memory resting, T cell CD4 memory activated, T cell follicular helper, T cell gamma delta, and macrophage M1 were higher in the PLOD3-low cohort than in the PLOD3-high cohort ( Figure 5A ). PLOD3 was significantly negatively correlated with StromalScore ( Figure 5B ) and immuneScore ( Figure 5C ). In addition, PLOD3 was significantly negatively correlated with multiple immune checkpoints ( Figure 5D ) and many other immune related genes, such as antigen-presentation, chemokines, interferons, and T cell inflamed genes ( Figure 5E ). What’s more, a significantly higher TIDE score was observed in PLOD3-high group ( Figure 5F ). In an immunotherapy-treated cohort, patients showed non-response to immunotherapy presented with higher PLOD3 expression ( Figure 5G ). Given that higher PLOD3 expression was associated with lower immunescore, infiltration of multiple immune cells and many immune-related genes, TMB, and MSI score, consistent with the higher TIDE score, we speculated that CRC patients with high PLOD3 expression may be resistant to immunotherapy, which was justified in an immunotherapy-treated cohort.

Two GEO datasets (GSE74602 and GSE113513) were used to validate the diagnostic value of PLOD3 using paired tumor and normal samples. PLOD3 expression was significantly higher in tumor than in normal samples ( Figures 6A, B ). PLOD3 protein expression was higher in CRC than normal tissues ( Figure 6C ). Figure 6D showed PLOD3 protein presented with good diagnostic value between CRC and normal tissues. To further confirm the diagnostic value of PLOD3 in CRC patients, the expression of PLOD3 were analyzed by TMA-based IHC. we compared 160 CRC tumor tissues with paired adjacent normal tissues in a microarray, the representative IHC images of positive PLOD3 expression in tumor tissue and negative PLOD3 expression in normal tissue were shown ( Figure 6E ). Grossly, PLOD3 was overexpressed in tumor parts comparing to normal specimens ( Figure 6F , p-value = 4.2e-14).

GSE17536 and GSE39582 were used to validate the prognostic value of PLOD3. In the GSE17536 dataset, patients with higher PLOD3 expression had a worse disease-free survival ( Figure 6G ). In GSE39582, patients with higher PLOD3 had a worse progression-free survival ( Figure 6H ). The proteomic dataset generated by Li et al. (30) was used to validate the prognostic value of PLOD3 protein expression, and the results revealed that patients with higher PLOD3 protein expression had a worse overall survival and progression-free survival ( Figures 6I, J ). Taken together, the results indicate that PLOD3 is a promising biomarker for the diagnosis, prognosis, and treatment of CRC.