Comprehensive multiomics analysis of the effect of ginsenoside Rb1 on hyperlipidemia

Confirmation of the successful construction of a hyperlipidemia mouse model

Compared with the control mice, the hyperlipidemic mice had significantly higher serum levels of triglyceride (TG), total cholesterol (TC) and low-density lipoprotein C (LDL-C), but significantly lower serum levels of high-density lipoprotein C (HDL-C). This indicates that our rat hyperlipidemic model was successful (Table 1 and Supplementary Figure 1A).

Microbiome sequencing characteristics and differential abundance of bacteria

10,091,546 total sequences composed of 4,159,619,833 total bases and an average length of 412.19 bp were obtained from the optimized sequences. Based on these, 3730 operational taxonomic units (OTUs) were obtained with an average size of 27.14 organisms. After flattening for cluster analysis, 40 phyla, 72 classes, 82 orders, 153 families, 335 genus and 461 species were obtained (Figure 1A). The PCA plot intuitively displayed the similarities and differences between samples. The model and control groups were obviously discrete on the two-dimensional graph (Figure 1B), and the Rb1 and model groups were relatively discrete (Figure 1C). After statistical analysis of the annotated species, we found that at the genus level, the relative abundances of Blautia and Allobaculum were greater, and Turicibacter was lower, in the model group than the control group (Figure 1D). However, when Rb1 was fed to the model group, the relative abundance profile of these bacteria was reversed. KEGG analysis showed that the level of enrichment in ether lipid metabolism, glycerolipid metabolism, glyoxylate and dicarboxylate metabolism, and type II diabetes mellitus significantly differed between the Rb1 and model groups (Figure 1E).

Figure 1. Characteristics of 16S rRNA sequencing. (A) TaxonTree revealed the distribution of each species in all samples. Each layer in the leaf node represents a classification level, from left to right is the boundary, gate, class, order, family, genus and species level, only genus annotated and with relative abundance large than 1% was shown in the figure. (B) PCA of OUT between Control and Model group. Each point in the graph represents a sample. The different colors of the points represent the grouping of the samples. (C) PCA of OUT between Rb1 and Model group. (D) Relative abundance of bacteria in genus level in each group. Blautia and Allobaculum increased in Model group and Turicibacter decreased in Model group. Notes: data was shown in mean +/- SD, * represent P<0.05, ** represent P<0.01, *** represent P<0.001. (E) KEGG analysis of abundance difference between Rb1 and model group. The result was performed by PICRUSt software with default parameters. Each color represents a group of samples, their average relative abundance is showed in X axis. The left of bar chart shows the pathways with significant differences in abundance between the two groups, and the right of chart shows P value of each pathway.

Characteristics of the lipidomics analysis and differential expression of lipids

In positive ion mode (POS), a total of 4119 peaks were detected, 1891 metabolites were identified using the interquartile range denoising method, 815 metabolites significantly differed between the control and model group, and 169 metabolites were significantly differed between the model and Rb1 group. The differences between the control and model groups, between the model and Rb1 groups were distinct in the OPLS-DA plot (Figure 2A, 2B). The original permutation test for the OPLS-DA model explained the difference between the two compared groups of samples (Figure 2C, 2D). Using a criterion of P<0.05 in the enrichment analysis, KEGG and MetaboAnalyst revealed that different metabolites in comparison between the control and model groups mainly enriched in glycerophospholipid metabolism and linoleic acid metabolism pathways (Figure 2E and Supplementary Table 1), while different metabolites in comparison between the Rb1 and model groups enriched in glycerophospholipid metabolism, linoleic acid metabolism and alpha-linolenic acid metabolism pathway (Figure 2F). Phosphatidylcholine, diacylglycerol and ceramide were the main enriched metabolites in the glycerophospholipid, linoleic acid and alpha-linolenic metabolism pathways, respectively.

