Expression profiles and gene set enrichment analysis of the transcriptomes from the cancer tissue, white adipose tissue and paracancer tissue with colorectal cancer

Tissue specimen collection

A total of 10 CRC patients were involved in this study. An experienced clinician removed WAT specimens from the mesentery 5 cm adjacent to the cancer or PCT of each patient. The collected tissue was rinsed with saline at 4 °C to remove bacteria and blood attached to the surface and immediately stored at −80 °C before sequencing and analysis. This study was approved by the Ethics Committee of the Affiliated Hangzhou Dermatology Hospital, Zhejiang University School of Medicine, and informed consent was obtained from all participants (No. 2021KA034). All patients participating in the study have signed informed consent.

High-throughput sequencing

We used TRIzol reagent (Invitrogen, Carlsbad, CA, USA) to extract RNA from primary cancer tissue, PCT, and WAT. The RNA of all samples with the Ribo-ZeroTM Magnetic Kit (Epicentre, San Diego, CA, USA) was fragmented using First Strand Master Mix (Invitrogen, Carlsbad, CA, USA). The first strand cDNA was generated by reverse transcription using random primers, after which the second strand cDNA was synthesized, the RNA template was elimated, and dTTP was replaced by dUTP. The linked products were amplified and purified. We completed RNA sequencing library construction by using BGI (Wuhan, Hubei, China).

Prediction of subcellular lncRNA Localization

