All gallstone gallbladder tissue samples and normal gallbladder tissue samples (which were surgically removed due to liver hemangioma) were procured from the Second Affiliated Hospital of Kunming Medical University (located in Kunming, Yunnan province, China) between March 2021 and July 2021. The criteria employed for participant selection were delineated as follows: (i) individuals within the age range of 20–55 years, devoid of any history of antibiotic usage within the preceding 3 months; (ii) individuals without a medical history of viral hepatitis, malignant tumors, adenomyomatosis, chronic cholecystitis, underlying diseases, and metabolic disorders; (iii) verification of the absence of evident inflammation through pathological confirmation. Individuals meeting these criteria were randomly chosen for inclusion in the study. Participants in the control group were definitively diagnosed with liver IV-V segments, lacked any gallbladder lesions, but the gallbladder was inevitably removed due to the surgical approach, rendering it an ideal control group. Gallbladder tissue samples were procured from fresh surgical gallbladder specimens, frozen in liquid nitrogen, and subsequently stored at −80°C. High-throughput sequencing was used to analyze four paired samples, with an additional 10 paired samples collected for RT-qPCR validation analysis.

Approval for this study was granted by the ethics committee of The Second Affiliated Hospital of Kunming Medical University (Protocol Number: PJ-2021-79). The participants provided written informed consent to participate in this study. Additionally, written informed consent was obtained from the participants for the publication of any potentially identifiable images or data included in this article. Our research protocol adheres to the principles outlined in the Helsinki Declaration.

RNA extraction, library construction, and sequencing procedures adhered to previously reported methodologies (Cui et al., 2022). Total RNA isolation and purification were conducted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) in accordance with the instructions of the manufacturer. The quantity and purity of each RNA sample were assessed using NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA). RNA integrity was assessed using Bioanalyzer 2,100 (Agilent, CA, USA) with a RIN number >7.0, and confirmation was obtained through electrophoresis with denaturing agarose gel. Approximately 2 μg of total RNA underwent ribosomal RNA removal as per the Epicentre Ribo-Zero Gold Kit manuscript (Illumina, San Diego, USA). Following purification, the ribo-minus RNA underwent fragmentation into small pieces using the Magnesium RNA Fragmentation Module (NEB, cat.e6150, USA) at 94°C for 5–7 min. Subsequently, the cleaved RNA fragments were reverse-transcribed to generate cDNA using SuperScript™ II Reverse Transcriptase (Invitrogen, cat. 1896649, USA). The cDNA was then employed to synthesize U-labeled second-stranded DNAs with E. coli DNA polymerase I (NEB, cat.m0209, USA), RNase H (NEB, cat.m0297, USA), and dUTP Solution (Thermo Fisher, cat.R0133, USA). An A-base was introduced to the blunt ends of each strand to facilitate ligation to indexed adapters. The above process is completed by the cooperation of T4 DNA polymerase, T4 polynucleotide kinase and Taq DNA polymerase. Each adapter featured a T-base overhang for ligating to the A-tailed fragmented DNA. Ligating single- or dual-index adapters to the fragments was followed by size selection using AMPureXP beads. After treatment with the heat-labile UDG enzyme (NEB, cat.m0280, USA) for the U-labeled second-stranded DNAs, the ligated products underwent PCR amplification under the following conditions: initial denaturation at 95°C for 3 min; eight cycles of denaturation at 98°C for 15 s, annealing at 60°C for 15 s, and extension at 72°C for 30 s; and a final extension at 72°C for 5 min. The average insert size for the final cDNA library was 300 ± 50 bp. Ultimately, 2 × 150 bp paired-end sequencing (PE150) was performed on an Illumina Novaseq™ 6000 (LC-Bio Technology CO., Ltd., Hangzhou, China) in accordance with the recommended protocol given by the vendor.

