Abstract

   Colorectal cancer (CRC) is a major global health concern, characterized
   by high morbidity and mortality rates. CRC progression involves
   intricate molecular networks that remain incompletely understood. In
   this study, we conducted an integrative multi-omics analysis of
   transcriptomic, proteomic, and metabolomic profiles from CRC tissues
   and matched normal adjacent tissues (NATs). Our analysis revealed 1,394
   differentially expressed long non-Coding RNAs (lncRNAs), 2,788 genes,
   548 proteins, and 91 metabolites. A significant interaction network
   comprising 22 lncRNAs, 14 mRNAs/proteins, and 9 metabolites was
   identified, among which lncRNA 60967.1 emerged as a pivotal regulator.
   Functional validation demonstrated that lncRNA 60967.1 is markedly
   downregulated in CRC cell lines and patient tissues. Overexpression of
   lncRNA 60967.1 restored expression of the tumor suppressor PLCD4 and
   increased levels of all-trans retinoic acid (ATRA). This modulation
   enhanced IFN-γ-induced apoptosis and increased expression of the IFN-γ
   receptor subunit IFNGR1, thereby partially reversing IFN-γ resistance.
   In murine models, lncRNA 60967.1 overexpression promoted immune cell
   infiltration and synergized with anti–PD-1 therapy to inhibit tumor
   growth. Collectively, our findings uncover a novel
   lncRNA-mRNA/protein-metabolite network, the lncRNA 60967.1-PLCD4-ATRA
   axis, that plays a critical role in CRC progression and immune
   modulation, offering promising therapeutic targets for improved
   treatment efficacy.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s12943-025-02359-x.

   Keywords: Colorectal cancer, Multi-omics, LncRNAs, Tumor
   microenvironment, Immunity

Significance

   Colorectal cancer (CRC) is characterized by substantial genetic and
   epigenetic heterogeneity, underscoring the need for novel therapeutic
   targets. While immunotherapy has led to significant advancements in
   cancer treatment, approximately 85% of CRC patients exhibit resistance
   due to different genetic and epigenetic features. Multi-omics
   approaches, which integrate data across genomic, proteomic, and
   metabolomic layers, have emerged as powerful tools for elucidating
   disease mechanisms. In this study, we conducted multi-omics analyses on
   tumor and adjacent normal tissues from 13 CRC patients. Complementary
   in vitro and in vivo experiments demonstrated that lncRNA 60967.1
   regulates the PLCD4/ATRA axis and modulates the immune response to
   anti-PD-1 therapy, thereby promoting CRC progression. Our findings
   reveal a novel regulatory network involving lncRNA, PLCD4, and ATRA,
   providing a potential new target for CRC therapy.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s12943-025-02359-x.

Introduction

   Colorectal cancer (CRC) ranks as the third most frequently diagnosed
   cancer worldwide and is the second leading cause of cancer-related
   deaths [[54]1–[55]3]. While the five-year overall survival rate for
   early-stage localized CRC is high, exceeding 90%, the prognosis
   significantly worsens for metastatic CRC (mCRC), where the five-year
   survival rate drops to around 14% [[56]4, [57]5]. Unfortunately, up to
   25% of CRC cases are already metastatic at diagnosis, and about 50% of
   localized CRC patients eventually develop mCRC. Conventional therapies
   such as surgery, chemotherapy, and radiotherapy offer limited long-term
   survival benefits [[58]6, [59]7]. The heterogeneity of CRC,
   characterized by complex genetic and epigenetic alterations,
   necessitates the identification of novel therapeutic targets to improve
   patient outcomes [[60]8, [61]9].

   In recent years, immunotherapy has transformed the treatment landscape
   for various cancers, including CRC. Immune checkpoint inhibitors (ICI),
   such as pembrolizumab and nivolumab, have shown promising results,
   especially in CRC patients with high microsatellite instability (MSI-H)
   or deficient mismatch repair (dMMR) tumors, which account for
   approximately 15% of stage I-III CRC cases and 3–6% of stage IV CRC
   patients [[62]10]. MSI-H/dMMR tumors harbor numerous somatic mutations
   encoding neoantigens, making them more responsive to immune checkpoint
   inhibitors [[63]3]. However, approximately 95% of mCRC patients are
   microsatellite stable (MSS) or mismatch repair-proficient (MMRp),
   rendering them resistant to immunotherapy [[64]11]. This underscores
   the need for further investigation into the molecular mechanisms
   regulating immune responses in CRC.

   Integrative multi-omics data analysis has emerged as a powerful
   approach to uncover the molecular networks that drive cancer
   progression [[65]12, [66]13]. By integrating high-throughput datasets
   such as transcriptomics, proteomics, and metabolomics, this approach
   provides a more comprehensive understanding of cancer pathogenesis than
   single omics layers alone [[67]14, [68]15]. Multi-omics analysis
   enables the classification of disease subtypes, identification of
   predictive biomarkers, and provides insights into the molecular
   mechanisms underlying diseases. For example, one study profiled 34 CRC
   cell lines using genomic, transcriptomic and proteomic analyses and
   confirmed that these cell lines accurately recapitulate the genomic
   profiles of primary carcinomas and can be stratified into distinct
   molecular subtypes [[69]16]. Similarly, Takei and colleagues conducted
   a multiomic characterization of MSS/MMRp CRC tumors from patients
   treated with ICI combinations, identifying molecular features such as
   upregulated epithelial-mesenchymal transition and G2M checkpoint
   pathways, POLE mutations, and tumor microenvironment factors as
   predictive biomarkers for treatment response, which could aid in the
   development of biomarkers for ICI therapy [[70]17]. In a large cohort
   study involving 100,204 CRC cases and 154,587 controls,
   Fernandez-Rozadilla and colleagues identified 50 novel risk association
   markers. Integrative genomic, transcriptomic, and methylomic analyses
   revealed 53 additional associations and 155 effector genes linked to
   CRC risk, many of which were previously unreported [[71]18].

   To investigate the molecular mechanisms driving CRC progression, we
   conducted a multi-omics analysis integrating transcriptomics,
   proteomics, and metabolomics of CRC tissues and matched normal adjacent
   tissues (NATs) from 13 patients. Although large-scale databases like
   The Cancer Genome Atlas (TCGA) provide a wealth of multi-omics data for
   CRC with a substantial sample size, TCGA primarily employs Reverse
   Phase Protein Array (RPPA) for proteomics analysis and lacks data on
   non-coding RNAs. Moreover, we incorporated metabolomics analysis,
   providing insights into metabolic changes in CRC, which is not covered
   in TCGA’s datasets. In contrast, the uniqueness of this study lies in
   its focus on non-coding RNAs, combined with mass spectrometry-based
   proteomics and metabolomics integration. Additionally, the use of
   consistent samples across all omics data in this study allows for a
   clearer understanding of the interactions between molecules at
   different layers [[72]19, [73]20]. Through hierarchical clustering and
   network analysis, we identified differentially expressed lncRNA
   (DELncs), genes (DEGs), proteins (DEPs), and metabolites (DEMs). Our
   analysis then examined the interactions between lncRNAs,
   mRNAs/proteins, and metabolites, revealing key lncRNAs strongly
   correlated with gene expression and metabolic alterations. Further
   validation using in vitro and in vivo experiments demonstrated the role
   of lncRNA 60967.1 in modulating the PLCD4/ATRA axis and anti-PD-1
   therapy, influencing immune responses and impacting CRC progression.

Materials and methods

Sample collection

   Thirteen pairs of human CRC tissue samples and matched NATs were
   collected immediately after surgical resection at Sir Run Run Shaw
   Hospital, Hangzhou, China (Supplementary Table [74]S1). Sample
   collection adhered to Institutional Review Board (IRB) guidelines and
   was approved by the Ethics Committee (20211204-30). Pathological
   examination confirmed that all specimens were CRC. Patients were
   excluded if they had received prior treatment, had other tumors,
   incomplete clinical data, coexisting blood disorders, or recent
   antibiotic use.

