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

Background Cholesterol gallstone (CG) is the most common gallstone disease, which is induced by biliary cholesterol supersaturation. The purpose of this study is to investigate the pathogenesis of CG. Methods Sixteen mice were equally and randomly divided into model group and normal control group. The model group was fed with lithogenic diets to induce CG, and then gallbladder bile lipid analysis was performed. After RNA-seq library was constructed, differentially expressed mRNAs (DE-mRNAs) and differentially expressed lncRNAs (DE-lncRNAs) between model group and normal control group were analyzed by DESeq2 package. Using the cluster Profiler package, enrichment analysis for the DE-mRNAs was carried out. Based on Cytoscape software, the protein-protein interaction (PPI) network and competing endogenous RNA (ceRNA) network were built. Using quantitative real-time reverse transcription-PCR (qRT-PCR) analysis, the key RNAs were validated. Results The mouse model of CG was suc cessfully established, and then 181 DE-mRNAs and 33 DE-lncRNAs between model and normal groups were obtained. Moreover, KDM4A was selected as a hub node in the PPI network, and lncRNA MEG3 was considered as a key lncRNA in the regulatory network. Additionally, the miR-107-5p/miR-149-3p/miR-346-3-MEG3 regulatory pairs and MEG3-PABPC4/CEP131/NUMB1 co-expression pairs existed in the regulatory network. The qRT-PCR analysis showed that KDM4A expression was increased, and the expressions of MEG3, PABPC4, CEP131, and NUMB1 were downregulated. Conclusion These RNAs might be related to the pathogenesis of CG.


INTRODUCTION
Gallstone disease is a kind of biliary tract diseases, in which cholesterol gallstone (CG) is the most frequent type (Pasternak et al., 2017). CG can be induced by dyslipidemia, overweight, insulin resistance, obesity, and the changes in cholesterol homeostasis (Chen,

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 Bioengineering 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 R Ultra TM RNA Library Prep Kit for Illumina R (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).
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.

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 TM SYBR Green qPCR Mix (2X) (Beyotime, Shanghai, China) and was used for qRT-PCR experiments. The amplification system included: 10 uL BeyoFast TM 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. Table 1 The primers used in quantitative real-time reverse transcription-PCR (qRT-PCR) experiments.

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.

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.

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.

DISCUSSION
Screening biomarkers in CG is beneficial for CG prevention and treatment, and a number of studies already reported several biomarkers that affect the development of CG. Joshi et al. investigated 4 novel susceptibility loci (SULT2A1, TM4SF4, GCKR, and CYP7A1) and confirmed one known locus (ABCG8) of CG (Joshi et al., 2016). Th'ng et al. (2018) reported that plasma miR-122, ull-length keratin-18 (flk-18) and caspase-cleaved keratin-18 (cck-18) concentrations were increased in patients with gallstones compared with those without. In the present study, after the mouse model of CG was successfully constructed, 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 were screened. The qRT-PCR experiments confirmed the increased expression of KDM4A, as well as the decreased expressions of MEG3, PABPC4, CEP131, and NUMB1. In the PPI network, KDM4A was selected as a hub node according to DC, BC, and CC. KDM4A expression is reduced during the activation of hepatic stellate cells and its knockdown induces the low expression of miR-29, which may provide potential therapeutic approaches for liver fibrosis (Kong et al., 2019). Through recruiting KDM4, SBP (S-ribonuclease binding protein) family protein (BRG1) activates β-catenin target genes and may contribute to hepatic homeostasis and liver repair (Li et al., 2019). These suggested that KDM4A might be correlated with the mechanisms of CG. After the regulatory network was built, the top eight up-regulated lncRNAs and the top eight down-regulated lncRNAs (including MEG3) were screened out as the key lncRNAs based on their degrees. MEG3 suppresses cell proliferation and promotes cell apoptosis in gallbladder cancer, and up-regulating MEG3 may be applied for inhibiting the deterioration of the tumor (Liu et al., 2016). MEG3 overexpression in mouse liver can destabilize Shp mRNA and induce cholestatic liver injury via interacting with polypyrimidine tract-binding protein 1 (PTBP1) (Zhang et al., 2017). Therefore, MEG3 might also play roles in the development of CG.
Moreover, the miR-107-5p/miR-149-3p/miR-346-3p-MEG3 regulatory pairs and MEG3-PABPC4/CEP131/NUMB1 co-expression pairs were found in the regulatory network. MiR-107 facilitates hepatic lipid accumulation, causes hyperglycemia and Figure 7 The enrichment results for up-regulated 4833418N02Rik, Gm8378, 1700020I14Rik, and Gm16973, and down-regulated Gm16576, Gm27216, and Gm12655. Triangles and circles represent up-regulated lncRNAs and down-regulated lncRNAs, respectively. The color changing from red to blue indicates that the significant p-value decreases. The bubble size represents the proportion of the genes involved in one term.
Full-size DOI: 10.7717/peerj.10803/ fig-7 damages glucose tolerance, and thus miR-107 plays important roles in hepatic lipid accumulation (Bhatia, Pattnaik & Datta, 2016;Joven et al., 2012). The miR-149 is upregulated in the HepG2 cells receiving the treatment of long-chain fatty acid (FFA) and contributes to lipogenesis in the HepG2 cells untreated with FFA, therefore, miR-149 serves as a promising target for treating non-alcoholic fatty liver disease (Xiao et al., 2016;An, Yang & An, 2017). The miR-346 expression in the peripheral blood mononuclear cells of primary biliary cirrhosis patients is down-regulated in relative to the healthy controls, which may be related to the pathogenesis of the disease (Tan et al., 2014). Zinc finger protein 664 (ZNF664) and PABPC4 variants have different correlations with the high density liptein cholesterol (HDL-C) in adolescents and adults, which may be induced by developmental changes or environmental differences (Middelberg  , 2012). The rs4660293 in PABPC4 is related to serum TC, HDL-C, low-density lipoprotein cholesterol (LDL-C) and apolipoprotein A-I (ApoAI ) levels in the Mulao and Han populations, and a gender-specific correlation is found in these populations . CEP131 overexpression promotes cell proliferation and migration in hepatocellular carcinoma through activating the phosphatidylinositol-3 kinase (PI3K)/AKT signaling pathway, therefore, CEP131 is an oncogene and a candidate prognostic marker in the disease (Liu et al., 2017). Numb in bile in liver mediates cholesterol reabsorption, and the G595D substitution of Numb damages NPC1 like intracellular cholesterol transporter 1 (NPC1L1)-associated cholesterol reabsorption in humans with low blood LDL-C (Wei et al., 2014). These indicated that the miR-107-5p/miR-149-3p/miR-346-3p-MEG3 and MEG3-PABPC4/CEP131/NUMB1 regulatory axises might be involved in the pathogenesis of CG.
In conclusion, 181 DE-mRNAs and 33 DE-lncRNAs between model and normal groups were identified. Besides, KDM4A was implicated in the mechanisms of CG. Furthermore, the miR-107-5p/miR-149-3p/miR-346-3p-MEG3 and MEG3-PABPC4/CEP131/NUMB1 regulatory axises might play roles in the development and progression of CG. There are some limitations in the present study. First, because the experimental budget is limited, only 3/16 control animals and 4/16 model animals used for RNA-Seq. Second, the immunohistochemical analysis in liver section of normal and control mice after week 1, 2, 3, 4, 5 was not performed. These will be the part of our future research work. Nevertheless, these results still need to be validated by more rigorous experiments.