Figure 2. Characteristics of lipidomics in positive ion model. (A) Score scatter plot of OPLS-DA model for group C vs M in POS model. The X axis t [1] P denotes the predicted principal component score of the first principal component, the Y axis t [1] O denotes the orthogonal principal component scores. The two groups of samples are very distinct in this model. C, control group. M, model group. POS, positive ion model. (B) Score scatter plot of OPLS-DA model for group M vs R in POS. M, model group. R, Rb1 group. POS, positive ion model (C) Permutation test of OPLS-DA model for group Control vs Model in POS. The X axis represents the replacement retention, the Y axis represents the value of R²Y or Q², the green dot represents the value of R²Y and the blue square represents the value of Q of the replacement test. The dotted lines represent the regression lines of R Y and Q, respectively. The original model can well explain the difference between the two groups of samples. (D) Permutation test of OPLS-DA model for group Rb1 vs Model in POS. (E) Pathway analysis for group Control vs Model in POS. The bubble size indicates the influencing factor in the topological analysis; the bubble color represents the P value of enrichment analysis. (F) Pathway analysis for group Rb1 ve Model in POS.

In the negative ion model (NEG), 2734 peaks and 1348 metabolites remained, 401 metabolites significantly differed between the control and model group, while 58 metabolites significantly differed between the model and Rb1 group. The same difference trend obtained with the POS model is illustrated in Supplementary Figure 2A, 2B, and the original permutation test of the OPLS-DA model well explains the difference (Supplementary Figure 2C, 2D). KEGG and MetaboAnalyst revealed that glycerophospholipid metabolism, linoleic acid metabolism, alpha-Linolenic acid metabolism, and glycosylphosphatidylinositol (GPI)-anchor biosynthesis were the main enriched pathways when comparing the control and model groups (Supplementary Figure 2E), while glycerophospholipid metabolism, linoleic acid metabolism, alpha-linolenic acid metabolism and arachidonic acid metabolism were the main enriched pathways when comparing the Rb1 and model groups (Supplementary Figure 2F). The main enriched metabolites in the NEG model were the same as in the POS model.

Characteristics of transcriptome sequencing and differential expression of transcripts

A total of 866.89 million total raw reads were obtained from the library, and 798.38 million (87.8%) clean reads with an average length of 144.67 nucleotides were screened for subsequent analysis. A total of 75.98% of clean reads mapped to the reference genome. We identified 118,009 transcripts, including 82,038 lncRNAs, 34,641 mRNAs, 418 circRNAs and 912 other RNAs. Only 8.64% of transcripts were known, and most assembled transcripts were novel. Fragments per kilobase of exon per million mapped reads (FPKM) were used to normalize the expression of transcripts. With replicate expression, a heat map was used to distinguish sample groups, and the differences between the control and model groups and between the model and Rb1 groups were distinct (Figure 3A, 3B). The average expression in the 9 samples was nearly 1.0 FPKM (Figure 3C). A total of 510 transcripts were differentially expressed between the model and control group, while 766 transcripts were differentially expressed between the Rb1 and the model group. The significant differences in transcript expression were depicted with clustering analysis (Figure 3D). GO analysis revealed that small molecule metabolic process, cellular carbohydrate metabolic process, regulation of lipid catabolic process, and L-serine metabolic process were the main enriched terms when comparing the model group with the control group (Supplementary Table 2). DNA metabolic process, DNA repair and alpha-amino acid metabolic process were the main enriched terms when comparing the Rb1 and model groups (Supplementary Table 3). The KEGG pathway analysis indicated that glycine, serine and threonine metabolism, mannose type O-glycan biosynthesis, amino sugar and nucleotide sugar metabolism, steroid biosynthesis and the PPAR signaling pathway are involved in the establishment of the high lipid model (Figure 3E). Rb1 may participate in lipid regulation by affecting lysine degradation, inositol phosphate metabolism, and glycerophospholipid metabolism (Figure 3F).