Multi-omics sample preparation and data analysis

   Total RNA was extracted from CRC tissues and NATs using Trizol Reagent
   (Invitrogen Life Technologies). To assess RNA quality, concentration,
   and integrity, we used the Thermo Scientific NanoDrop
   spectrophotometer. Subsequently, ribosomal RNA (rRNA) was removed from
   the total RNA using the Epicentre Ribo-Zero rRNA Removal Kit
   (Human/Mouse/Rat), following the manufacturer’s protocol. Next,
   first-strand cDNA was synthesized using random oligonucleotides as
   primers, and RNA degradation was carried out with RNase H. For the
   second strand, DNA polymerase I was used to synthesize cDNA, with dTTP
   substituted by dUTP. Afterward, the cDNA was purified, end-repaired,
   A-tailed, and adapter-ligated. To eliminate the second cDNA strand
   containing uracil (U), the USER enzyme (NEB, USA) was added.

   To enrich the synthesized double-stranded DNA, a 15-cycle PCR
   amplification was performed using the Illumina PCR Primer Cocktail. The
   PCR products were then purified with the AMPure XP Bead-Based Reagent
   (Beckman Coulter, Beverly, CA, USA). The quality and quantity of the
   library was assessed using the Agilent High Sensitivity DNA Assay on
   the Agilent 2100 Bioanalyzer System and then sequenced on the Illumina
   NovaSeq 6000 platform.

   For gene expression analysis, we used HTSeq (v0.9.1) to calculate the
   read count values for each gene. These values were then normalized
   using fragments per kilobase of transcript per million mapped reads
   (FPKM). In order to assemble the transcripts, we utilized StringTie
   (v2.2.1) along with Hisat2 (v2.1.0) alignment results. Candidate long
   non-coding RNAs (lncRNAs) were identified from the assembled
   transcripts based on characteristic splicing patterns and structural
   features typical of lncRNAs. To further validate the non-coding nature
   of the transcripts, we used three software tools-PLEK (v1.2), CNCI
   (v1.0), and Pfamscan (v1.6.4)-to predict their coding potential. Only
   those transcripts for which all three tools indicated non-coding
   characteristics were classified as lncRNAs. RNA sequencing was carried
   out by Suzhou Panomics Biopharmaceutical Technology, with transcript
   quantification performed using HTSeq (v0.9.1). DEGs were identified
   using DESeq, considering both q-values and fold changes [[75]21].
   Finally, detailed information regarding the HE staining results and
   tumor purity, assessed by pathologists, can be found in Supplementary
   Fig. [76]S1.

   For proteomic analysis, CRC tissues were homogenized using the MP
   FastPrep-24 instrument and lysed in SDT buffer (4% SDS, 100 mM DTT, 150
   mM Tris-HCl, pH 8.0), supplemented with an appropriate amount of
   Protease and Phosphatase Inhibitor Cocktail (EDTA-Free, 10x in
   ddH[2]O). The protein concentration was normalized, and ultrafiltration
   was employed to remove detergents, DTT, and other small molecules. To
   block reduced cysteine residues, iodoacetamide was added, and the
   reaction was incubated in the dark for 30 min. The samples were then
   washed with UA buffer (8 M urea, 150 mM Tris-HCl, pH 8.0), followed by
   NH4HCO3 buffer to further clean the preparations. Subsequently, trypsin
   (Promega) was added to the samples, and protein digestion was carried
   out overnight at 37 °C. The resulting peptides were collected and
   prepared for fractionation using the Scientific Pierce high-pH
   reversed-phase kit. Retention time correction was performed using
   iRT-Kits (Biognosys) according to the manufacturer’s guidelines. For
   Data-Dependent Acquisition (DDA) library generation, all peptide
   fractions were injected into a Thermo Scientific Q Exactive HF-X mass
   spectrometer connected to a Thermo Scientific Easy nLC 1200
   chromatography system. The DDA library data were searched against the
   FASTA sequence database using Spectronaut (v14.4.200727.47784,
   Biognosys).or Data-Independent Acquisition (DIA), peptides from each
   sample were analyzed by liquid chromatography-tandem mass spectrometry
   (LC-MS/MS). A spectral library was generated by importing the original
   raw files and DDA search results into Spectronaut Pulsar X
   (v12.0.20491.4, Biognosys). The DIA data were processed in Spectronaut,
   referencing the constructed spectral library to ensure accurate peptide
   identification and quantification.

   DEPs were identified using a combination of statistical tests. A t-test
   was applied for normally distributed data, while the Wilcoxon rank-sum
   test was used for non-normally distributed data. Principal component
   analysis (PCA) was performed to assess data variability, while Pearson
   correlation and weighted gene co-expression network analysis (WGCNA)
   were used to explore the relationships between mRNA and protein
   expression levels [[77]22].

   For metabolomics analysis, the liquid chromatography (LC) system
   (Thermo Vanquish) was coupled with a mass spectrometer (MS) (Thermo
   QE-HF-X). Initially, the sample was precisely weighed, and an
   extraction mixture composed of 75% methanol and a 9:1
   chloroform-to-water ratio (25%) was added. Steel beads were introduced,
   and the sample was homogenized using a high-throughput tissue
   homogenizer at 50 Hz for 60 s. The mixture was then subjected to
   ultrasonication at room temperature for 30 min, followed by a 30-minute
   incubation on ice. After centrifugation at 12,000 rpm for 10 min at
   4 °C, the supernatant (approximately 850 µL) was collected and
   transferred into a 2 mL centrifuge tube. The solution was dried using a
   vacuum concentrator, and reconstitution was performed with 50%
   acetonitrile solution containing 4 ppm 2-chlorophenylalanine. The
   solution was filtered through a 0.22 μm membrane into a detection vial.
   To create a quality control (QC) sample, 20 µL from each sample was
   pooled. LC-MS analysis was subsequently performed on the remaining test
   samples.

   Metabolomics data analysis was conducted using OmicStudio
   ([78]https://www.omicstudio.cn) and OECloud
   ([79]https://cloud.oebiotech.com) platforms [[80]23]. Partial least
   squares discriminant analysis (PLS-DA) was used to generate the PLS
   score plot and model validation plot in R (v4.0.3). A t-test was
   employed to identify differentially expressed metabolites (DEMs),
   followed by functional pathway enrichment analysis using MetaboAnalyst
   (27). Stamp plot analysis was generated using OmicStudio, and receiver
   operating characteristic (ROC) curves were constructed to assess the
   potential of DEMs for CRC identification. Pearson correlation was
   applied to explore the relationships between DEMs and mRNAs/proteins.

Cell lines and reagents

   CRC cell lines, including DLD1, HCT8, HCT116, and SW480, were obtained
   from Prosai, while NCM460 and MC38 cell lines were sourced from Warner
   Bio. All cell lines were cultured in either DMEM or RPMI-1640 medium,
   each supplemented with 10% fetal bovine serum (FBS) and 1%
   penicillin/streptomycin. Authentication of all cell lines used in this
   study was performed through STR (short tandem repeat) profiling to
   ensure their identity and integrity.

Isolation of RNA, reverse transcription, and quantitative RT-PCR (qRT-PCR)
analysis.

   Total RNA was isolated from samples using TRIzol reagent (Invitrogen),
   according to the manufacturer’s instructions. cDNA was synthesized from
   1 µg of RNA using the PrimeScript RT Master Mix (Takara). qRT-PCR was
   performed with TB Green Premix Ex Taq II (Takara) following the
   manufacturer’s protocols. The amplification conditions were optimized
   for each target gene. Primer sequences used in this study are provided
   in Supplementary Table [81]S2.

Western blot analysis

   Cells were lysed in RIPA buffer (Sigma-Aldrich) containing protease
   inhibitors to extract total protein. Protein samples were separated by
   SDS-PAGE and transferred to PVDF membranes (Millipore). Membranes were
   blocked with 2.5% nonfat milk and incubated overnight at 4 °C with
   primary antibodies. After washing, membranes were incubated with
   appropriate secondary antibodies. Protein bands were detected using an
   enhanced chemiluminescence (ECL) detection system (Thermo Fisher
   Scientific). Details of the antibodies used, including their sources
   and dilutions, are provided in Supplementary Table [82]S3.

Hematoxylin and Eosin (H&E) and immunohistochemistry (IHC) analysis

   CRC cell lines (DLD1, HCT8, SW480) and normal colon cells (NCM460) were
   fixed in 4% paraformaldehyde. The fixed cells were subjected to
   standard H&E staining or IHC staining using the SP Rabbit & Mouse HRP
   Kit (CW2069S, Beijing, China) with hematoxylin as a counterstain.
   Primary antibodies against PLCD4, CD4, and CD8 were applied according
   to the manufacturer’s instructions. Antibody details, including
   dilution factors, are provided in Supplementary Table [83]S3. For IHC
   analysis, the staining intensity and the percentage of positively
   stained cells were evaluated semi-quantitatively. The percentage of
   positive cells was categorized into the following four grades: 0
   (< 5%), 1 (6-25%), 2 (26-50%), 3 (51-75%), and 4 (> 75%). The intensity
   of staining was classified into four levels: 0 (negative), 1 (weak), 2
   (moderate), and 3 (strong). For each tissue section, ten
   high-magnification fields (x400) were analyzed, and the average score
   was calculated based on the intensity and percentage of positive cells
   [[84]24, [85]25].

Enzyme-linked immunosorbent assay (ELISA)

   ELISA assays were performed using culture supernatants from cell
   cultures, and detection was carried out using commercially available
   kits from Shanghai Enzyme-linked Biotechnology Co., Ltd (Shanghai). For
   lentiviral transfection, cells were seeded at a density of 3–5 × 10⁴
   cells/mL in 12-well plates. When cells reached 30–50% confluence,
   polybrene (7 µg/mL) was added to enhance viral infection, and cells
   were infected with lentiviral particles from General Biology Anhui Co.,
   Ltd. Transfection efficiency was assessed 72 h post-infection by
   measuring the expression of the transgene.

RNA fluorescence in situ hybridization (RNA FISH)

   For RNA FISH (GimaGen), cells were seeded on confocal dishes and fixed
   with 4% paraformaldehyde. After fixation, the cells were blocked to
   reduce non-specific binding. Hybridization with specific probes was
   performed overnight at 37 °C. Nuclei were counterstained with DAPI.
   Images were captured using a laser scanning confocal microscope to
   visualize RNA localization and expression.

Analysis of PLCD4 expression and tumor microenvironment (TME) infiltration
atterns

   The relationship between PLCD4 expression and TME infiltration patterns
   in CRC was systematically investigated using data from the TCGA CRC
   cohort. Standardized gene expression profiles (TPM) for both TCGA Colon
   Adenocarcinoma (TCGA-COAD) and TCGA Rectal Adenocarcinoma (TCGA-READ)
   were downloaded via the “GDCquery” function of the “TCGAbiolinks” R
   package. Prior to merging the mRNA expression matrices of TCGA-COAD and
   TCGA-READ, potential batch effects were assessed using the “draw_pca”
   function of the “tinyarray” R package. No significant batch effects
   were detected (Supplementary Fig. [86]S7A). The two datasets were then
   combined, resulting in a unified mRNA expression matrix for the
   TCGA-CRC cohort, which included 650 CRC samples.

   To assess the TME, the Estimation of Stromal and Immune Cells in
   Malignant Tumors Using Expression Data (ESTIMATE) algorithm (1) was
   applied to calculate the ImmuneScore, StromalScore, and Tumor Purity
   for each sample. Pearson’s correlation analysis was performed to
   examine the relationship between PLCD4 mRNA expression levels
   (log2(TPM + 1)) and these scores. Additionally, to validate the
   findings from the ESTIMATE algorithm, the xCell algorithm (2) was
   employed using the “xCell” R package. The xCell-based microenvironment
   score represents the sum of all immune and stromal cell types within
   the TME.

RNA immunoprecipitation (RIP) assay for PLCD4-RNA complexes

   RIP was performed using the RNA Immunoprecipitation Kit (Cat. No.
   P0102, Geneseed, Guangzhou) according to the manufacturer’s
   instructions. NCM460 cells were collected and lysed, and 100 µL of the
   lysate was saved as the input sample. To capture the PLCD4-RNA complex,
   human PLCD4 antibody (30133-1-AP, Proteintech) and IgG antibody were
   incubated with Protein A + G beads overnight at 4 °C. The
   immunoprecipitated RNA-protein complexes were then used for downstream
   analysis, including qRT-PCR to quantify RNA levels and Western blot
   analysis to detect protein expression.

In vivo tumor xenograft and treatment models

   Immunocompetent and nude BALB/c mice (4–6 weeks old) were obtained from
   Jiangsu Jicui Yaokang Biotechnology Co., Ltd. All animal experiments
   were approved by the Institutional Animal Care and Use Committee
   (IACUC) of the respective institutions. For the human CRC xenograft
   model, 6-week-old BALB/c nude mice were subcutaneously injected with
   5 × 10⁶ HCT116 cells transfected with either a control or
   PLCD4-overexpressing vector (n = 6 per group). Tumor size was measured
   every 2–3 days starting on day 5, and tumor volume was calculated using
   the formula: (length × width²)/2.

   For the mouse CRC model, 5 × 10⁶ CT26 cells, transfected with either
   GFP or lncRNA 60967.1 vectors, were injected subcutaneously. Tumor size
   was measured every 4 days using the same formula. After 20 days, tumors
   were harvested, weighed, and embedded in paraffin for further analysis.

   For PD-1 and lncRNA-targeting treatment, 250 µg of anti–PD-1 monoclonal
   antibody (RMP1-14, BioXCell) was administered via intraperitoneal
   injection once a week, starting 5 days after the establishment of the
   CT26 mouse model.

Statistical analysis and data interpretation

   Statistical significance was assessed using GraphPad Prism 8 software
   and the R programming language. All assays were performed in
   triplicate, with results presented as the mean ± SEM unless otherwise
   indicated. P-values were calculated using two-tailed Student’s t-tests.
   Pearson’s correlation coefficient was used for correlation analysis to
   assess the relationship between mRNA and protein module eigengenes and
   traits, with a correlation coefficient threshold of ≥ 0.3 and a
   p-value < 0.05 for selection. We performed WGCNA on mRNA and proteins,
   setting the power value to 1 for mRNA and 14 for proteins based on the
   distribution and connectivity of the scale-free network. PLS-DA was
   conducted in R (v4.0.3) to evaluate the PLS score and model validation
   plots. ROC analysis of metabolites was performed using OmicStudio or
   OECloud, with p-values < 0.05 considered statistically significant.
   ELISA data were similarly analyzed using two-tailed Student’s t-tests.
   Further details on statistical methods are provided in the figure
   captions.

Results

Integrative multi-omics analysis reveals novel lncRNA-mRNA/protein-metabolite
crosstalk driving CRC progression

   To investigate the mechanisms underlying CRC progression, we performed
   a comprehensive multi-omics analysis, including transcriptomics,
   proteomics, and metabolomics, on CRC tissues and matched NATs from 13
   patients. We analyzed and validated datasets consisting of 26
   transcriptomes, 22 proteomes (11 patients), and 66 metabolomes (11
   patients, triplicate per sample). Initially, 2,788 DEGs and 818 DEPs
   were identified (Fig. [87]1A and [88]B). Unsupervised hierarchical
   clustering of mRNA and protein expression data revealed distinct
   expression patterns between CRC and NAT samples (Fig. [89]1C and
   [90]D).

Fig. 1.

   [91]Fig. 1
   [92]Open in a new tab

   Identification of hub genes using WGCNA and validation against external
   databases. (A) Volcano plot of mRNA expression in 13 CRC tissues and
   paired NATs. (B) Volcano plot of protein expression in 11 CRC tissues
   and paired NATs. (C) Hierarchical heatmap of DEGs in 13 CRC tissues and
   paired NATs. (D) Hierarchical heatmap of DEPs in 11 CRC tissues and
   paired NATs. (E) Venn diagram illustrating the overlap between mRNAs
   and proteins with consistent expression patterns. (F) WGCNA
   co-expression module analysis of mRNAs. (G) GSEA of 52 hub genes in the
   turquoise module. (H) Survival analysis of the gene set in the
   turquoise module. Abbreviations: DEGs, differentially expressed genes;
   DEPs, differentially expressed proteins; WGCNA, weighted gene
   co-expression network analysis; GSEA, gene set enrichment analysis

   By integrating transcriptomic and proteomic data, we identified 106
   upregulated and 124 downregulated genes in CRCs compared to NATs,
   demonstrating a consistent expression trend between mRNA and protein
   expression (Fig. [93]1E). After excluding genes with low and unstable
   protein expression, 27 upregulated and 52 downregulated genes with|log2
   FC| >2 were selected for further analysis. These 79 genes formed two
   co-expression modules based on WGCNA analysis: M1 (52 genes highlighted
   in turquoise) and M2 (27 genes labeled in blue for mRNAs and gray for
   proteins) (Fig. [94]1F and Supplementary Fig. [95]S2A) [[96]22].
   Pearson correlation analysis was performed to assess the relationship
   between mRNA and protein module eigengenes and traits, selecting
   significant associations with a correlation coefficient ≥ 0.3 and
   p-value < 0.05, with mRNA and protein modules selected separately (see
   Supplementary Fig. [97]S2B and [98]S2C). The M1 module, which showed
   strong co-expression at both mRNA and protein levels (Supplementary
   Fig. [99]S2D, and [100]S2E), was found to be associated with tumor cell
   proliferation, metastasis, amino acid metabolism, and signal
   transduction, as determined by Gene Set Variation Analysis (GSVA)
   [[101]26] (Fig. [102]1G). Further analysis using the GEPIA2 [[103]27]
   and CPTAC2 [[104]28] databases confirmed the decreased RNA and protein
   expression in CRC tissue compared to normal colon tissue for 42 genes
   in the M1 module (Supplementary Fig. [105]S2F). Additionally, survival
   analysis of Disease-Free Survival (RFS) using the GEPIA2 database
   showed that altered mRNA expression levels of 51 genes in the M1 module
   were associated with poor prognosis in CRC patients, with 50% high and
   50% low expression as cutoff values for grouping (Fig. [106]1H).

   We further explored the role of lncRNAs in regulating the mRNA/protein
   co-expression module (M1). Of the 1,394 DELncs, 436 showed|Log2FC| >5
   (Supplementary Fig. [107]S3A). Then, we examined the co-expression
   correlation between 436 lncRNAs and the 52 mRNA/protein expressions in
   the M1 module (Supplementary Fig. [108]S3B and [109]S3C). Correlation
   analysis between these lncRNAs and the 52 mRNA/protein expressions in
   the M1 module revealed a significant association for 174 lncRNAs
   (P < 0.05) (Fig. [110]2A). Among them, 22 candidate lncRNAs, including
   2 known and 20 novel lncRNAs, were identified as potential regulators
   based on a correlation threshold (r > 0.9) between lncRNA and mRNA
   [[111]29], as well as a significant correlation with protein expression
   (Fig. [112]2B).

Fig. 2.

   [113]Fig. 2
   [114]Open in a new tab

   Integrated multi-omics analysis reveals network associations among
   lncRNAs, mRNAs, proteins, and metabolites. (A) Upset plot of lncRNAs
   showing statistically significant associations with mRNAs and proteins.
   (B) Network diagram depicting lncRNA interactions with target genes.
   (C) PLS-DA analysis of metabolomic data. (D) Evaluation of the PLS-DA
   model for metabolomic data, including model performance metrics and the
   separation of sample groups. (E) ROC analysis identifying two
   metabolites with AUC > 0.8. (F) Venn diagram illustrating the overlap
   among mRNA–metabolite pairs, protein–metabolite pairs, the top 30 AUC
   metabolites, and the top 30 STAMP metabolites. (G) Bioinformatic
   analysis reveals predicted associations among lncRNAs, mRNAs/proteins,
   and metabolites. Top: data analysis pipeline; Bottom: network of
   lncRNAs, mRNAs/proteins, and metabolites. Abbreviations: ROC, receiver
   operating characteristic; AUC, area under the curve; PLS-DA, partial
   least squares discriminant analysis; STAMP, statistical analysis of
   taxonomic and functional profiles

   To better understand the metabolic changes in CRC, we performed
   metabolomic profiling on 11 CRC tissues and matched NATs, each with
   three replicates. PLS-DA analysis revealed that metabolic features
   could effectively distinguish CRC tissues from NATs (Fig. [115]2C),
   with good predictive ability (R²=0.685, Q²=0.508) (Fig. [116]2D)
   [[117]30]. We identified 91 DEMs (P < 0.05 and|log2FC|>0)
   (Supplementary Fig. [118]S3D), which were enriched in pathways related
   to amino acid metabolism, according to functional analysis using
   MetaboAnalyst (Supplementary Fig. [119]S3E) [[120]31]. Through ROC
   analysis, we found 30 DEMs with AUC values > 0.7; among these, ATRA and
   13,14-dihydro-15-keto-PGE2 were significantly reduced in CRC tissues
   compared to NATs (AUC = 0.832 and 0.818, respectively) (Fig. [121]2E).
   Using the STAMP approach [[122]32], we identified and selected the top
   30 DEMs (Supplementary Fig. [123]S3F), and then explored the
   correlation between these two sets of DEMs and gene or protein
   expression in the M1 module (Supplementary Fig. [124]S3G and H). Our
   findings revealed that 17 DEMs were significantly correlated with 40
   genes in the M1 module both at mRNA and protein levels, 15 of which
   were among the top 30 DEMs identified by both STAMP and ROC analysis
   (Fig. [125]2F). Furthermore, we identified 10 metabolites from 15
   candidate markers that showed significant correlations with 29
   molecular entities at both the mRNA and protein levels, suggesting
   their potential as reliable biomarkers for CRC.

   Taken together, by linking lncRNAs to metabolites through mRNA and
   protein interactions, we uncovered a complex interaction network
   involving 22 lncRNAs, 14 mRNAs/proteins, and 9 metabolites
   (Fig. [126]2G).

In vitro validation of multi-level molecular expression

   To validate the associations among lncRNAs, mRNAs/proteins, and
   metabolites, we focused on a specific network branch involving a novel
   lncRNA 60967.1, its potential target genes, and related metabolites
   (Fig. [127]2G). This lncRNA exhibited significant associations with 12
   out of 14 genes and 6 out of 9 metabolites within the network,
   suggesting its potential roles in CRC. Sequencing data revealed that
   lncRNA 60967.1 is 370 bp in length and located on chromosome 14
   (Supplementary Table [128]S4). We evaluated its expression by qRT- PCR
   in three CRC cell lines, a human colon mucosal epithelial cell line
   (NCM460), and tumor and adjacent normal tissues from 12 additional CRC
   patients (Supplementary Table [129]S5). Notably, the expression of
   lncRNA 60967.1 was significantly downregulated in CRC cell lines
   compared to the normal control (Fig. [130]3A). This finding was further
   corroborated by qRT-PCR analysis of the 12 CRC tissue samples and
   matched NATs (Fig. [131]3B).

Fig. 3.

   [132]Fig. 3
   [133]Open in a new tab

   Experimental validation of mRNA, lncRNA, and metabolite levels. (A)
   Expression of lncRNA 60967.1 in CRC cell lines (DLD1, HCT8, SW480)
   relative to the normal colon cell line (NCM460), measured by qRT-PCR.
   (B) qRT-PCR validation of lncRNA 60967.1 expression in 12 paired CRC
   tumor and normal tissues. (C) qRT-PCR analysis of PLCD4 mRNA levels in
   the three CRC cell lines (DLD1, HCT8, SW480) and the normal colon cell
   line. (D) qRT-PCR validation of PLCD4 mRNA levels in 12 paired CRC
   tumor and normal tissues. (E-F) Western blot analysis of PLCD4 protein
   expression in the three CRC cell lines (DLD1, HCT8, SW480) and the
   normal colon cell line (NCM460) (E) and 12 paired tumor (T#) and normal
   (N#) tissues (F). (G) IHC analysis of PLCD4 protein expression in the
   three CRC cell lines (DLD1, HCT8, SW480) and the normal colon cell line
   (NCM460). (H) ELISA-based quantification of ATRA levels in culture
   supernatants from the three CRC cell lines (DLD1, HCT8, SW480) and the
   normal colon cell line (NCM460). ***P < 0.01 (t test)

   We then performed qRT-PCR analysis of the mRNA levels of
   lncRNA-targeted genes and identified that the mRNA levels of three
   target genes, CAP2, GNAO1, and PLCD4, were significantly decreased in
   the three CRC cell lines (DLD1, HCT8, SW480) compared to NCM460
   (Fig. [134]3C, Supplementary Fig. [135]S4A and [136]S4B), as well as in
   12 CRC tissues relative to matched NATs (Fig. [137]3D, Supplementary
   Fig. S4C and S4D). Next, we examined the protein expression of CAP2,
   GNAO1, and PLCD4 in both cultured CRC cells and tissue samples using
   Western blotting. Among the three, only PLCD4 consistently showed a
   significant reduction in protein levels in both CRC cells and tissue
   samples (Fig. [138]3E and [139]F). Immunohistochemistry (IHC) analysis
   further confirmed the decreased expression of PLCD4 protein in CRC
   tumor cells (Fig. [140]3G).

   Notably, PLCD4, a member of the phosphatidylinositol-specific
   phospholipase C (PLC) family, hydrolyzes
   phosphatidylinositol-4,5-bisphosphate (PIP[2]) into inositol
   1,4,5-trisphosphate (IP[3]), a second messenger that triggers calcium
   (Ca²⁺) mobilization, and diacylglycerol, which activates protein kinase
   C (PKC), a key component in signaling pathways [[141]33–[142]35]. We
   then performed a pathway analysis of PLCD4 using the KEGG database
   ([143]https://www.kegg.jp). To evaluate the levels of PIP[2], IP[3],
   and PKC, we performed IHC assays on CRC cell lines and the normal
   intestinal epithelial cell line NCM460. The expression levels of PIP[2]
   and its hydrolyzed product, IP[3], were not significantly downregulated
   in the CRC cell lines compared to the control. Similarly, no
   significant differences in PKC protein levels were observed between the
   CRC cell lines and the control cells. These results suggest that PLCD4
   may exert its function through alternative signaling pathways
   (Supplementary Fig. [144]S5A and [145]S5B).

   Previous studies have reported significant downregulation of ATRA
   synthesis enzymes in CRC, suggesting reduced ATRA levels in CRC tumors
   [[146]36, [147]37]. To explore this, we assessed the levels of the
   metabolite ATRA in the cell-free supernatant of cultured CRC cell lines
   and NCM460 cells using an ELISA assay. Consistent with previous
   findings, we observed a significant reduction in ATRA levels in the CRC
   cell lines compared to normal control cells (Fig. [148]3H).

   Overall, our cellular and tissue-based validation experiments revealed
   a critical regulatory axis comprising lncRNA 60967.1, PLCD4, and the
   metabolite ATRA. These findings reinforce the reliability of the
   identified network and suggest a potential role for these biomolecules
   in CRC pathogenesis.

LncRNA 60967.1 modulates the PLCD4/ATRA axis and reverses IFN-γ resistance in
CRC cells

   To investigate the regulatory role of lncRNA 60967.1 on the PLCD4/ATRA
   axis, we first examined its physical localization relative to PLCD4
   mRNA using fluorescence in situ hybridization (FISH). As shown in
   Fig. [149]4A, lncRNA 60967.1 and PLCD4 mRNA were colocalized in the
   cytoplasm. Next, we evaluated whether lncRNA 60967.1 interacts with
   PLCD4 protein by performing an RIP assay in NCM460 cells, followed by
   qRT-PCR and Western blot analysis. The results demonstrated that lncRNA
   60967.1 binds to PLCD4 protein significantly more than to the IgG
   control (Fig. [150]4B), and Western blotting confirmed the enrichment
   of PLCD4 in the immunoprecipitated complex (Fig. [151]4C). These
   findings suggest a direct interaction between lncRNA 60967.1 and PLCD4.
   Subsequently, we performed lentiviral expression of lncRNA 60967.1 in
   CRC cells. Overexpression of lncRNA 60967.1 led to increased levels of
   PLCD4 at both the mRNA (Fig. [152]4D) and protein (Fig. [153]4E)
   levels, as well as elevated concentrations of the metabolite ATRA
   (Fig. [154]4F).

Fig. 4.

   [155]Fig. 4
   [156]Open in a new tab

   Molecular and functional effects of lncRNA 60967.1 on PLCD4 expression,
   ATRA levels, and CRC cell proliferation. (A) Fluorescence in situ
   hybridization (FISH) showing the localization of lncRNA 60967.1 and
   PLCD4 mRNA in NCM460 cells. (B) RIP assay followed by qRT-PCR analysis
   of lncRNA 60967.1 in immunoprecipitated (IP) and IgG control groups.
   (C) RIP assay followed by Western blot analysis of PLCD4 in the input,
   immunoprecipitated (IP), and IgG control groups. (D) qRT-PCR analysis
   of PLCD4 mRNA in CRC cells (DLD1, HCT8, SW480) transfected with GFP
   control or lncRNA 60967.1 vector. (E) Western blot analysis of PLCD4
   protein expression in CRC cells (DLD1, HCT8, SW480) with GFP or lncRNA
   60967.1 vector. (F) ELISA measurement of ATRA levels in CRC cells
   (DLD1, HCT8, SW480) with GFP or lncRNA 60967.1 vector. (G-I) CCK-8
   assay for CRC cell proliferation (DLD1, HCT8, SW480) transfected with
   lncRNA 60967.1 vector, followed by IFN-γ stimulation. (J) qRT-PCR
   analysis of IFNGR1/IFNGR2 mRNA in CRC cells transfected with GFP or
   lncRNA 60967.1 vector. (K) Western blot analysis of IFN-γRa protein in
   CRC cells (DLD1, HCT8, and SW480) transfected with GFP or lncRNA
   60967.1 vector. (L) qRT-PCR evaluation of IFNGR1/IFNGR2 mRNA in CRC
   cells (DLD1, HCT8, and SW480) stimulated by ATRA or DMSO. Abbreviation:
   CCK-8, Cell Counting Kit-8. NS: not significant (t-test); **P < 0.05,
   ***P < 0.01 (t-test)

   ATRA is a well-established metabolite that regulates key cellular
   processes, including cell growth, differentiation, and apoptosis and
   modulates immune responses [[157]38–[158]40]. However, our in vitro
   functional assays showed that overexpression of lncRNA 60967.1, and the
   resulting increase in ATRA levels, did not significantly affect CRC
   cell proliferation or apoptosis (data not shown). Given that DLD1 and
   SW480 cell lines are known to be resistant to IFN-γ, a critical
   cytokine involved in activating immune responses [[159]41], and that
   synergistic effects of ATRA and IFN-γ have been reported in animal
   models and cultured cells [[160]42–[161]44], we investigated whether
   the lncRNA 60967.1-regulated PLCD4/ATRA axis could modulate IFN-γ
   resistance in CRC cells. Notably, overexpression of lncRNA 60967.1
   significantly reduced CRC cell proliferation in response to IFN-γ
   stimulation (Fig. [162]4G-[163]I) and increased IFN-γ-induced apoptosis
   (Supplementary Fig. [164]S6A-[165]F). While CRC cells transfected with
   control lentivirus only had no significant impact on IFN-γ stimulation
   (Supplementary Fig. [166]S6G).

   We next investigated the expression levels of the IFN-γ receptor
   subunits, IFNGR1 and IFNGR2. Overexpression of lncRNA 60967.1 in CRC
   cells resulted in a significant increase in IFNGR1 mRNA levels. In
   contrast, IFNGR2 mRNA expression showed variability across the three
   CRC cell lines tested (Fig. [167]4J). This upregulation of IFNGR1 was
   further confirmed at the protein level (Fig. [168]4K). To explore the
   potential difference in IFNGR1 expression between CRC and normal cells,
   we compared IFNGR1 expression in CRC cells and NCM460 cells. The
   results indicated that IFNGR1 expression was significantly lower in CRC
   cells compared to NCM460 cells (Supplementary Fig. [169]S6H-[170]I).
   Additionally, we observed that ATRA stimulation in CRC cells induced
   the upregulation of IFNGR1 expression, although IFNGR2 mRNA expression
   remained variable across the different CRC cell lines (Fig. [171]4L).

   These findings suggest that lncRNA 60967.1 regulates the PLCD4/ATRA
   axis in CRC cells and contributes to the partial reversal of IFN-γ
   resistance by upregulating the IFN-γ receptor subunit IFNGR1, thereby
   enhancing the response to immune stimulation.

PLCD4 acts as a tumor suppressor in CRC by inhibiting migration, invasion,
and tumor growth

   PLCD4 is a critical component of the lncRNA/PLCD4/ATRA axis. To
   investigate its role in CRC development, we generated a PLCD4
   overexpression vector and transfected it into two CRC cell lines, DLD1
   and HCT116. In vitro migration and invasion assays revealed that PLCD4
   overexpression significantly inhibited the invasive and migratory
   capacities of CRC cells compared with cells transfected with the empty
   vector control, suggesting that PLCD4 reduces cell motility and may
   impede tumor progression (Fig. [172]5A). Furthermore, HCT116 cells
   overexpressing PLCD4 were injected into Balb/c nude mice. Analysis of
   tumor volume and weight showed a significant reduction in the PLCD4
   overexpression group compared to the vector control group, confirming
   the tumor suppressive effect of PLCD4 (Fig. [173]5B and [174]C, and
   [175]5D).

Fig. 5.

   [176]Fig. 5
   [177]Open in a new tab

   PLCD4 acts as a tumor suppressor in CRC. (A) Migration and invasion
   assays using two CRC cell lines (HCT116 and DLD1) transfected with
   either GFP (NC) or PLCD4 expression (PLCD4 OE) vectors. (B) Tumors from
   Balb/c nude mice injected with HCT116 cells transduced with either GFP
   (Top) or PLCD4 (Bottom) expression vectors. (C) Tumor growth over time
   in Balb/c nude mice injected with HCT116 cells transduced with GFP
   (blue) or PLCD4 (pink) expression vectors. (D) Measurement of tumor
   weight in Balb/c nude mice injected with HCT116 cells transduced with
   GFP (blue) or PLCD4 (pink) expression vectors. (E) Western blot
   analysis of PLCD4 expression in CRC cell lines (DLD1, SW480, HCT116)
   following siRNA-mediated knockdown (siPLCD4#1 or siPLCD4#2) or
   non-targeting control siRNA (siControl). (F) Migration and invasion
   assays using CRC cell lines (DLD1, SW480, HCT116) transfected with two
   PLCD4-targeting siRNAs (siPLCD4#1 or siPLCD4#2) or non-targeting
   control siRNA (siControl)

   Conversely, two distinct siRNAs were used to knock down PLCD4 in three
   CRC cell lines. Western blot analysis confirmed a significant reduction
   in PLCD4 expression in all three cell lines (Fig. [178]5E). Notably,
   migration and invasion assays showed that, compared to cells
   transfected with the empty vector control, PLCD4 knockdown
   significantly enhanced the migratory and invasive potential of CRC
   cells (Figs. [179]5F).

   Lastly, we calculated the immune score, stromal score, and tumor purity
   for 650 TCGA-CRC samples using the ESTIMATE algorithm, which allowed us
   to gain insights into the infiltration patterns within the CRC TME.
   Pearson’s correlation analysis revealed a significant positive
   correlation between PLCD4 mRNA expression and both the ImmuneScore and
   StromalScore, indicating enhanced infiltration of immune and stromal
   cells within the TME. In contrast, a negative correlation was observed
   between PLCD4 expression and Tumor Purity, suggesting a reduction in
   tumor cell content as immune and stromal components increase
   (Supplementary Fig. [180]S7B). To validate these results, a similar
   analysis was conducted using the xCell algorithm, which corroborated
   the findings from the ESTIMATE algorithm (Supplementary Fig. [181]S7C).
   These findings suggest that PLCD4 may play a crucial role in modulating
   the TME by promoting the infiltration of stromal and immune cells,
   potentially influencing disease progression and therapeutic outcomes.

   Taken together, these findings demonstrate that PLCD4 acts as a tumor
   suppressor in CRC by inhibiting cell motility and tumor growth.

LncRNA 60967.1 suppresses tumor development by modulating the tumor
microenvironment via PLCD4 /ATRA regulation

   To assess the impact of lncRNA 60967.1 overexpression on the tumor
   microenvironment, we transplanted CT26 and MC38 murine CRC
   cells—transfected with either an lncRNA 60967.1 expression vector or a
   control vector—into immunocompetent BALB/c or C57BL/6 mice.
   Overexpression of lncRNA 60967.1 led to a significant increase in both
   PLCD4 and IFNGR1 expression at the mRNA and protein levels in the mouse
   tumor tissues (Fig. [182]6A1, [183]6A2, [184]6B1 and [185]6B2). In
   parallel, lncRNA 60967.1 overexpression also resulted in elevated ATRA
   levels (Fig. [186]6C1 and [187]6C2). In addition, mice bearing tumors
   derived from cells overexpressing lncRNA 60967.1 exhibited significant
   reductions in tumor volume (Fig. 6D1 and 6D2), proliferation rate
   (Fig. 6E1 and 6E2), and tumor weight (Fig. [188]6F1 and [189]6F2)
   compared to control groups.

Fig. 6.

   [190]Fig. 6
   [191]Open in a new tab

   LncRNA 60967.1 overexpression promotes PLCD4 and ATRA expression,
   enhances immune cell infiltration, and reduces tumor growth in mouse
   CRC models. (A1–A2) qRT-PCR analysis of lncRNA 60967.1, PLCD4, and
   IFNGR1/2 expression in tumor tissues from mice injected with mouse CRC
   cells (CT26, A1; MC38, A2) transfected with lncRNA 60967.1 or GFP
   vectors. (B1–B2) Western blot analysis of PLCD4 and IFN-γRα proteins in
   tumor tissues from mice injected with CT26 cells (B1) or MC38 cells
   (B2) transfected with lncRNA 60967.1 or GFP vectors. (C1–C2) ELISA
   analysis of ATRA in tumor tissues from mice injected with CT26 cells
   (C1) or MC38 cells (C2) transfected with lncRNA 60967.1 or GFP vectors.
   (D1–D2) Tumors harvested from Balb/c mice that were subcutaneously
   injected with either CT26 cells (D1) or MC38 cells (D2). Each cell line
   was transfected with GFP or lncRNA 60967.1 vectors prior to injection,
   allowing comparison of tumor growth and size between the control and
   experimental groups. (E1–E2) Tumor growth curves in Balb/c mice
   injected with CT26 cells (E1) or MC38 cells (E2) transfected with
   either GFP or lncRNA 60967.1 vectors. Tumor sizes were measured at
   regular intervals to compare progression between the control (GFP) and
   lncRNA 60967.1 overexpression groups. (F1–F2) Measurement of tumor
   weight in Balb/c nude mice injected with CT26 cells (F1) or MC38 cells
   (F2), transfected with either GFP or lncRNA 60967.1 expression vectors.
   (G) IHC analysis of CD4 and CD8 expression in tumor tissues from Balb/c
   mice subcutaneously injected with either CT26 or MC38 cells. Each cell
   line was transfected with GFP or lncRNA 60967.1 vectors before
   injection, enabling a comparison of immune cell infiltration between
   control and lncRNA 60967.1-overexpressing groups. (H) IHC analysis of
   IFN-γ, PD-1, PD-L1, and Ki67 expression in tumor tissues from Balb/c
   mice subcutaneously injected with either CT26 or MC38 cells. Each cell
   line was transfected with GFP or lncRNA 60967.1 vectors before
   injection, allowing for a comparison of the functional states or
   activation/exhaustion profiles of immune cells between the control and
   lncRNA 60967.1-overexpressing groups. NS indicates no significant
   difference (t-test); **P < 0.05, ***P < 0.01 (t-test)

   To further investigate the impact on immune cell infiltration, IHC was
   performed to evaluate CD4⁺ and CD8⁺ cell populations in the tumor
   tissues. As expected, a significant increase in both CD4⁺ and CD8⁺ T
   cells (Fig. [192]6G) was observed in the CT26-Lnc group compared to the
   CT26-GFP group, suggesting that lncRNA 60967.1 promotes immune
   infiltration within the tumor microenvironment. Furthermore, we
   performed additional characterization of immune cell activation,
   exhaustion, and proliferation in the tumor microenvironment using IHC
   assays. Specifically, we analyzed the effects of lncRNA 60967.1 on the
   following markers: activation markers (IFN-γ and Granzyme B),
   exhaustion markers (PD-1 and TIM-3), and proliferation marker (Ki67) in
   CT26 and MC38 mouse tumor tissues. Our results showed that
   overexpression of lncRNA 60967.1 significantly increased IFN-γ
   expression in the tumor microenvironment, suggesting enhanced immune
   activation, as shown in Fig. [193]6H. In contrast, PD-1 and Ki67 levels
   were significantly reduced, suggesting a potential suppression of
   immune exhaustion and decreased cell proliferation (Fig. [194]6H).
   However, we did not observe statistically significant changes in
   Granzyme B and TIM-3 expression, suggesting that the effects of lncRNA
   60967.1 may be more selective in its modulation of immune cell
   activation. These findings indicate that lncRNA 60967.1 overexpression
   enhances key immune and signaling components, ultimately suppressing
   tumor growth and promoting a more active immune response in vivo.

Overexpression of LncRNA enhances the efficacy of anti–PD-1 therapy in a
murine CRC model

   LncRNAs have emerged as key regulators of immune evasion and
   therapeutic resistance in CRC [[195]45]. Building on the observed
   tumor-suppressive and immune-modulatory effects of lncRNA 60967.1 in
   the tumor microenvironment, we next explored whether its overexpression
   could enhance the efficacy of immune checkpoint blockades in a murine
   CRC model. Overexpression of lncRNA 60967.1 via lentiviral transduction
   led to a significant upregulation of the immune-related checkpoint
   molecule PD-L1, while the expression patterns of CD47 and CD276 were
   inconsistent across three CRC cell lines, as measured by qRT-PCR
   (Fig. [196]7A1, [197]7A2 and [198]7A3). Western blot analysis confirmed
   that the observed changes in PD-L1 were consistent at protein levels
   across three CRC cell lines (Fig. [199]7B).

Fig. 7.

   [200]Fig. 7
   [201]Open in a new tab

   LncRNA 60967.1 upregulates immune checkpoint gene expression and
   enhances anti–PD-1 efficacy in CRC mouse models. (A1-A3) qRT-PCR
   analysis of PD-L1, CD47, and CD276 mRNA levels in three CRC cell lines
   (DLD1, A1; HCT8, A2; and SW480, A3) transfected with either a GFP
   control or an lncRNA 60967.1 expression vector. (B) Western blot
   analysis of PD-L1 levels in DLD1, HCT8, and SW480 cells transfected
   with either a GFP control or a lncRNA 60967.1 expression vector. (C)
   qRT-PCR analysis of lncRNA 60967.1, PLCD4, IFNGR1, and IFNGR2
   expression in HCT8 cells transfected with an siRNA vector targeting
   lncRNA 60967.1 (HCT8-siRNA-1) or a non-targeting control vector
   (HCT8-NC). (D) qRT-PCR validation of PD-L1, CD47, and CD276 mRNA
   expression in HCT8 cells transfected with siRNA targeting lncRNA
   60967.1 (HCT8-siRNA-1) or a non-targeting control (HCT8-NC). (E)
   Western blot analysis of PD-L1 expression in HCT8 cells transfected
   with siRNA targeting lncRNA 60967.1 (HCT8-siRNA-1) or a non-targeting
   control (HCT8-NC). (F1–F2) qRT-PCR measurement of immune checkpoint
   gene mRNA levels (PD-L1, CD47, and CD276) in tumor tissues from Balb/c
   mice injected with CT26 cells (F1) or MC38 cells (F2). Cell lines were
   transfected with GFP (control) or lncRNA 60967.1 vectors before
   injection to compare checkpoint gene expression. (G) Tumors harvested
   from Balb/c mice that were subcutaneously injected with CT26 cells
   transfected with either GFP (Top) or lncRNA 60967.1 vectors (Bottom),
   followed by intraperitoneal administration of anti–PD-1 monoclonal
   antibody. (H) Tumor growth over time in Balb/c nude mice injected with
   CT26 cells transfected with either GFP or lncRNA 60967.1 vectors,
   followed by intraperitoneal administration of anti-PD-1 monoclonal
   antibody. (I) Tumor weight in Balb/c nude mice injected with CT26 cells
   transfected with either GFP or lncRNA 60967.1 vectors, followed by
   intraperitoneal administration of anti-PD-1 monoclonal antibody. NS
   indicates no significant difference (t-test); **P < 0.05, ***P < 0.01
   (t-test)

   To further dissect the functional importance of lncRNA 60967.1, we used
   three different siRNAs in HCT8 cells, ultimately selecting siRNA-1 for
   its higher specificity (Fig. [202]7C). Knockdown of lncRNA 60967.1
   resulted in reduced expression of PLCD4, IFNGR1, and PD-L1 at both the
   mRNA and protein levels (Fig. [203]7D-[204]E). These observations
   suggest that lncRNA 60967.1 may play a regulatory role in immune
   checkpoint pathways, potentially affecting tumor immune evasion
   mechanisms.

   Based on these findings, we hypothesized that lncRNA 60967.1
   overexpression could synergize with immunotherapy. To test this
   hypothesis, we first examined whether overexpression of lncRNA 60967.1
   increases PD-L1 levels in the murine CRC cell lines CT26 and MC38. qPCR
   analysis confirmed that PD-L1 expression was elevated in these cell
   lines following lncRNA 60967.1 overexpression (Fig. [205]7F1 and
   [206]7F2). Next, we established a murine CRC model by transplanting
   CT26 cells engineered to overexpress lncRNA 60967.1 into
   immunocompetent mice. When treated with anti-PD-1 monoclonal antibody
   (clone RMP1-14), the lncRNA overexpression group exhibited a
   significantly enhanced anti-tumor response with significant reductions
   in tumor volume (Fig. [207]7G), tumor weight (Fig. [208]7H), and
   proliferation rate (Fig. [209]7I), outperforming the lncRNA
   overexpression group alone. These data support the notion that
   combining lncRNA-targeted therapy with PD-1 blockade could offer a more
   effective therapeutic approach to CRC.

   Together, these results indicate that lncRNA 60967.1 may suppress CRC
   cell proliferation by modulating the tumor microenvironment and
   regulating key immune checkpoints. While additional investigations are
   needed, our findings point toward a similar inhibitory effect of lncRNA
   60967.1 on human CRC cells through its capacity to influence
   immune-related pathways, highlighting its potential as a combinatorial
   target for immunotherapy.

Discussion

   Multi-omics analysis has become a powerful tool for uncovering
   pathogenic mechanisms in human diseases, particularly in cancer
   [[210]46, [211]47]. In this study, we integrated transcriptomic,
   proteomic, and metabolomic data of CRC tissues and matched NATs
   collected from 13 patients. The analysis revealed a complex interaction
   network involving 22 lncRNAs, 14 mRNAs/proteins, and 9 metabolites
   (Fig. [212]2G). Among these, lncRNA 60967.1 emerged as a key regulator.
   Validation in CRC cell lines and an additional patient cohort confirmed
   the significant downregulation of both lncRNA 60967.1 and PLCD4.
   Furthermore, a concurrent reduction in the metabolite ATRA in CRC cells
   supports the involvement of the lncRNA-PLCD4-ATRA axis in CRC
   pathogenesis. Interestingly, lncRNA 60967.1 and PLCD4 mRNA were found
   to colocalize in the cytoplasm, and overexpression of lncRNA 60967.1
   increased both PLCD4 mRNA and protein levels, suggesting a regulatory
   role. However, the precise mechanism by which lncRNA 60967.1 modulates
   PLCD4 expression remains unclear. Given that lncRNAs can influence gene
   expression through multiple mechanisms—such as altering mRNA stability,
   modulating translation, or acting as competing endogenous RNAs—further
   investigation is warranted to elucidate the underlying mechanism
   [[213]45, [214]48–[215]50].

   PLCD4 encodes phospholipase C delta 4, a member of the
   phosphoinositide-specific phospholipase C (PLC) family, which
   hydrolyzes phosphatidylinositol 4,5-bisphosphate into inositol
   trisphosphate and diacylglycerol, key mediators of calcium signaling,
   cell proliferation, and differentiation [[216]51–[217]53]. The PLC
   family consists of 13 isoforms grouped into six subfamilies with
   distinct regulatory mechanisms and tissue-specific roles [[218]53,
   [219]54]. Although research on PLCD4 is limited, mutations in the gene
   have been observed in 2.18% of cancers, including melanoma,
   glioblastoma, skin squamous cell carcinoma, and adenocarcinoma of
   unknown primary origin [[220]55]. Leung et al. demonstrated that PLCD4
   overexpression in the MCF7 breast cancer cell line promoted cell
   proliferation by upregulating ErbB1/2 expression and activating the Erk
   signaling pathway, suggesting an oncogenic role [[221]56]. In contrast,
   our data showed that ectopic expression of PLCD4 in CRC cells
   significantly inhibited their migratory capacity and reduced tumor size
   in Balb/c nude mice xenografted with PLCD4-overexpressing CRC cells.
   Similarly, overexpression of lncRNA 60967.1 increased PLCD4 expression,
   which led to reduced tumor sizes in immunocompetent Balb/c mice
   injected with CRC cells overexpressing the lncRNA. These findings
   suggest that PLCD4 functions as a tumor suppressor in CRC but may have
   an oncogenic role in other cancers like breast cancer. The distinct
   roles of PLCD4 in various cancer types may be attributed to its
   involvement in different signaling pathways. In breast cancer, the DAG
   generated by PLCD4 overexpression can activate specific PKC isoforms,
   such as PKC-φ. This activation subsequently enhances the transcription
   of EGFR and ErbB2, leading to the activation of the Erk signaling
   pathway, which promotes a more proliferative cell growth state
   [[222]56]. In contrast, reduced PLCD4 expression in CRC likely engages
   alternative pathways that influence cell growth. For instance,
   overexpression of PLCD4 in rhabdomyosarcoma (RMS) RD cells upregulates
   cyclin B1, while knockdown of PLCD4 in A204 RMS cells reduces cyclin B1
   expression, suggesting that nuclear PLCD4 regulates cell proliferation
   via cyclin B1 [[223]57]. Similarly, an independent study demonstrated
   that PLCD4 knockdown affects cell proliferation and induces senescence
   in mesenchymal stromal stem cells [[224]58]. However, the precise
   mechanisms by which PLCD4 functions in CRC remain under intensive
   investigation.

   ATRA is a biologically active derivative of vitamin A and plays a key
   role in cellular differentiation, growth regulation, and apoptosis by
   activating retinoic acid receptors [[225]37, [226]59]. Notably, ATRA is
   the standard therapy for acute promyelocytic leukemia as a
   differentiation inducer [[227]60, [228]61]. While its effectiveness in
   treating solid tumors remains limited, several studies have
   demonstrated its anti-tumor effects in thyroid, breast, ovarian, liver,
   and lung cancers by inhibiting cancer cell proliferation and invasion
   [[229]62–[230]64]. In this study, we observed reduced ATRA levels in
   both CRC patients and cell lines. Although overexpression of lncRNA
   60967.1 increased ATRA levels, ATRA alone did not impact CRC cell
   proliferation or apoptosis. However, when combined with IFN-γ
   stimulation, ATRA significantly inhibited cell growth, suggesting a
   synergistic effect between ATRA and IFN-γ or a potential role for ATRA
   in modulating IFN signaling [[231]65–[232]67]. Additionally, ATRA
   treatment in IFN-γ-resistant CRC cells upregulated IFNGR1 expression
   (Fig. [233]4J), a key component of the IFN-γ receptor complex, further
   supporting ATRA’s role in modulating IFN signaling. A recent study
   confirmed that ATRA restores IFN-γ sensitivity by regulating bisecting
   N-glycosylation and stabilizing IFNγRa [[234]41]. Moreover,
   downregulation of IFNGR1 in murine CRC cells conferred resistance to
   anti-PD-1 therapy [[235]68].

   Recent studies indicate that ATRA also regulates both the innate and
   adaptive immune systems, playing a critical role in modulating
   intestinal immunity [[236]69, [237]70]. In CRC, microbiota-induced
   intestinal inflammation disrupts ATRA metabolism, leading to a colonic
   ATRA deficit that exacerbates carcinogenesis. Conversely, ATRA
   supplementation has been shown to inhibit cancer progression through a
   CD8 + T cell-dependent mechanism [[238]71]. Consistent with these
   findings, we observed that ectopic expression of lncRNA 60967.1
   upregulates the PLCD4/ATRA axis, resulting in increased infiltration of
   both CD8⁺ and CD4⁺ T cells and a corresponding reduction in tumor sizes
   in CRC xenograft models compared to controls. Notably,
   tumor-infiltrating CD8⁺ T cells are prolific producers of IFNγ, a
   cytokine critical for antitumor immunity [[239]72, [240]73]. Given the
   increased expression of the key IFNγ pathway gene IFNGR1 in CRC
   xenograft models (Fig. [241]6), our results suggest that lncRNA 60967.1
   suppresses tumor development by modulating tumor microenvironment
   through ATRA, enhancing immune infiltration, activating IFNγ signaling,
   and inhibiting tumor proliferation. However, the precise mechanisms by
   which lncRNA 60967.1 regulates ATRA metabolism and signaling remain to
   be elucidated and warrant further investigation.

   Overexpression of lncRNA 60967.1 significantly enhanced PD-L1
   expression in murine CRC cells and synergized with anti–PD-1 therapy to
   improve tumor control (Fig. [242]7G), suggesting that lncRNA 60967.1
   modulates immune checkpoint pathways and can serve as a promising
   target for combination immunotherapy (Fig. [243]8). Furthermore,
   knockdown experiments confirmed its regulatory role in immune-related
   genes (Fig. [244]7C and [245]D), further underscoring its potential as
   a key modulator of the tumor microenvironment. Nevertheless, the
   molecular pathways through which lncRNA 60967.1 governs PD-L1
   expression have yet to be fully defined. Whether these effects occur
   via a single, shared pathway or through multiple distinct signaling
   axes remains a critical question for future research [[246]48, [247]74,
   [248]75]. Furthermore, lncRNA delivery via viral vectors, lipid-based
   nanoparticles, or viral-like particles holds potential for improving
   therapeutic results in cancer treatment [[249]76, [250]77].

Fig. 8.

   [251]Fig. 8
   [252]Open in a new tab

   A proposed model suggests that lncRNA 60967.1 plays a regulatory role
   in modulating the PLCD4/ATRA axis and anti-PD-1 therapy, influencing
   immune responses and affecting CRC progression

   In summary, by analyzing multi-omics data including transcriptomics,
   proteomics, and metabolomics, we identified a novel regulatory network
   involving lncRNA, PLCD4, ATRA, and PD-L1 in CRC. The lncRNA
   60967.1-PLCD4-ATRA axis plays a role in regulating immune responses and
   tumor progression in CRC. Overexpression of lncRNA 60967.1 increases
   PLCD4 expression and ATRA levels, which enhances immune cell
   infiltration and promotes IFN-γ-induced apoptosis. This modulation may
   help reverse immune resistance, particularly in response to immune
   checkpoint inhibitors like anti–PD-1 therapy. The clinical relevance of
   this axis lies in its potential to improve treatment outcomes,
   especially in patients resistant to current therapies. Targeting this
   axis could serve as a therapeutic strategy to enhance the efficacy of
   immunotherapy in CRC, and further studies may reveal its utility as a
   biomarker for treatment response. However, certain limitations must be
   acknowledged, such as the relatively small sample size and the
   inability to fully elucidate the regulatory mechanisms among lncRNA
   60967.1, PLCD4, ATRA, and PD-L1 due to the complexity of the involved
   networks [[253]40, [254]41, [255]48, [256]75, [257]78]. Moreover, a
   total of 13 patients were enrolled in the study; however, tissue
   samples from two patients were insufficient for proteomic and
   metabolomic analysis. Due to limited tissue availability from two
   patients, the final multi-omics datasets include 26 transcriptomes, 22
   proteomes, and 66 metabolomes (triplicates per sample).

Electronic supplementary material

   Below is the link to the electronic supplementary material.
   [258]Supplementary Material 1^ (6.8MB, docx)
   [259]Supplementary Material 2^ (39.1KB, docx)

Author contributions

   SD, NZ, HG, FZ, LX, and YC conceived and designed the study. NZ, YC,
   LX, YW, XJ, CL, YS, FW, and YL collected samples and performed the
   experiments. NZ, YC, LX, JH, and XL conducted the data analysis. SD,
   NZ, HG, XL, YC, and FZ drafted the manuscript. All authors reviewed the
   manuscript.

Funding

   SD is supported by grants from the 2022 Zhejiang Provincial Key R&D
   Program (No. 2022C03032), the 2023 Provincial-Ministerial Joint Key
   Project (No. WKJ-ZJ-2313), and the National Natural Science Foundation
   of China (No. 82072628). HG is supported by the Noncommunicable Chronic
   Diseases–National Science and Technology Major Project (No.
   2024ZD0520100) and the CASHIPS Seed Grant. FZ received the Hundred
   Talents Program Award from the Chinese Academy of Sciences. LX is
   supported by the National Natural Science Foundation of China (No.
   32200593) and the Natural Science Foundation of Zhejiang Province (No.
   LQ23H160040).

Data availability

   No datasets were generated or analysed during the current study.

Declarations

Ethics approval and consent to participate

   The study received ethical approval from Sir Run Run Shaw Hospital,
   Hangzhou, China. All patients provided written informed consent. The
   Institutional Animal Care and Use Committee (IACUC) of the respective
   institutions approved all experiments.

Consent for publication

   All the authors have read and approved the final manuscript for
   publication.

Competing interests

   The authors declare no competing interests.

Footnotes

   Publisher’s note

   Springer Nature remains neutral with regard to jurisdictional claims in
   published maps and institutional affiliations.

   Yiyi Chen, Ningning Zhao, and Lingna Xu contributed equally to the
   work.

Contributor Information

   Ningning Zhao, Email: nnzhao@cmpt.ac.cn.

   Fan Zhang, Email: fzhang@cmpt.ac.cn.

   Hongcang Gu, Email: gu_hongcang@cmpt.ac.cn.

   Sheng Dai, Email: daimd@zju.edu.cn.

References