We obtained lncRNA sequences from Dr. Tom Data Visualization (BGI).The lncRNA subcellular localization was predicted using input transcripts in lncLocator (Shanghai Jiao Tong University, http://www.csbio.sjtu.edu.cn/bioinf/lncLocator/) (Cao et al., 2018). GraphPad Prism 9 software (version 8.0; GraphPad Software, San Diego, CA, USA) was used to show and analyze the result.

Gene set enrichment analysis

Gene set enrichment analysis (GSEA) was performed using the Dr. Tom system with the Broad Institute software package (http://biosys.bgi.com, BGI), based on the Kyoto Encyclopedia of Gene and Genome (KEGG) pathway database.

Protein–protein interaction network analysis

To illustrate the functional interactions of the screened mRNAs, the Dr. Tom system was used to construct and visualize a protein–protein interaction (PPI) network. The scores ranged from 0 to 1,000. Significant interactions were considered with a required combined score of >950.

Construction of the mRNA–lncRNA co-expression network

The Pearson coefficient (r) >0.9 and p < 0.05 mRNA-LncRNA pairs were used to construct the co-expression network, which was visualized using the Cytoscape software (version 3.8.0, San Francisco, California, USA).

Gene ontology and KEGG pathway analysis

We further analyzed the potential biological processes, cellular components, molecular functions, and pathways of the target genes. The Dr. Tom system included Gene Ontology (GO) and KEGG pathway enrichment analyses. Predominately enriched GO terms were defined as those with a p-value < 0.05. A heatmap was created and enriched using R studio.

Expression profiles and GSEA of CRC and PCT transcriptomes

We obtained tissue samples from CRC tissue, PCT, and WAT and performed high-throughput sequencing. In order to avoid error caused by taking samples, we took CRC from the primary cancer, PCLT and WAT were taken from the normal tissue and mesangial fat, which was 5 cm from the edge of primary cancer (Fig. 1A). We identified 19,514 and 107,009 genes as mRNAs and lncRNAs, respectively. We identified 4,959 DEGs between WAT and CRC tissue (WAT/CRC), 2,848 between PCT and WAT, and 1,509 between CRC tissue and PCT (CRC/PCT).

We then used GSEA to observe the function of genes contained in CRC tissue and PCT. We used strict inclusion criteria (i.e., normalized enrichment score (NES) >1.0, nominal p-value <5%, FDR <25%). We found that 25 of the 162 detected gene sets were enriched in CRC tissue, compared to PCT. These included the ubiquitin-mediated proteolysis, Fanconi anemia pathway, steroid biosynthesis, p53 signaling pathway, and basal transcription factors. Table S1 shows the specific information of 26 significant enrichment pathways. These pathways are related to metabolism, genetic information processing, and environmental information processing. The alteration of steroid biosynthesis regulates lipid metabolism, and the change in local steroid concentration can cause abnormal adipose tissue function, such as dyslipogenesis or abnormal production of adipokines, eventually leading to lipid metabolism disorder (Ostinelli et al., 2022; Wawrzkiewicz-Jałowiecka, Lalik & Soveral, 2021). Twenty candidate genes of interest were enriched in steroid biosynthesis, and the leading edge included 14 genes (Fig. 1B). The p53 signaling pathway influences the tumorigenesis of CRC; TP53 mutations in CRC often affect the biological function of wild-type p53, such as promoting cancer cell stemness, proliferation, invasion, and metastasis, all of which contribute to cancer progression (Giacomelli et al., 2018). There were 73 candidate genes of interest enriched in the p53 signaling pathway, and the leading-edge subsets included 29 genes (Fig. 1C). In order to avoid error casued by taking samples, we took CRC from the primary cancer, PCLT and WAT were taken from the normal tissue and mesangial fat, which was 5 cm from the edge of primary cancer.

Protein–protein interaction network construction and analysis

To identify the key genes enriched in these two pathways and study the mechanism of crosstalk between lipid metabolism and CRC tumorigenesis, a network analysis of protein–protein interaction (PPI) was conducted. Genes with an overall interaction score of 950 (0–1,000) were considered important and likely core genes. Possible core genes in steroid biosynthesis and the p53 signaling pathway are shown in Figs. 2A and 2B (Table S2). Two functional clusters were constructed for the steroid biosynthesis and p53 signaling pathway. The KEGG pathway enrichment analysis of functional clusters in Fig. 2A show that genes in this functional cluster were involved in steroid synthesis, cholesterol metabolism, and parathyroid hormone synthesis (Fig. 2C). Figure 2B shows that the clusters were involved in the p53 signaling pathway and other cancer tumorigenesis-related pathways, including apoptosis (Fig. 2D).

Analysis of differential genes between CRC tissue and PCT

Using DESeq2 R package and setting strict inclusion criteria (log2 fold change ≥ 1, FDR ≤ 0.001), we obtained DEGs between CRC tissue and PCT from the RNA-seq database. We identified 1,228 genes with altered expression levels between CRC tissue and PCT (Table S3). Figure 3A depicts the genes differentially expressed between the CRC tissue and PCT. We found that, compared with PCT, 190 differentially expressed lncRNAs (DElncRNAs) and 435 differentially expressed mRNAs (DEmRNAs) were upregulated in CRC tissue, while 54 DElncRNAs and 191 DEmRNAs were downregulated. We visualized the top 80 DElncRNAs as expressive heatmaps using normalized raw z-score (Fig. 3B). Similarly, cluster analysis was performed on the top 100 DEmRNAs (Fig. 3C).

Numerous studies have highlighted the ability of lncRNAs to bind to microRNA (miRNA) sites as competing endogenous RNAs (ceRNAs), thereby affecting and regulating the expression of mRNAs and target genes (Xu et al., 2022). Based on the ceRNA network constructed by the Dr. Tom system, we identified 476 target mRNAs between CRC tissue and PCT and performed KEGG pathway enrichment analysis. The results indicated that the PPAR, TNF, NF-kappa B, MAPK, and NOD-like receptor signaling pathways were associated with the tumorigenesis of CRC and lipid metabolism (Fig. 3D).

To identify the potential functions of the DEGs between the CRC tissue and PCT, we analyzed 626 DEmRNAs using the KEGG pathway list. The results showed that the upregulated expression of 435 in CRC tissue resulted in significant enrichment of DEmRNA in 13 of the 255 pathways.The 13 KEGG pathways belong to the following level 1 categories: organic systems, cellular processes and environmental information processing, genetic information processing, metabolism, and human diseases. Table S4 shows the detailed information of the level 2 KEGG terms, including eukaryotic ribosome biogenesis, steroid biosynthesis, the IL-17 signaling pathway, cytokine-cytokine receptor interaction, homologous recombination, and nucleocytoplasmic transport. Among these 13 pathways, steroid biosynthesis was associated with lipid metabolism (Fig. 3E and Table S3). We found that 46 pathways were downregulated in CRC tissue (p < 0.05) (Table S4). According to the level 1 KEGG categories, we identified 20 pathways in organismal systems, cellular processes and environmental information processing, genetic information processing, metabolism, and human diseases. The level 2 KEGG terms include the regulation of lipolysis in adipocytes, which is associated with lipid metabolism. We also observed that the B cell receptor signaling pathway (Fig. 3F), which is related to immunity, critically regulates the tumorigenesis of cancer (Tanaka & Baba, 2020).

Analysis of DEGs between CRC tissue and PCT and between PCT and WAT, but not between WAT and CRC tissue

To further analyze the correlation between CRC progression and lipid metabolism, we selected 974 genes (log2 fold change >1, FDR <0.05) that revealed differences between CRC tissue and PCT and between PCT and WAT but not between WAT and CRC tissue. The 974 genes, 671 mRNAs, 214 lncRNAs, and other types of RNAs are presented in Table S7. We performed an expression clustering analysis (Fig. 4A). The KEGG pathway analysis revealed that the 671 mRNAs were associated with human diseases, organismal systems, environmental information processing, metabolism, cellular processes, and genetic formation processing (Fig. 4B). The GO analysis of the 671 mRNAs indicated that the genes were enriched in cytosolic ribosomes, large ribosomal subunits, large ribosomal subunits, and extracellular exosomes (Fig. 4C). When normal tissue become cancerous, they often produce more exosomes, and as exosomes are excreted, cancer cells change their immediately surrounding and distant microenvironment (Zhang & Yu, 2019). Therefore, we analyzed the genes enriched in the extracellular exosomes by using the KEGG pathway. The result suggested that the ECM-receptor interaction is related to the regulatory network and its prognostic role in CRC (Fig. 4E).

Based on the Dr.Tom ceRNA network database, we identified 2,676 mRNAs regulated by 214 lncRNAs (Fig. 4D and Table S5). A KEGG pathway enrichment analysis was performed on the genes in the identified networks, and several pathways were predominately enriched. The results revealed that the HIF-1α signaling, arginine, and proline pathways enriched in the target genes were related to the tumorigenesis of several cancers, including CRC (Fig. 4F). The absolute values of log2 (CRC/PCT) and (WAT/PCT) were ≥ 2, and 75 lncRNAs were screened. To further clarify the differences in the expression levels of these lncRNAs among the groups, we performed a cluster analysis (Fig. 4G). According to these criteria, the expression levels of both WAT and CRC genes were higher than those of PCT. The following lncRNAs were screened: MIR222HG, DUXAP8, DUXAP9, LINC00920, MIR503HG, URS0001BF3125, DNAH17-AS1, lncRNA-ATB, HSD52, FAM222A-AS1, URS00008B328B, URS0000D584C2, URS0000D5BAAB, and URS0001BF711B. Based on several studies that found that the subcellular localization of lncRNA directly affects its function, we used lncLocator to analyze the locations of these 14 lncRNAs. Database prediction results indicated that lncRNA-ATB localized to the nucleus, MIR503HG localized to the exosomes, and the others localized to the cytoplasm or cytosol (Fig. 4H). We predicted the mRNA targeting MIR503HG based on the ceRNA network database and performed a GO analysis (Fig. 4I and Table S6).

Co-expression analysis and functional prediction

To further explore the interaction of MIR503HG with CRC progression and tumorigenesis, we determined the correlation between MIR503HG and 3,988 mRNA expression levels in the DEGs between the CRC tissue and PCT by calculating the Pearson correlation coefficient (r) (Fig. 5A). We also found that MIR503HG exhibited a strong correlation with 20 mRNAs (r > 0.9, Fig. 5B). We then performed KEGG pathway enrichment analysis to identify the genes with strong correlations for functional analysis. Genes co-expressed with MIR503HG were enriched in the PI3K-Akt, Rap1, and calcium signaling pathways (Fig. 5C).