Figure 3. Characteristics of transcriptomics. (A) Heat map of intersample correlation between model and control group. The correlation was evaluated by Pearson’s correlation coefficient, different color represents differences among different samples. (B) Heat map of intersample correlation between Rb1 and model group. (C) FPKM distribution boxes of transcripts in each sample. Y axis is shown in log₁₀ scale, X axis one sample per region. (D) Clustering map of differential genes. The depth of the color represents the level of gene expression. (E) KEGG pathway analysis of differentially expressed genes in control group versus model group. Gene ratio represents enrichment factor, bubble scale represents number of different genes, depth of bubble color represents P value. (F) KEGG pathway analysis of differentially expressed genes in Rb1 group versus model group.

Comprehensive multiomics analysis reveals that Rb1 affects glycerophospholipid metabolism

We found that phosphatidylcholine levels differed between the control and model groups and between the Rb1 and model groups in both the POS and NEG model. When mice were fed Rb1, the level of phosphatidylcholine increased (Figure 4A). Phosphatidylcholine is a key metabolite in the glycerophospholipid metabolism pathway, and its levels in the liver were affected by related genes. Using RNA-seq and qRT-PCR with liver tissue from the Rb1 and model groups, we found that Cdipt, Dgkg, Pcyt2, Phospho1, Pla2g4b and Plpp5 were enriched in glycerophospholipid metabolism pathway, and they were both dysregulated (Figure 4B). This means that Rb1 affects key enzymes involved in regulating the glycerophospholipid metabolic pathway. Dgkg, Phospho1 and Pla2g4b regulated the level of phosphatidylcholine, which indicates that Rb1 increases level of phosphatidylcholine by modulating the activities of enzymes in the glycerophospholipid metabolic pathway. In addition, comparison of the microbiomes in the Rb1 and model groups revealed that Rb1 alters ether lipid and glycerolipid metabolism in the mouse gut. These two pathways are co-related to glycerophospholipid metabolism. We also found that TC levels were lower in model mice fed Rb1 (Table 1).

Figure 4. Comprehensive analysis of microbiology, lipidomics and transcriptomics. (A) Lipidomics showed content of PC increase in normal mice of Rb1 group comparing with that of Model group. Y axis is shown in foldchange, X axis shows different electron cut-off points of PC. All the metabolites with VIP >1 and P <0.05.PC, phosphatidylcholine. POS, positive ion model. NEG, negative ion model. (B) RNA-seq and qRT-PCR showed expression of RNAs associated with glycerophospholipid metabolism. RNA-seq was performed once, PCR were repeated thrice, n=3, data was shown in mean +/- SD. (C) Serum level of hyperlipidemia associated index showed Rb1 reduced total cholesterol level. TG, triglycerides; TC: total cholesterol; LDL, low density lipoprotein; HDL, high density lipoprotein;** represents P<0.01, n=6. (D) Lipidomics showed content of PC increase in fecal transplanted germ-free mice of Rb1 group comparing with that of Model group. Y axis is shown in foldchange, X axis shows different electron cut-off points of PC. All the metabolites with VIP >1 and P <0.05.PC, phosphatidylcholine. (E) QRT-PCR showed that Pmet increased and Dgkg, Phospho1, Pla4g4b decreased in fecal transplanted germ-free mice of Rb1 group. This meant Rb1 increased synthesis and decelerate decomposition of phosphatidylcholine.

Based on the results summarized above, we predicted that to suppress hyperlipidemia, Rb1 modulates gut microbiota to increase the level of phosphatidylcholine in glycerophospholipid metabolism, decreasing TC. To test the idea, three groups of germ-free mice were administered by gavage feces from control, model and Rb1 mice. We subsequently found that serum TC levels were significantly reduced in mice receiving feces from the Rb1 mice (Figure 4C), while phosphatidylcholine was increased (Figure 4D). In addition, levels of Pemt mRNA were increased and levels of Phospho1 and Pla2g4b were decreased (Figure 4E). These findings support our prediction that, in these mice, Rb1 modulates gut microbiota to regulate the synthesis and decomposition of phosphatidylcholine in glycerophospholipid metabolism, which in turn decreases serum TC, suppressing hyperlipidemia (Figure 5).

