The long non-coding RNA MEG3 plays critical roles in the pathogenesis of cholesterol gallstone

Animal modeling and sample collection

Totally, 16 C57 male mice were purchased from Beijing Vital River Laboratory Animal Technology Co., Ltd. (Beijing, China). The mice were fed with chow diets in specific pathogen free (SPF) laboratory animal room for one week. Then, the mice were randomly divided into model group (n = 8) and normal control group (n = 8). Mice were housed at 22 ± 2 °C and 60 ± 10% relative humidity in a specific pathogen-free environment, with a 12:12 h light: dark cycle. The model group was fed with lithogenic diets (containing 15% fat, 1% cholesterol, and 0.5% sodium cholate) (Jiangsu Xietong Pharmaceutical Bio-engineering Co., Ltd., Jiangsu, China) for 5 weeks. Meanwhile, the normal control group was fed with chow diets (Wang et al., 2018; Tanaka et al., 2018). During the 5 weeks, food and water were ad libitum. After an overnight fasting, but free access to water, mice were anesthetized with 4% chloral hydrate by intraperitoneal injection. The liver, gallbladder and bile were subsequently isolated, photographed, and kept at −80 °C. The experiments were conducted in accordance with the National Institutes of Health guide for the care and use of laboratory animal, and also approved by the Animal Care and Use Committee in The Second Department of Biliary Surgery, Eastern Hepatobiliary Surgery Hospital.

Gallbladder bile lipid analysis

According to the manufacturer’s instructions, the changes of total cholesterol (TC), total bile acid (TBA), total bilirubin (TBL), and direct bilirubin (DBL) in bile were detected by corresponding kits (Nanjing Jiancheng Bioengineering Institute, Nanjing, China). Besides, the ratios of TC, phospholipids (PL), and TBA in the model and normal groups were calculated according to the previous reported methods (Carey, 1978). The critical Carey tables were used to calculate the cholesterol saturation index (CSI) to evaluate the cholelithiasis (Carey, 1978).

RNA-seq library construction and data preprocessing

Using Trizol reagent (Invitrogen, Shanghai, china), total RNA was extracted from four liver tissues from the model group and three liver tissues from the normal control group following the manufacturer’s manual. Then, the integrity and purity of the total RNA were detected by Agarose Gel Electrophoresis and spectrophotometer (Merinton, Beijing, China), respectively. After the RNA-seq library was established using the NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (New England Biolabs, Beverly, MA, USA), library purification, library detection, library quantitation, and cBOT automatic clusters successively were conducted. Furthermore, sequencing was performed using the TruSeq SBS kit v4-HS (Illumina, San Diego, CA, USA).

