Identification of Long Noncoding RNAs Expression Profiles Between Gallstone and Gallbladder Cancer Using Next-Generation Sequencing Analysis

Introduction

Gallbladder cancer (GBC) is the most common malignant tumor of the biliary system (accounting for 80–95%).¹ It is the sixth most common malignant tumor of the digestive tract and occurs three to four times as frequently in women as in men.² Gallbladder cancer is relatively high in Latin America and Asia.³ GBC is not sensitive to radiotherapy and chemotherapy.^1,4 There is currently no recognized chemotherapy regimen that can effectively treat gallbladder cancer.⁵ Once the tumor metastasizes or recurs after surgery, doctors are helpless and cannot effectively achieve effective treatment through comprehensive treatment. Gallbladder cancer patients have an extremely poor prognosis.⁶ Therefore, the patients would benefit from an early and accurate diagnosis of gallbladder cancer.

Due to the lack of specific clinical manifestations of gallbladder cancer, and no other highly specific auxiliary examination methods, most patients are already at the advanced stage of gallbladder cancer due to tumor progression when they are first diagnosed, and they lose the opportunity for radical surgery.⁷ Mostly, GBC is diagnosed incidentally in patients undergoing cholecystectomy for the treatment of cholecystitis or cholelithiasis.⁸ Epidemiological and pathological evidence confirms that gallstone (GS) patients are more susceptible to GBC than healthy individuals. Based on this, exploring the molecular signatures between GBC and GS could help to understand the molecular mechanism of the progression of GS to GBC, and to find highly specific biomarkers and effective drug treatment targets for GBC, which is beneficial to the development of therapeutic drugs for gallbladder cancer.

The occurrence and progression of tumors is a complex process involving multiple genes and multiple factors.⁹ Genetic factors have always been considered to play an important role.¹⁰ Studies in recent years have found that the function of RNA is not only to translate into protein. Many RNAs do not have the function of encoding protein, so-called non-coding RNA (ncRNA).¹¹ Among them, those with a length> 200 nt are named long non-coding RNA (lncRNA).¹² At the earliest, this type of molecule was mistakenly regarded as a by-product of RNA polymerase II transcription, a “noise” in the genome transcription process, without biological function, and did not attract widespread attention from scholars.¹³ However, accumulating evidence shows that lncRNAs play an important role in the development of organisms and diseases,¹⁴ including the occurrence and progression of malignant tumors.^15,16 At present, the related research of lncRNA is still in its infancy. For most lncRNAs discovered by sequencing, their specific biological function is still unclear. Therefore, lncRNA as an emerging field has broad research prospects and is worthy of in-depth research by researchers.

In this study, we used next-generation RNA-seq and bioinformatics analysis to identify differentially expressed lncRNAs (DElnRNAs) and differentially expressed genes (DEGs), which may drive the GS process to GBC. Then, the differential expression of these biomarkers was validated in patients. The objective is to provide new insights into the mechanisms of GBC onset and progression.

Materials and Methods

Ethics Statement

All aspects of the present study were approved by the Ethics Committee of The First Affiliated Hospital of Kunming Medical University and accordance with the Declaration of Helsinki. Six GBC patients without gallstones (GBC-nGS), nine GBC patients with gallstones (GBC-GS) and nine gallstone patients (GS) were admitted. The paratumor (normal tissue away from the tumor) tissues from GBC-nGS were considered NS. All patients provided written informed consent for the use of their surgical specimens in the present study in accordance with the Committee’s regulations. The clinical characteristics of the patients are presented in Table 1. Tumor and paratumor (normal tissue adjacent to the tumor) tissues from GBC patients with or without GS were collected and immediately stored in liquid nitrogen until RNA-seq (Figure 1).

Table 1 Clinical Characteristics of Cases in This Study

Figure 1 A Flowchart of the study process. Abbreviations: NS, normal gallbladder; GS, gallstone; GBC-nGS, gallbladder carcinoma without gallstone; GBC-GS, gallbladder carcinoma with gallstone; DEGs, differentially expressed genes; PPI, protein-protein interaction.

