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. The results of gallbladder bile lipid analysis showed that TCH, TBA, TBL, and DBL levels in the model groups were all higher than those in the normal control group (Figure 1). These 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 Figure 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 (Figure 3A). For the down-regulated mRNAs, 229 GO_BP terms (such as lipid hydroperoxide transport (involving sterol carrier protein 2, SCP2)), 44 GO_CC terms, and 54 GO_MF terms, and seven pathways were enriched (Figure 3B).
PPI network analysis
After the PPI pairs for the DE-mRNAs were predicted, PPI network (involving 101 nodes and 116 edges) was constructed (Figure 4). According to DC, BC, and CC, protein tyrosine phosphatase receptor type C (PTPRC), lysine demethylase 4A (KDM4A), and spectrin alpha, non-erythrocytic 1 (SPTAN1) all were among the top 15 network nodes and were taken as the hub nodes (Table 1).
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 (Figure 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 (Figure 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, 16973, Gm16973) and the top eight down-regulated lncRNAs (predicted gene 15270, Gm15270; maternally expressed 3, MEG3; RIKEN cDNA C730036E19 gene, C730036E19Rik; predicted gene 16576, Gm16576; predicted gene 27216, Gm27216; predicted gene 12655, Gm12655; predicted gene 11695, 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 Figure 7.