Quality assessment of the sequencing data was performed by FastQC software (Brown, Pirrung & McCue, 2017) (version 0.10.1, https://github.com/pnnl/fqc). Using Cutadapt software (Chen et al., 2014) (version 1.9.1, https://pypi.org/project/cutadapt/), the adapter sequences, the bases with mass value less than 20 or containing N at the 5′ or 3′ ends, and the reads with length less than 75 bp were filtered out. Subsequently, the clean data was compared with the reference genome using Hisat2 software (Keel & Snelling, 2018) (version 2.0.1, http://daehwankimlab.github.io/hisat2/, default parameters).

Differential expression analysis and enrichment analysis

Using DESeq2 package (Love, Huber & Anders, 2014) (http://www.bioconductor.org/packages/release/bioc/html/DESeq.html) in R, differential analysis between model group and normal control group was carried out. The mRNAs with the adjusted p-value <0.05 and —log₂ fold change (FC)— >1 were selected as the differentially expressed mRNAs (DE-mRNAs). The lncRNAs with p-value <0.05 and —log₂ FC—>1 were taken as the differentially expressed lncRNAs (DE-lncRNAs). Using pheatmap package (Zhang et al., 2016) (https://cran.r-project.org/web/packages/pheatmap/) in R, hierarchical clustering analysis was performed and clustering heatmap was drew.

Based on clusterProfiler package (Yu et al., 2012) (http://bioconductor.org/packages/release/bioc/html/clusterProfiler.html) in R, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses for the DE-mRNAs were implemented. The threshold for selecting the significant results was the p-value <0.05.

Protein–protein interaction (PPI) network analysis

Under the threshold of PPI score (medium confidence) >0.4, PPI network analysis for the DE-mRNAs was conducted using STRING database (Szklarczyk et al., 2017) (http://string-db.org). Combined with Cytoscape software (Kohl, Wiese & Warscheid, 2011) (http://www.cytoscape.org/), the PPI network was constructed. The CytoNCA plug-in (Tang et al., 2015) (parameter: without weight; http://apps.cytoscape.org/apps/cytonca) in Cytoscape software was used to analysis the topology properties of network nodes. The hub nodes (Hsing, Byler & Cherkasov, 2008) were selected according to Degree Centrality (DC), Betweenness centrality (BC), and Closeness centrality (CC) of network nodes.

Co-expression analysis and prediction of the genes targeted by miRNAs

Pearson correlation coefficients (Schober, Boer & Schwarte, 2018) of the DE-lncRNAs and the DE-mRNAs were calculated. The r > 0.95 and p-value <0.05 were utilized for screening the significant results. Using miRanda database (Joshi et al., 2016) (http://www.microrna.org), the miRNAs targeting the DE-lncRNAs and the DE-mRNAs were predicted. Under the thresholds of score >1200 and energy <-150, the significant miRNA-lncRNA pairs and miRNA-mRNA pairs were selected.

Competing endogenous RNA (ceRNA) network analysis and selection of key lncRNAs

Combined with the lncRNA-mRNA co-expression pairs, the miRNA-lncRNA regulatory pairs, and the miRNA-mRNA regulatory pairs, the mRNA-miRNA-lncRNA regulatory relationships were obtained. Using Cytoscape software (Kohl, Wiese & Warscheid, 2011), the ceRNA regulatory network was visualized.

According to the degrees of the lncRNAs in ceRNA regulatory network, the top 8 up-regulated lncRNAs and the top 8 down-regulated lncRNAs separately were selected as the key lncRNAs. Combined the lncRNA-mRNA co-expression pairs, the mRNAs co-expressed with the key lncRNAs were considered as the potential target genes of the key lncRNAs. To obtain the underlying functions of the key lncRNAs, enrichment analysis for these potential target genes was conducted using clusterProfiler package (Yu et al., 2012).

qRT-PCR analysis

After total RNA was extracted and its integrity and purity were detected, reverse transcription-PCR was performed to synthesize the first-strand cDNA. Then, BeyoFast™ SYBR Green qPCR Mix (2X) (Beyotime, Shanghai, China) and was used for qRT-PCR experiments. The amplification system included: 10 uL BeyoFast™ SYBR Green qPCRMix (2X), 1 uL cDNA template, 0.5 uL forward primer (3 uM), 0.5 uL reverse primer (3 uM), and 8ul RNase-free water. The reaction condition was: 95 °C for 2 min; 95 °C for 15 s, 60 °C for 15 s, 40 cycles; 60 to 95 °C melting curve. Using 2^−ΔΔct method, the expression of the targeted genes in relative to GADPH was calculated. The primers used in qRT-PCR experiments were listed in Table 1. Each PCR reaction had three repeats.

Animal modeling and gallbladder bile lipid analysis

In the normal control group, the bile in the gall bladder of the mice was transparent and yellow. After five weeks of lithogenic diets, the gall bladder of the mice in the model group was larger than that in the normal group, and the bile was cloudy and darker (Fig. 1A).The results of Fig. 1B showed that TC, TC/TBA, TC/PL and CSI were increased in the Model group. The TBA was decreased in the Model group (Fig. 1B). These results suggested that the mouse model of CG was successfully established.

Identification of DE-lncRNAs and DE-mRNAs

There were 181 DE-mRNAs (including 104 up-regulated mRNAs and 77 down-regulated mRNAs) and 33 DE-lncRNAs (including 17 up-regulated lncRNAs and 16 down-regulated lncRNAs) between model and normal groups. The clustering heatmaps for the DE-lncRNAs and the DE-mRNAs are shown in Fig. 2.

For the up-regulated mRNAs, 419 GO_biological process (BP) terms, 86 GO_cellular component (CC) terms, and 134 GO_molecular function (MF) terms, and eight pathways were enriched (Fig. 3A). For the down-regulated mRNAs, 229 GO_BP terms, 44 GO_CC terms, and 54 GO_MF terms, and seven pathways were enriched (Fig. 3B).

PPI network analysis

After the PPI pairs for the DE-mRNAs were predicted, PPI network (involving 101 nodes and 116 edges) was constructed (Fig. 4). According to DC, BC, and CC, protein tyrosine 166 phosphatase receptor type C (PTPRC), lysine demethylase 4A (KDM4A), and pectrin alpha, non- 167 erythrocytic 1 (SPTAN1) all were among the top 15 network nodes and were taken as the hub nodes (Table 2).

Co-expression analysis and prediction of the genes targeted by miRNAs

A total of 173 lncRNA-mRNA co-expression pairs were obtained, involving 24 lncRNAs and 96 mRNAs. For each lncRNA implicated in the co-expression pairs, enrichment analysis was conducted for its co-expressed mRNAs. Finally, 457 GO_BP terms, 80 GO_CC terms, and 137 GO_MF terms, and 11 pathways were enriched (Fig. 5). Based on miRanda database, 9320 miRNA-lncRNA pairs (involving 1754 miRNAs and 19 lncRNAs) and 86 miRNA-mRNA pairs (involving 49 miRNAs and 10 mRNAs) were predicted.

CeRNA network analysis and selection of key lncRNAs

Combined with the mRNA-miRNA-lncRNA regulatory relationships, the ceRNA regulatory network (involving 24 up-regulated mRNAs, 53 down-regulated mRNAs, 11 up-regulated lncRNAs, 11 down-regulated lncRNAs, and 47 miRNAs) was built (Fig. 6). There were 42 miRNA-mRNA regulatory pairs, 127 miRNA-lncRNA regulatory pairs, and 115 lncRNA-mRNA co-expression pairs in the ceRNA regulatory network.

Based on the degrees of the lncRNAs in the regulatory network, the top eight up-regulated lncRNAs (RIKEN cDNA 4933407K13 gene, 4933407K13Rik; RIKEN cDNA 4833418N02 gene, 4833418N02Rik; predicted gene 8378, Gm8378; RIKEN cDNA F730311O21 gene, F730311O21Rik; RIKEN cDNA A530020G20 gene, A530020G20Rik; Opa interacting protein 5, opposite strand 1, 1700020I14Rik; RAB10, member RAS oncogene family, opposite strand, Rab10os; and predicted gene, 16,973, Gm16973) and the top eight down-regulated lncRNAs (predicted gene 15,270, Gm15270; maternally expressed 3, MEG3; RIKEN cDNA C730036E19 gene, C730036E19Rik; predicted gene 16,576, Gm16576; predicted gene 27,216, Gm27216; predicted gene 12,655, Gm12655; predicted gene 11,695, Gm11695; and predicted gene 6135, Gm6135) were screened out as the key lncRNAs. In the regulatory network, the miR-107-5p/miR-149-3p/miR-346-3p—MEG3 regulatory pairs and MEG3—PABPC4/CEP131/NUMB1 co-expression pairs existed. To predict the potential functions of these key lncRNAs, enrichment analysis for their co-expressed mRNAs was performed. Moreover, the enrichment results for four up-regulated lncRNAs (4833418N02Rik, Gm8378, 1700020I14Rik, and Gm16973) and three down-regulated lncRNAs (Gm16576, Gm27216, and Gm12655) are presented in Fig. 7.

qRT-PCR analysis

Based on qRT-PCR experiments, the expression levels of key genes differentially expressed between model and normal groups were examined. As shown in Fig. 8, the expression of KDM4A was increased, and the expressions of MEG3, PABPC4, CEP131, and NUMB1 were decreased in the model group compared with the normal group. The results of qRT-PCR analysis further supported the results of differential expression analysis.