Figure 5. RB1 increase phosphatidylcholine in glycerophospholipid metabolism to prevent hyperlipidemia. Through multi-omics analysis, we found Glycerophospholipid metabolism was a common pathway between Rb1 and model group. By using germ-free mice, we found Rb1 increased the abundance of Turicibacter and decreased the abundance of Blautia and Allobaculum to increase content of phosphatidylcholine by raising Pemt and reducing Pla2g4b and Phospho1, and resulting the downregulation of total cholesterol to play a prevention role in hyperlipideamia. Notes: The green box in mocrobe and metabolism box represents the content reduced, the orange box represents the content increased. Pemt promotes synthesis of phosphatidylcholine, Pla2g4b and Phospho1 accelerate decomposition of phosphatidylcholine. Notes: Experiments were repeated twice, data was shown in mean +/- SD, * represent P<0.05, ** represent P<0.01, *** represent P<0.001.

Mice and samples

The Ethics Committee of Liaoning University of Traditional Chinese Medicine approved and supervised the research protocol. Twenty-one mice were divided into 3 equal groups: a control group fed a regular diet, a model group fed a high-fat diet for 10 weeks, and a Rb1 group fed a high-fat diet and Rb1 (200 mg·kg^-1·d^-1). Stool, serum and hepatic tissues were collected. The samples were stored at -80° C immediately after separation.

RNA extraction

Total RNA was extracted from hepatic tissues with TRIzol (Invitrogen, USA) according to the manufacturer’s instructions. The RNA concentration and purity were assessed with an Agilent 2100 Bioanalyzer (Agilent Technologies, USA) and an RNA 6000 Nano LabChip Kit (Agilent Technologies, USA).

Microbiome sequencing process and analysis

Stool samples from 3 mice in the control group, 5 mice in the model group and 5 mice in the Rb1 group were used for microbiome sequencing. After the amplification of the V3 and V4 regions, a library was quantified and pooled, after which Agilent 2100 Bioanalyzer was used to evaluate the library. Microbiome rRNA sequencing was performed with a MiSeq instrument. TrimGalore was used remove base pairs of less than 20 base pairs and the possible adapter sequences. Then short sequences less than 100 bp in length were removed, and FLASH2 was used to obtain the merge sequence and remove low-quality sequence. Mother was used to research and remove the primers in the sequence. Usearch was used to remove sequences with a total base error rate greater than 2 and sequences with lengths less than 100 bp in order to get clean, high-quality, reliable reads for subsequent bioinformatics analysis. Operational taxonomic units (OTUs) were used to measure the similarity between two sequences, and the threshold was set at 97%. OTU taxonomic annotation was obtained after sequence mapping to a database. TaxonTree in the R software package was used to reveal the genus distribution of each species in all samples; only relative abundances greater than 1% are shown. Based on the OTU abundance table, variance decomposition performed through principal component analysis (PCA) using R software was determined to reflect the differences between samples on a two-dimensional graph. The PICRUSt analysis tool was used to predict and analyze species functions based on amplicon sequencing data.

Lipidomics experiment and analysis