The removal of reads containing adaptor contamination, low-quality bases, and undetermined bases was performed using fastp (Chen et al., 2018). Subsequent verification of sequence quality was conducted using fastp. Bowtie2 and Tophat2 were employed to map reads to the Homo sapiens GRCh38 genome (Langmead and Salzberg, 2012; Kim et al., 2013). The mapped reads for each sample underwent assembly using StringTie (Gupta et al., 2016). Subsequently, the transcriptomes from all samples were merged to reconstruct a comprehensive transcriptome, achieved through gffcompare (https://github.com/gpertea/gffcompare/). After generating the final transcriptome, StringTie was used to estimate the expression levels of all transcripts. Transcripts were annotated with known mRNAs, known lncRNAs, and transcripts shorter than 200 nt were discarded. CPC and CNCI were then used to predict transcripts with coding potential (Kong et al., 2007; Sun et al., 2013). Transcripts with CPC scores < −1 and CNCI scores <0 were removed. The remaining transcripts with class codes (I, j, o, u, x) were considered as lncRNAs. StringTie was used to determine the expression levels of mRNAs and lncRNAs by calculating Fragments Per Kilobase Million (FPKM = [total_exon_fragments/mapped_reads(millions) × exon_length(kB)]). Differentially expressed mRNAs (DEmRNAs) and lncRNAs (DElncRNAs) between gallstone and control samples were selected based on |log2 (fold change)| >1 and a p-value <0.05 using the “limma” R package. The location of DElncRNAs was visualized using the R package “circlize".

The investigation into the function of DElncRNAs involved the use of the Multi Experiment Matrix online tool (http://biit.cs.ut.ee/mem/index.cgi). This tool was used to predict the target genes of DElncRNAs by selecting the dataset “A-AFFY-44: Affymetrix GeneChip Human Genome U133 Plus 2.0 [HG-U133_PLUS_2]" encompassing 2,811 datasets. Subsequently, the “ClusterProfiler” R package was applied to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses for both DElncRNAs and differentially expressed messenger RNAs (DEmRNAs). Enriched GO and KEGG pathways were considered significant with an adjusted p-value threshold of <0.05. Additionally, Ingenuity Pathway Analysis (IPA) was employed to identify diseases and functions significantly associated with DElncRNAs and DEmRNAs, with a significance threshold set at p-value <0.05.

To establish the co-expression network between DElncRNAs and DEmRNAs, Pearson’s correlations were computed. The criteria for co-expression were defined as |cor| ≥ 0.97 and p-value <0.0001. Subsequently, Cytoscape software was used for the construction and visualization of the co-expression network. For the prediction of miRNAs targeting DEmRNAs, the miRWalk, TargetScan, and miRDB databases were used. The starbase database was used to predict DElncRNAs interacting with miRNAs. Subsequently, miRNA-DEmRNA and miRNA-DElncRNA pairs were extracted based on their expression patterns to construct the ceRNA network.

For cis-regulation, the top five upregulated and top five downregulated DElncRNAs with the highest degrees in the co-expression network were selected. Following the definition of cis-regulation, we identified adjacent genes located within 300 kb upstream or downstream of the aforementioned 10 lncRNAs (Font-Cunill et al., 2018). Subsequently, these adjacent genes were intersected with DEmRNAs in the co-expression network. The resulting overlapped DEmRNAs were designated as cis-regulated genes. Regarding trans-regulation, the focus was on understanding how lncRNAs function through transcription factors (TFs). The TFcheckpoint database (http://tfcheckpoint.org) provided a list of 3,479 TFs, which was then intersected with DEmRNAs. Subsequently, the TRRUST database (https://www.grnpedia.org/trrust/) was used to predict the targets of the overlapped TFs. Finally, the lncRNA-TF-mRNA network was constructed and visualized using Cytoscape software.

Total RNA was extracted using Nuclezol LS RNA Isolation Reagent (ABP Biosciences Inc, China) following the instructions provided by the manufacturer. The concentration of RNA and the ratios of OD260/280 and OD260/230 were determined using a Nanophotometer N50 (OD260/280: 1.8-2.2, RNA concentrations ˃ 100 ng/ul). After determining the concentration and purity of the RNA, qualified samples underwent reverse transcription using the SureScript-First-strand-cDNA-synthesis-kit (GeneCopoeia, USA). Subsequently, qPCR was conducted on a CFX96 Real-time PCR System (Bio-Rad, USA) using BlazeTaq™ SYBR^® Green qPCR Mix 2.0 (GeneCopoeia, USA). The thermal cycling conditions comprised 40 cycles at 95°C for 30 s, 95°C for 10 s, 60°C for 20 s, and 72°C for 30 s. The 2^-△△Ct method was used to calculate gene expressions (Schmittgen and Livak, 2008).

All data are presented as means ± standard deviation. Comparisons between two groups were conducted using Student’s t-test. A p-value <0.05 was considered statistically significant.

A total of 304 DElncRNAs, comprising 103 upregulated and 201 downregulated lncRNAs in gallstone samples compared to control samples were identified (Figure 1A; Supplementary Table S1). The Circos image illustrated the widespread distribution of these DElncRNAs across all chromosomes (Figure 1B). A comprehensive list of 926 target genes was obtained (Supplementary Table S2), and they exhibited significant enrichment in 60 biological processes (BP), 37 cellular components (CC), and 23 molecular functions (MF) (Supplementary Table S3). The top 30 GO terms are depicted in Figure 1C, highlighting the involvement of potential target genes of DElncRNAs in various cellular functions. Furthermore, the potential target genes of DElncRNAs were notably enriched in KEGG pathways related to tyrosine metabolism, metabolism of xenobiotics by cytochrome P450, and drug metabolism-cytochrome P450 (Figure 1D). Additionally, Ingenuity Pathway Analysis (IPA) indicated that DElncRNAs may play roles in cell cycle regulation, cardiovascular system development and function, cell death and survival, and cellular development, among other processes (Figure 1E).

A total of 934 DEmRNAs, comprising of 396 upregulated and 538 downregulated mRNAs in gallstone samples were identified (Figure 2A, Supplementary Table S4). The DEmRNAs exhibited significant enrichment in 723 biological processes (BPs), 48 CCs, 96 MFs, and 26 KEGG pathways (Supplementary Table S5, Supplementary Table S6). These DEmRNAs were primarily associated with metabolic processes and extracellular matrix organization, including hormone metabolic processes, terpenoid metabolic processes, and diterpenoid metabolic processes (Figure 2B). Consistent with the GO enrichment results, the DEmRNAs were predominantly enriched in metabolism and extracellular matrix-related pathways, such as retinol metabolism, porphyrin and chlorophyll metabolism, and ECM-receptor interaction (Figure 2C). Moreover, the DEmRNAs were closely related to the metabolic effects of cytochrome P450 on exogenous substances, biosynthesis of steroid hormones, bile secretion, metabolism and protein digestion and absorption, biosynthesis of cofactors, and pentose and glucuronic acid conversion metabolic pathways, among others (Figure 2D).

Through Pearson’s correlation analysis, 597 DEmRNAs were found to be strongly correlated with 194 DElncRNAs. Subsequently, a co-expression network was constructed and visualized using Cytoscape software (Figure 3). After excluding lncRNAs with no neighboring mRNAs, the top 10 lncRNAs (5 upregulated and five downregulated) with correlations in the lncRNA-mRNA network were considered as target lncRNAs. These included RP11-399H11.2, VCAN-AS1, LINC01564, RP11-499O7.7, LINC02306, ROR1-AS1, NAV2-AS2, LINC01725, AC005160.3, and RP11-708L7.6. A total of 51 mRNAs were found to be cis-regulated with 10 target lncRNAs, and the number of adjacent mRNAs for each lncRNA varied (e.g., RP11-399H11.2 had 21 adjacent mRNAs, while VCAN-AS1 had only two adjacent mRNAs) (Figure 4A). After overlapping adjacent mRNAs with mRNAs in the co-expression network, a cis-regulation network composed of 32 nodes and 34 edges (Figure 4B) was constructed using Cytoscape software. For example, the downregulated BCO2 and KLHL31 were cis-regulated by LINC02306.

Simultaneously, the trans-regulatory function of lncRNAs was investigated. Initially, 50 TFs were obtained by overlapping DEmRNAs with TFs in the Tfcheckpoint database. The corresponding TFs (DEmRNAs)-lncRNAs pairs were then extracted from the lncRNA-mRNA co-expression network (Figure 5). Subsequently, using the TRRUST database, 22 TFs out of 55 were predicted to regulate the transcription of 150 targets, with 13 targets being in the co-expression network. Finally, using Cytoscape software, the TF-mRNA and lncRNA-mRNA co-expression networks were integrated to construct the lncRNA-mRNA-TF regulatory network, comprising 927 nodes and 1,927 edges (Figure 6).

A total of 245 miRNAs were identified to target 289 DEmRNAs based on information from the miRWalk, TargetScan, and miRDB databases. Additionally, 220 miRNAs interacting with 24 DElncRNAs were predicted through the starbase database. Considering the expression patterns of mRNA-miRNA and miRNA-lncRNA pairs, a comprehensive network was constructed, consisting of 24 DElncRNAs, 201 DEmRNAs, and 120 predicted miRNAs, to form the ceRNA network (Figure 7). The 24 DElncRNAs involved in the ceRNA network encompassed eight upregulated DElncRNAs and 16 downregulated DElncRNAs.

Furthermore, we conducted RT-qPCR to validate the findings obtained from the sequencing results. The top two upregulated and downregulated DEmRNAs and DElncRNAs between gallstone and control samples were selected for validation. The top two upregulated lncRNAs were AC004692.4 and HECW1-IT1, while the top two downregulated lncRNAs were RP1-27K12.2 and LINC01564. Similarly, the top two upregulated mRNAs were SFRP4 and COMP, and the top two downregulated DEmRNAs were SLC26A3 and GSTA2. A t-test was used to compare the expression differences of these four lncRNAs and four mRNAs between gallstone and control samples. The results indicated that AC004692.4, HECW1-IT1, SFRP4, and COMP were significantly higher expressed in gallstone samples, whereas LINC01564, SLC26A3, RP1-27K12.2, and GSTA2 were markedly lower expressed in gallstone samples compared to control samples (Figure 8). These findings were consistent with the results obtained from the sequencing, reinforcing the reliability of the sequencing data.