RNA Sequencing

The tissue was broken into small blocks in liquid nitrogen, then total RNAs were isolated using TRIzol (Invitrogen, Carlsbad, CA). RNA was quantified and assessed for integrity using the Agilent Bioanalyzer 2100 (Agilent, Santa Clara, CA). Ribosomal RNA was removed by performing 1 to 2 washings of RNA annealed to poly-T oligo beads (Invitrogen, Life Technologies, Grand Island, NY). RNA was fragmented and reverse transcribed using random hexamers (Invitrogen). Double-stranded cDNA synthesis was performed using Pol I and RNase H. Short fragments were purified with QiaQuick PCR extraction kit (Qiagen, Hilden, Germany) and resolved with elution buffer for end reparation and poly(A) addition followed by ligation with sequencing adaptors for cluster generation and sequenced on the Illumina HiSeq 3000 (Illumina, San Diego, CA).

The Identification of Differential Expressed Genes and lncRNAs

For each sample, the obtained sequencing reads were cleaned by 5’-end quality trimming and clipping of the sequencing adapters by SOPAnuke (http://soap.genomics.org.cn/). Then, we performed prealignment quality control of the cleaned sequencing reads using FastQC (v 0.11.5). Using STAR (version 2.5.0a), we aligned the RNA-Seq reads against human genome reference sequence GRCh37 under default parameters, allowing one or two base mismatches. Furthermore, read summarization of only reads spanning introns (intron spanning) was performed with HTSeq (version 0.6.0, https://htseq.readthedocs.io/). We included protein-coding and noncoding RNAs during read mapping, summarization, and subsequent analyses. Notably, we excluded the genes encoded on mitochondrial DNA and the Y chromosome from the analyses. In addition, the genes not expressed in 5% of the samples were excluded. Moreover, gene expression clustering analysis using k-mean clustering, principal component analysis, and subsequent statistical analyses were performed using R (version 3.0.3) and R-studio (version 0.98.1091).

RNA profile was normalized using counts per million. Using the edgeR package, we performed differential gene expression analysis for the paired samples. The DEGs and DElncRNAs were identified by using a t-test, first in GBC-nGS vs NS; GS vs NS; and second, GBC-GS vs GBC-nGS. DEGs and DElncRNAs (FDR filtering (fold change ≥ 2.0) and p < 0.05) were selected for further study. Hierarchical clustering was performed to identify DEGs and DElncRNAs expression patterns, using Cluster 3.0 (Stanford, CA, USA). Subsequently, we performed GO functional enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of DEGs using the R package ClusterProfiler.¹⁷

Functional Analysis

The gene ontology (GO; http://www.geneontology.org) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.kegg.jp/) analysis were used for understanding advanced functions and biological systems. The differential GO terms and KEGG pathways were identified using Fisher’s exact test with p < 0.05 and false discovery rate (FDR) <0.05.

Co-Expression Correlation Network Analysis

A co-expression correlation network was constructed for DEGs and DElncRNA across GBC-nGS and GBC-GS, to identify gene modules associated with the GS process to GBC. The Search Tool for Retrieval of Interacting Genes database (STRING version 10.5; https://string-db.org)¹⁸ and protein-protein interaction (PPI) network analysis (PPI score = 0.9) were used. The network was visualized using Cytoscape version 3.6.1 (http://www.cytoscape.org/).¹⁹

RT-qPCR Verification

For gene expression validation, cDNA was generated from total RNA using lnRcute lncRNA cDNA kit (Tiangen Biotech Co., Ltd., China). cDNA (250 ng) was used for qPCR to determine the lncRNA expression levels using lnRcute lncRNAqPCR kit (SYBR Green) (Tiangen Biotech Co., Ltd., China) with custom PCR primers (Table S1). Samples were run on QuantStudio 12K Flex Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.) and were subject to 3 min at 95°C before undergoing 40 cycles of denaturing at 95°C for 3 sec followed by annealing/extension at 60°C for 20 sec. Three replicate PCRs were performed. Relative expression of each gene was quantified using the gene’s 2 ^−ΔΔCt.

Statistical Analysis

All statistical analyses were performed using SPSS software (version 22.0; SPSS, Inc., Chicago, IL, USA). For comparison, a two-tailed Student’s t-test was performed, when appropriate. Data were presented as the mean ± standard error. All histograms were constructed using Prism for Windows (version 6.0; GraphPad Software, Inc., La Jolla, CA, USA). P<0.05 was considered to indicate a statistically significant difference.

Results

Characteristics of Patients and Sequencing Data

Clinical characteristics of GBC-nGS, GBC-GS and GS patients are presented in Table 1. The acquisition of RNA-seq data consists of an average of 8.08 G clean data per sample with 94.83% mean Q30. HISAT2 showed that 97.30% of reads were mapped into GRCh37. A total of 53,261 lncRNAs were detected, of which 839 lncRNAs were newly predicted. After quality control of sequencing data, there was no significant difference between the groups of patients. The quality control demonstrated that RNA sequencing produced high-quality data.

Identification of DEGs and DElncRNA Expression Profiles Between NS and GBC-nGS

Differential gene expression analysis using DESeq2, including 947 DEGs (517 upregulated and 432 downregulated) and 2871 DElncRNA (1559 upregulated and 1312 downregulated) (Figure 2A and D). The top 20 DEGs and DElncRNAs are shown in Table 2.

Table 2 Next-Generation Sequencing Analysis of the 20 Most Aberrantly Regulated mRNA and lncRNAs Between Gallbladder Carcinoma Without Gallstone (GBC-nGS) and Normal Gallbladder Tissue (NS)

Figure 2 Differences in the lncRNA and mRNA expression profiles between GBC-nGS and NS tissues. (A) Volcano plot of the differential expression analysis of DEGs (fold change ≥ 2.0 and p < 0.05). The top 15 items of DEGs of GBC-nGS were identified by (B) GO function enrichment and (C) KEGG enrichment. (D) Volcano plot of the differential expression analysis of DElncRNAs (fold change ≥ 2.0 and p < 0.05). The top 15 items of DElncRNAs of GBC-nGS were identified by (E) GO function enrichment and (F) KEGG enrichment. Abbreviations: GBC-nGS, gallbladder cancer tissues without gallstone; NS, normal gallbladder; DEG, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes pathway.

For mRNA, GO analysis suggested that 256 GO items were functionally enriched involved in complement activation, antigen binding, nuclear nucleosome, cell adhesion, angiogenesis, and extracellular matrix organization (Figure 2B). The 16 pathways (p<0.05) associated with DEGs were also enriched in Systemic lupus erythematosus, alcoholism, Transcriptional misregulation in cancers, Protein digestion and absorption, and ECM-receptor interaction signaling pathway (Figure 2C).

For lncRNAs, 1522 target genes were predicted by trans- and cis-regulation, of which 191 showed significant differences in expression. GO enrichment analyses demonstrated that 41 GO items were distinctly related to focal adhesion, protein transport, ribosome, cell-cell adherens junction, protein kinase binding et al (Figure 2E). There are 29 KEGG items with p-value<0.05, including microRNA in cancer, Ribosime, ECM-receptor interaction, Glycolysis signaling pathway, and Steroid hormone biosynthesis (Figure 2F).

DEGs and DElncRNAs (nGS vs GBC-nGS) were enriched in Transcriptional misregulation in cancers, Protein digestion and absorption, ECM-receptor interaction and microRNA in cancer, which contributed to GBC progress.

Identification of DEGs and DElncRNA Expression Profiles Between NS and GS

Differential gene expression analysis using DESeq2, including 372 DEGs (150 upregulated and 252 downregulated) and 2823 DElncRNA (1323 upregulated and 1500 downregulated) (Figure 3A and D). The top 20 DEGs and DElncRNAs are shown in Table 3. Compared to DEGs (nGS vs GBC-nGS), the number of DEGs (nGS vs GS) groups is significantly reduced. We hypothesized GS tissue was more similar to para-tumor tissues in transcriptome level and lower risk of cancer.

Table 3 Next-Generation Sequencing Analysis of the 20 Most Aberrantly Regulated mRNAs and lncRNAs Between Gallstone (GS) and Normal (NS)

Figure 3 Differences in the lncRNA and mRNA expression profiles between GS and NS tissues. (A) Volcano plot of the differential expression analysis of DEGs (fold change ≥ 2.0 and p < 0.05). The top 15 items of DEGs of GBC-GS were identified by (B) GO function enrichment and (C) KEGG enrichment. (D) Volcano plot of the differential expression analysis of DElncRNAs (fold change ≥ 2.0 and p < 0.05). The top 15 items of DElncRNAs of GBC-nGS were identified by (E) GO function enrichment and (F) KEGG enrichment. Abbreviations: NS, normal gallbladder; GS, gallstone; DEG, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes pathway.

GO analysis (p<0.05) was performed to DEGs enrichment in biological processes, cellular components and molecular functions. The results suggested that 239 items were functionally enriched and involved in the apical plasma membrane, serine-type endopeptidase activity, heparin binding, ion transport, extracellular matrix organization (Figure 3B). The 96 pathways (p<0.05) associated with DEGs (nGS vs GS) were closely associated with metabolism, including cell adhesion molecules, Protein digestion and absorption, bile secretion, mineral absorption, Drug metabolism - cytochrome P450 (Figure 3C), which contributes to the pathogenesis of gallstones (GS).

For DElncRNAs, 195 target genes were predicted by trans- and cis-regulation, and showed significant differences in expression. GO enrichment analyses demonstrated 34 GO items were distinctly related to focal adhesion, protein kinase binding, cell-cell adherens junction, ribosome et al (Figure 3E). There are 81 KEGG items with p-value<0.05, including Tuberculosis, Purine metabolism, Ribosome, Apoptosis, Shigellosis, Glycolysis signaling pathway, VEGF signaling pathway (Figure 3F). Compared to DElncRNAs (nGS vs GBC-nGS), DElncRNAs (nGS vs GS) are less enriched in cancers, such as MicroRNAs in cancers and ECM-receptor interaction.

These results suggest that the expression profile of GS may be more similar to para-tumor tissues in transcriptome level and lower risk of cancer.

Identification of lncRNA and mRNA Expression Profiles Between GBC-GS and GBC-nGS

A total of 915 DEGs (408 upregulated and 507 downregulated) and 3383 DElncRNAs (1824 upregulated and 1559 downregulated) were identified between GBC-GS and GBC-nGS tissue (Figure 4A and D). The top 20 DEGs and DElncRNAs are shown in Table 4. Although DEGs were from GBC, the number of DEGs (GBC-GS vs GBC-nGS) remains relatively high, similar to DEGs (nGS vs GS). It suggests that GBC with and without gallstones may exhibit higher differences. Interestingly, several DEGs (GBC-GS vs GBC-nGS) overlap with DEGs (nGS vs GS), rather than DEGs (nGS vs GBC-nGS). This suggests that gallstones may act independently of GBC.

Table 4 Next-Generation Sequencing Analysis of the 20 Most Aberrantly Regulated mRNAs and lncRNAs Between Gallbladder Carcinoma Without Gallstone (GBC-GS) and Withgallstone (GBC-nGS)

Figure 4 Differences in the lncRNA and mRNA expression profiles between GBC-nGS and GBC-GS tissues. (A) Volcano plot of the differential expression analysis of DEGs (fold change ≥ 2.0 and p < 0.05). The top 15 items of DEGs of GBC-GS were identified by (B) GO function enrichment and (C) KEGG enrichment. (D) Volcano plot of the differential expression analysis of DElncRNAs (fold change ≥ 2.0 and p < 0.05). The top 15 items of DElncRNAs of GBC-nGS were identified by (E) GO function enrichment and (F) KEGG enrichment. Abbreviations: GBC-GS, gallbladder cancer tissues with gallstone; GBC-nGS, gallbladder carcinoma without gallstone; DEG, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes pathway.

For DEGs, there were 333 GO items (p<0.05), which were distinctly related to complement activation, antigen binding, serine-type endopeptidase activity and receptor-mediated endocytosis (Figure 4B). 37 KEGG items (p<0.05) were enriched (Figure 4C). The unique DElncRNAs were associated with 33 GO items, including focal adhesion, protein transport, ribosome, cell-cell adherens junction et al (Figure 4E). 31 KEGG items (p<0.05) was observed, including microRNA in cancer, Endocytosis, Spliceosome, Ribosime, Protein digestion and absorption, Glycolysis signaling pathway et al (Figure 4F).

Comparing DEGs (nGS vs GBC-nGS) and DEGs (nGS vs GS), as well as DElncRNAs, KEGG pathways (GBC-nGS vs GBC-GS) take into account the characteristics of both. Some pathways are related to cancers, like Transcriptional misregulation in cancers, Protein digestion and absorption, ECM-receptor interaction and microRNAs in cancers. Meanwhile, some pathways are related to gallstone metabolism, such as Bile secretion, Rennin secretion, and Tryptophan metabolism. This implies that there may be two aspects of differences between GBC-GS and GBC-nGS. We are more inclined to believe that the gallstone may not be associated with progression of GBC.

Co-DEGs and Co-DElncRNAs Among Three Lists

Our study aims to identify molecular signatures between GBC-GS and GBC-nGS to understand the possible mechanisms through which GS progresses to GBC. There were 29 overlapping DEGs and 552 overlapping DElncRNAs identified among three lists (Figure 5A and D).

Figure 5 Differential gene and lncRNA expression profiles of overlapping and unique signatures. Venn diagram showing the number of unique and overlapping (A) DEGs and (B) DElncRNA between four groups. The top 25 items of co-DEGs were identified by (C) GO function enrichment and (D) KEGG enrichment. The top 25 items of co-DElnRNAs were identified by (E) GO function enrichment and (F) KEGG enrichment. (G) Venn diagram showing the number of overlapping between co-DEGs and co-DElncRNAs predicted genes. (H) FXYD2, MPZL1, PAH were determined to be differentially expressed in four groups compared with NS group using next-generation sequencing analysis and validated using RT-qPCR in tissues. Abbreviations: RT-qPCR, reverse transcription-polymerase chain reaction. GBC-nGS, gallbladder cancer tissues without gallstone; GBC-GS, gallbladder cancer tissues with gallstone; NS, normal gallbladder; GS, gallstone; DEG, differentially expressed genes.

The Gene Ontology analysis demonstrated co-DEGs were enriched biological processes like Cell cycle, DNA packaging, DNA packaging, Protein binding, and Cytoskeleton (Figure 5B). KEGG pathways analysis enriched pathways like alcoholism, systemic lupus erythematosus, and homologous recombination (Figure 5C). Co-DElncRNAs were also enriched in Negative regulation of gene expression, epigenetic, rDNA heterochromatin assembly, nucleosome and rDNA heterochromatin assembly and rDNA heterochromatin assembly, whereas the DElncRNAs were mainly linked with Systemic lupus erythematosus and alcoholism pathway (Figure 5E and F). Interestingly, DEGs are enriched in immune-related pathways, while DElncRNAs are enriched in cancer-related pathways, such as pancreatic cancers, Non-small cell lung cancer, AMPK signaling pathway, etc.

Next, we mainly focused on DEGs and DElncRNAs with enlarged expression fold-change trends in GBC-GS group versus GBC-nGS group. We found three genes were co-expressed in lncRNAs and mRNAs, including FXYD2, MPZL1, PAH (Figure 5G). Then, we used qPCR to examine their levels in NS, GS, GBC-GS and GBC-nGS groups. The results showed that FXYD2 and PAH were upregulated, while MPZL1 was downregulated between GBC-GS and NS or GBC-nGS group. It is consistent with RNA-sequencing data (Figure 5H).

These results suggest that FXYD2, MPZL1, and PAH may be important mediators between immunity and gallbladder cancer. Additionally, these factors may also serve as bridges between immunity and gallstones.

Discussion

Although without the function of encoding protein, lncRNA was reported to play an important role in the progression of tumors, including gastric cancer,²⁰ hepatocellular cancer,²¹ colorectal cancer²² and pancreatic cancer.²³ However, few studies have focused on the role of lncRNAs in the progression of gallstones to gallbladder cancer. In this study, the differential expressed lncRNA and mRNA were screened in gallbladder cancer patients with gallstone and gallbladder carcinoma without gallstone. Under the premise of good quality control, we obtained a total of 29 DEGs and 500 DElncRNAs. Using the obtained DElncRNAs and DEGs, a co-expression network of DEGs was constructed. The DEGs were significantly enriched in GO function, including focal adhesion cell-cell adhesion, ribosome, spliceosome, while enriched in KEGG pathway, including Transcriptional misregulation in cancers, Protein digestion and absorption, ECM-receptor interaction, bile secretion signaling pathways.

Previous studies have used microarray chips to study the expression profiles of DEGs and DElncRNA in GBC.²⁴ Compared with matched para-cancerous tissues, the level of more than 1000 lncRNA expression changed.²⁵ Our data screened more known DElncRNAs and novel DElncRNAs. The discrepancy between the data of these studies may be due to different kinds of tumor tissues and different platforms used. The platform that was used as the Human lncRNA array V4.0 (4×180K) in Kong et al’s study,²⁴ whereas the platform for the GEO data was the Affymetrix Human Gene 2.0 ST Array in Wang et al’s study.²⁵ Our study is based on next-generation sequencing data. High-throughput ncRNA sequencing screened more new lncRNAs, which avoided the omission of lncRNA which plays an important role in GBC.

We identified 947 DEGs (NS vs GBC-nGS), 372 DEGs (NS vs GS), and 915 DEGs (GBC-nGS vs GBC-GS) respectively. The numbers reflect that gallstones may be more inclined towards para-tumor tissues and contribute less to gallbladder cancer. Gallstones are very common but gallbladder cancer is very rare. Most people with an inflamed gallbladder or gallstones do not get gallbladder cancer.²⁶ Our study indirectly supports the conclusions regarding gallstones not related to gallbladder cancer. Trans-regulation and cis-regulation are important regulatory mechanisms for lncRNAs.²⁷ We utilized this mechanism to predict target genes of DElncRNAs, together with DEGs followed by GO and KEGG function analysis. Our data indicates that lncRNA transcription undergoes significant changes during the process of GBC, consistent with previous studies.^28,29 DEGs and DElncRNAs (NS vs GBC-nGS) were enriched in cancer-related pathways, while DEGs and DElncRNAs (NS vs GS) were enriched in metabolism-related pathways. Studies reported a close association between gallstones and processes such as bile secretion and transport.³⁰ It reflects the dual spectrum characteristics of DEGs and DElncRNAs (GBC-nGS vs GBC-GS), which exhibit both tumor and stone-related profiles. This suggests that gallstones and gallbladder cancer may exist independently rather than in a promoting relationship.

Several shared signaling pathways between gallstones and gallbladder cancers are worth our attention, including ECM-receptor interaction pathways, transcription misregulation, and Protein digestion and absorption. ECM-receptor interaction pathway plays an important role in the process of tumor shedding, adhesion, degradation, movement and hyperplasia.³¹ Accumulated evidence has confirmed ECM-receptor interaction pathway is associated with several cancers, including breast cancer,³¹ kidney cancer,³² and colorectal cancer.³³ In gallbladder cancer, Kong et al showed mRNAs and lncRNAs may be involved in gallbladder cancer progression via ECM-receptor interaction pathways and the complement and coagulation cascades.²⁴ Chen et al compared the glandular epithelial between normal and tumor tissues and found that extracellular matrix (ECM)–receptor interaction plays an important role in gallbladder cancer progression.³⁴ In our study, upregulated DEGs and DElncRNAs were significantly enriched in the ECM-receptor interaction pathway, which possibly promoted the development of gallbladder cancer. Transcription misregulation arises in cancer from disease-defining genetic alterations either indirectly, via mutation of signaling factors converging on transcriptional control, or directly, via genetic alterations in gene control factors themselves.³⁵ Nepal et al’ study showed that in gallbladder cancer, transcriptional levels appear to be affected by mutations, as indicated by the enrichment of mutations for the ‘Transcriptional misregulation in cancer’ pathway.³⁶ We agree with this viewpoint. In this study, we also found that transcription misregulation is significantly upregulated in GBC-GS. However, direct or indirect effects requires further validation.

Interestingly, the pathway “Protein digestion and absorption” is presented in both cancer and gallstones. Protein digestion and absorption have been consistently found to be involved in digestive system cancers from RNA-sequencing studies.³⁷ However, it has been identified in gallbladder cancer for the first time, including the FXYD2 and the Collagen family. Collagen is the major component of the tumor microenvironment and participates in cancer fibrosis. In gallstone, the total and insoluble collagen contents were significantly high in gallstone patients. The activity levels of collagenase and protease were significantly low.³⁸ In this study, Protein digestion and absorption participated in the progress of both gallstone and gallbladder cancer. In addition, DElncRNAs (NS vs GBC-nGS) and DElncRNAs (GBC-GS vs GBC-nGS), but not DElncRNAs (NS vs GS). were enriched in the pathway “MicroRNAs in cancers”. It indicated gallbladder cancer was related to miRNAs. Previous studies have shown that lncRNAs can act as competitive endogenous RNAs (ceRNAs), affecting the function of miRNAs through sponge effects and participating in tumor progression. Our results indirectly support this view.

Finally, we analyzed DEGs and DElncRNAs related to both gallstones and tumors. We initially filtered out co-DEGs and co-DElncRNAs from all four groups. Subsequently, we overlapped the genes predicted by lncRNA cis-regulation with DEGs, and identified three genes that exhibited significant changes. For instance, FXYD2 showed a significant decrease in expression in GBC-nGS compared to NS, and the decrease was even more pronounced in GBC-GS compared to GBC-nGS. Similar patterns were observed for MPZL1 and PAH. Three genes were reported to relate to gallstone or GBC. Downregulation of FXYD2 is associated with poor prognosis in clear cell renal cell carcinoma.³⁹ Prevalence ratios (PR) of gallbladder disease are significantly higher for the Phenylalanine hydroxylase (PAH) deficiency cohort.⁴⁰ MPZL1 is highly expressed in advanced gallbladder carcinoma and promotes the aggressive behavior of human gallbladder carcinoma GBC-SD cells.⁴¹ Although the dynamic changes in these genes did not directly prove the role of gallstones in GBC, we found that these coDEGs are more closely associated with immunity. Previous studies have shown that immunity may play an important role in both GS and GBC,^42,43 potentially serving as a bridge between the two.

This study still has some limitations. (1) Due to the nature of gallbladder disease, we were unable to complete follow-up processes from gallstones to GBC in the same patient. Therefore, the role of GS in GBC cannot be confirmed. (2) Due to the low incidence of GBC, it is also not possible to obtain sufficient GBC patients with different periods of follow-up gallstones. (3) Although we have used reliable bioinformatics tools to identify some DEGs and DElncRNAs that may be related to gallbladder cancer, these have not been confirmed by experimental evidence in this study. We will further verify the lncRNA-mRNA network in the process of GBC in vitro and in vivo studies.

Conclusions

In summary, we profiled the transcriptome landscape of GBC with/without GS using RNA sequencing and identified three genes (FXYD2, MPZL1 and PAH) and Transcriptional misregulation in cancers, Protein digestion and absorption, ECM-receptor interaction signaling pathways, which may contribute to the development of gallbladder cancer. The differential expressed trend and diversity of DEGs and DElncRNAs indicated GS profile was more similar with para-tumor tissues in transcriptome level and lower risk of cancer. Further exploration is necessary from GBC patients with different periods of follow-up gallstone.