For non-target lipidomics, all 21 mice were used for the lipidomics experiment. Serum samples were thawed at 4° C, after which 30-μl aliquots were mixed with CH2Cl2-CH3OH (containing 0.1 g/L 2,6-tert butyl-4-methyphenol, 0° C), and 600 μL of mixed solvent was added. The mixture was then vortexed, incubated at room temperature and, after 200 μL of water were added, centrifuged. The organic phase was removed, dried with nitrogen, mixed with acetonitrile-isopropanol (CH3CN-IPA, 1:1, V/V), redissolved in 500 μL of mixed solvent, and centrifuged at 10,000 rpm for 10 min. The supernatant was then collected, and 10 μL was collected from each sample and mixed thoroughly to obtain the QC samples. Ultrahigh-performance liquid tandem chromatography quadrupole time of flight mass spectrometry (UHPLC-QTOFMS) was used for metabolite detection in the 3 groups. Data were obtained in positive ion mode (POS) and negative ion mode (NEG). The three-dimensional data, including the peak number, sample name, and normalized peak area were imported into the SIMCA14.1 software package (V14.1, Sartorius Stedim Data Analytics AB, Umea, Sweden) for orthogonal projections to latent structures-discriminate analysis (OPLS-DA). A permutation test was performed to further validate the OPLS-DA model. The first principal component of variable importance in the projection (VIP) was obtained from the OPLS-DA model. The criteria VIP >1 and P value< 0.05 (assessed by Student’s t-test) were used to screen differential metabolites between two comparison groups. KEGG (http://www.genome.jp/kegg/) and MetaboAnalyst (http://www.metaboanalyst.ca/) were utilized to search for the pathways of metabolites.

For target lipidomics, 10 μL of the samples were mixed with 190 μL water, after which 480 μL of extract solution (MTBE:MeOH = 5:1) containing an internal standard was added. The samples were sonicated, centrifuged, and the supernatants were transferred for drying. The dried samples were reconstituted in resuspension buffer (DCM:MeOH:H2O= 60:30:4.5) by sonication. The samples were then centrifuged for LC/MS analysis. UHPLC separation was carried out using a SCIEX ExionLC series UHPLC System. An AB Sciex QTrap 6500+ mass spectrometer was applied for assay development. Typical ion source parameters were: IonSpray Voltage: +5500/-4500 V, Curtain Gas: 40 psi, Temperature: 350° C, Ion Source Gas 1:50 psi, Ion Source Gas 2: 50 psi, DP: ±80V.

Transcriptomic sequencing process and analysis

Hepatic tissues from 3 mice from each group were used for whole transcriptome analysis. Libraries were constructed, and deep sequencing was performed using a HiSeq X instrument. The read counts of each transcript were normalized by FPKM. A fold-change >2 and P<0.05 between two samples were the criteria used to identify differentially expressed transcripts. GO (http://www.geneontology.org) and KEGG analyses were performed to construct annotations for differentially expressed transcripts. P values were used to test the reliability of the analysis.

Quantitative real-time PCR validation

qRT-PCR was performed according to the manufacturer’s instructions. β-Actin was used as an endogenous reference. The 2^−∆∆CT method was used to calculate relative RNA expression levels. Student’s t-tests were used to assess differences; P<0.05 was considered significant. The primer sequences are shown in Supplementary Table 4.

Germ-free mice

Eighteen 6-week-old germ free (GF) male mice were fed in a sterile isolator in the experimental animal center at the Shanghai Shrek Experimental Animal Co., Ltd. The temperature was maintained at 20-22° C, and the humidity was maintained at 50-60%. The water, feed and bedding materials were irradiated with 50 kGy radiation, and the drinking bottle and cage were sterilized under high temperature and pressure at 121° C for 60 min. The day/night cycle was 12/12.

Blood lipid detection

The night before the material was taken, the mice fasted and couldn't help water. On the day of sampling, the mice were anesthetized using 1% pentobarbital (0.2 ml/mouse). Blood was collected from abdominal aorta and centrifuged at 3000 rpm for 30 min, after which the serum was collected and stored at -80° C. Lipid detection was then carried out according to the instructions of the biochemical kit (HDL-C, LDL-C, TG, TC; FuJIFILm Corporation) and the operation process of the instrument (automatic biochemical analyzer; FuJIFILm FDC-7000I).

Hematoxylin-eosin staining

After fixing samples in 4% paraformaldehyde, they were dehydrated through an ethanol gradient, embedded in paraffin in the embedding machine, and cut into 4-μm-thick sections. For HE staining, the sections were deparaffinized with xylene, hydrated through an ethanol gradient, stained with HE, dehydrated in alcohol and xylene, and mounted on slides. Images were acquired using a Nikon Eclipse CI microscope.

Oil red O staining

The tissues were dehydrated in 15% and 30% sucrose and then embedded in paraffin to prepare 8- to 10-μm-thick frozen sections. The frozen sections were fixed, soaked in oil red O dye, counterstained with hematoxylin and sealed in glycerin gelatin. Images were acquired with a Nikon E100 microscope.