Core biomarkers analysis benefit for diagnosis on human intrahepatic cholestasis of pregnancy

Abstract

Background

The pregnant women with intrahepatic cholestasis were at high risk of fetal distress, preterm birth and unexpected stillbirth. Intrahepatic cholestasis of pregnancy (ICP) was mainly caused by disorder of bile acid metabolism, whereas the specific mechanism was obscure.

Methods

We performed proteomics analysis of 10 ICP specimens and 10 placenta specimens from patients without ICP through data-independent acquisition (DIA) technique to disclose differentially expressed proteins. We executed metabolomic analysis of 30 ICP specimens and 30 placenta specimens from patients without ICP through UPLC-MS/MS to identify differentially expressed metabolites. Enrichment and correlation analysis was used to obtain the direct molecular insights of ICP development. The ICP rat models were constructed to validate pathological features.

Results

The heatmap of proteomics analysis showed the top 30 up-regulated and 30 down-regulated proteins. The metabolomic analysis revealed 20 richer and 4 less abundant metabolites in ICP samples compared with placenta specimens from patients without ICP, and enrichment pathways by these metabolites included primary bile acid biosynthesis, cholesterol metabolism, bile secretion, nicotinate and nicotinamide metabolism, purine metabolism and metabolic pathways. Combined analysis of multiple omics results demonstrated that bile acids such as Glycohyocholic acid, Glycine deoxycholic acid, beta-Muricholic acid, Noncholic acid, cholic acid, Gamma-Mercholic Acid, alpha-Muricholic acid and Glycochenodeoxycholic Aicd were significantly associated with the expression of GLRX3, MYL1, MYH7, PGGT1B, ACTG1, SP3, LACTB2, C2CD5, APBB2, IPO9, MYH2, PPP3CC, PIN1, BLOC1S1, DNAJC7, RASAL2 and ATCN3 etc. The core protein ACAT2 was involved in lipid metabolic process and animal model showed that ACAT2 was up-regulated in placenta and liver of pregnant rats and fetal rats. The neonates had low birth weight and Safranin O-Fast green FCF staining of animal models showed that poor osteogenic and chondrogenic differentiation of fetal rats.

Conclusion

Multiple metabolites-alpha-Muricholic acid, beta-Muricholic acid, Glycine deoxycholic acid and Glycochenodeoxycholic Acid etc. were perfect biomarkers to predict occurrence of ICP. Bile acids were significantly associated with varieties of protein expression and these proteins were differentially expressed in ICP samples. Our study provided several biomarkers for ICP detection and potential therapeutic targets for ICP development.

Introduction

Intrahepatic cholestasis of pregnancy (ICP) was mainly caused by the disorder of bile acid metabolism. The accumulating toxic bile acids in liver and maternal-fetal circulation led to gestational pruritus and adverse perinatal effects including preterm labour, fetal distress, meconium-stained amniotic fluid and intrauterine death [1]. The risk of fetal stillbirth increased with total bile acid level ≥ 100 µmol/L [2]. The classification of bile acid covered primary bile acid–cholic acid (CA), chenodeoxycholic acid (CDCA), α-muricholic acid (α-MCA), β-muricholic acid (β-MCA), secondary bile acid–deoxycholic acid (DCA), lithocholic acid (LCA), hyodeoxycholic acid (HDCA), ursodeoxycholic acid (UDCA) and conjugated bile acid–taurocholic acid (TCA), glycocholic acid (GCA), taurolithocholic acid (TLCA), glycholithocholic acid (GLCA), glychodeoxycholic acid (GDCA), taurodeoxycholic acid (TDCA), taurochenodeoxycholic acid (TCDCA), glycochenodeoxycholic acid (GCDCA), tauroursodeoxycholic acid (TUDCA), glycoursodeoxycholic acid (GUDCA) [3, 4]. Shao et al. indicated that DCA, TDCA, TCA, GDCA and GLCA were elevated in ICP patients by analyzing 15 types of serum bile acid profiles [5], moreover, GCA was found to be biomarker for ICP diagnosis [6]. However, the mechanism about bile acid disturbance and how protein affected metabolism were still unclear.

Bile acids, the products of cholesterol catabolism, were excreted from liver into intestine and broken down further by the gut microbiota [7, 8]. Many proteins were involved in the process of bile acids conversion directly or indirectly. Cholesterol 7α-hydroxylase was the rate-limiting enzyme in the classical bile acid synthesis pathway, similarly sterol-27 hydroxylase had an important role in the alternative pathway. Sterol 12α-hydroxylase encoded by CYP8B1 was the key enzyme of CA and CDCA generation [9]. Bile acid-mediated signaling pathway activation were implicated in tissue remodeling, energy balance, glucose/lipid homeostasis, intestinal function, biliary physiology, gut microbiome, aging, immune homeostasis and liver regeneration [10]. Farnesoid X receptor (FXR), pregnane X receptor (PXR), vitamin D receptor (VDR), Takeda G protein-coupled receptor 5 (TGR5) and sphingosine-1-phosphate receptor 2 (S1PR2) were responsive receptors of bile acids [11,12,13]. In vivo, bile acid transporter-sodium-dependent taurocholic co-transporting polypeptide (NTCP), organic anion transporting polypeptides (OTAP), fatty acid binding protein (FABP), multidrug resistance-associated protein1/2/3/4 (MDR1/2/3/4), bile salt export pump (BSEP), apical sodium-dependent bile acid transporter (ASBT) and ileal bile acid binding protein (IBABP) were responsible for uptake and efflux of bile acid [14]. Impaired bile acid transporters in liver decreased bile acid flow and consumption which caused the elevated bile acid concentrations in ICP patients, and women with ICP were more prone to have metabolic dysfunction and developed gestational diabetes mellitus (GDM) [15]. In our study, we reported many proteins and metabolites of placentas through proteomic and metabolomic analysis in ICP patients, and these products were closely related to ICP development and hardly been studied.

Women with ICP had many adverse perinatal outcomes, and bile acid as a toxic substance played a predominant role in fetal injury by cellular membrane damage, mitochondrial dysfunction and reactive oxygen species (ROS) generation [16]. Wu et al. reported that bile acid had a damage and oxidative stress on mouse placenta which led to the edema and apoptosis of trophoblasts, and the FXR agonist W450 protected the placenta from impairment [17]. Moreover, obeticholic acid (OCA) reduced the expression of placental proinflammatory genes by activating FXR signaling pathway to alleviate lipopolysaccharide-induced intrauterine growth restriction and fetal death in mice model [18]. Nonetheless, there was still no effective treatment for ICP. Ursodeoxycholic acid (UDCA), as a first-line drug in clinical treatment, obviously relieve maternal pruritus, but had fewer benefits on neonates which was confirmed by randomized placebo-controlled study [19]. Our study substantiated many ICP-related proteins and metabolites and this research opened a new field that might fuel novel opportunities for ICP therapy.

In this study, we analyzed proteins expression and metabolites variation of ICP placenta samples. Our results showed that compared with placenta form patients without ICP, proteins encoded by PSAT1, KRT72, RABGAP1L and ACAT2 etc. were up-regulated in ICP placenta, and MYH7, F7, BLOC1S1, GPLD1, SP9 and TSFM-encoded proteins were down-regulated. Metabolomics analysis showed that glycine deoxycholic acid, GCDCA, α-MCA, β-MCA, LPC (0:0/16:1), LPC (16:1/0:0), Arg-Thr, 3-N-Methyl-L-Histidine, IlE-Leu, Cyclo (His-Pro), Noncholic acid, CA, Pterine, Glycohyocholic acid, B-Nicotinamide Mononucleotide, S-methyl-L-thiocitrulline, Gamma-Mercholic Acid, AG-183, Abietic acid and 2-Furanoic Acid were more abundant in ICP placenta, Carnitine C7:0, Oleamide, (E, Z)-2-Amino-3,14-octadecadien-1-ol and Guanine the opposite. Multi-omics results showed that metabolic pathways, purine metabolism, bile secretion and glyoxylate and dicarboxylate metabolism were the common pathway. The results of clinicopathologic parameters showed that ICP was associated with early delivery gestational age, low fetal weight and gestational diabetes mellitus, moreover, PIN1, ACAT2, PPIL3, INTS10, DCD, TMEM258, KRT72, RABGAP1L, PPP3CC, MYH7, BLOC1S1, CDC42EP4, CRYAB and MYH2 were related to ICP development. Animal model showed that fetal rats of ICP group were much smaller, and Safranin O-Fast green FCF staining showed that poor osteogenic and chondrogenic differentiation of fetal rats. On the whole, our study on placenta of ICP might provide promising therapeutic targets to improve adverse perinatal outcomes of ICP patients.

Materials and methods

Placental specimens

The diagnostic threshold of ICP was the total serum bile acid values was non-fasting 10umol/L in clinical practice [20]. The inclusion criteria of placental specimens: 1) Diagnosed cases of ICP; 2) Availability of placental specimens from ICP pregnancies; 3) Appropriate control placental samples from patients without ICP. The exclusion criteria of placental specimens as follows: 1) Placental samples from pregnancies with other liver or biliary tract disorders; 2) Placental samples from pregnancies with congenital abnormalities; 3) Incomplete clinical data or missing key details about the pregnancy. In our study, we collected 60 placenta specimens with clinicopathological parameter derived from 30 patients without ICP and 30 ICP patients (Supplementary file 1). All patients had signed informed consent forms before samples collection, and this research was approved by the Guangzhou Women and Children’ Medical Center Ethics Committee. The collected samples were fixed in formalin and paraffin-embedded or maintained at -196 °C (liquid nitrogen) for further study. In this study, all procedures were followed in accordance with the declaration of Helsinki. All methods were carried out in accordance with relevant regulations and ARRIVE guidelines.

Quantitative proteomics of Data-independent acquisition (DIA) and metabolomics

In this study, 10 samples from patients without ICP and 10 ICP samples were randomly selected for quantitative proteomics of DIA of the 60 placental samples collected, and all collected samples were analyzed by metabolomics based on broadly targeted metabolome techniques. In proteome analysis, the Q Exactive Plus mass spectrometer (MS) coupled with tandem EASY-nLC 1200 liquid phase was used for data acquisition, and the accurate parameters were as follows: analysis column (50 μm*15 cm, C18, 2 μm, 100Å), mobile phase (A: 0.1% formic acid, B: 0.1% formic acid and 80% ACN), velocity (300 nL/min), MS1 (R = 70 K, AGC = 3e6, Max IT = 30ms, scan range = 350–1250 m/z), MS2 (R = 17.5 K, AGC = 1e6, Max IT = 50ms) and collisional energy (value = 28). In metabolome analysis, the data was collected by Ultra Performance Liquid Chromatography (UPLC) and Tandem mass spectrometry (MS/MS), and the T3 liquid phase acquisition parameters were as follows: chromatographic column (Waters ACQUITY UPLC HSS T3 C18 1.8 μm, 2.1 mm*100 mm), mobile phase (A: ultrapure water and 0.1% formic acid, B: ACN and 0.1% formic acid), gradient of elution (water/ACN; V/V; 10:90 23 min; 95:5 26.1 min), velocity (0.4 ml/min), column temperature (40 °C), sample injection volume (2ul). The Amide liquid phase acquisition parameters were as follows: chromatographic column (Waters ACQUITY UPLC BEH Amide 1.7 μm, 2.1 mm*100 mm), mobile phase (A: ultrapure water, 20mM ammonium formate and 0.4% ammonium hydroxide, B: pure ACN), gradient of elution (water/ACN; V/V; 40:60 9 min; 60:40 21 min; 10:90 26.1 min), velocity (0.4 ml/min), column temperature (40 °C), sample injection volume (2ul). The condition of T3 and Amide acquisition: electrospray ionization (500 °C), ion spray voltage 5500 V (positive), 4500 V (negative), ion source gas (I:55psi; II:60psi), curtain gas (25psi), collision-activated dissociation (high). The differentially expressed proteins were identified according to the fold change (FC) of protein expression (FC > 1.2 or FC < 0.83, P < 0.05). Significantly regulated metabolites between groups were determined by variable importance in projection (VIP) ≥ 1 and absolute Log2FC (fold change) ≥ 1. VIP values were extracted from OPLS-DA results. The mass spectrometry proteomics data had been deposited to the ProteomeXchange Consortium via the iProX partner repository with the dataset identifier PXD031472.

Different maps of differential metabolites

Dynamic distribution and VIP maps of metabolite were drawn according to FC and VIP value respectively. The Volcano Plot was mainly used to demonstrate the relative content difference of metabolites and the significance of the statistical difference in different groups. Z-score standardization map was used to standardize the relative contents of different metabolites by calculating Z-values, and its calculation formula was: z = (x-µ) / σ, x: the relative content of metabolites, µ: the average value, σ: the standard deviation. The violin plot was a combination of boxplot and density plot, mainly used to display data distribution and its probability density.

Enrichment analysis and KEGG annotation

Heatmaps were drawn by R packages-heatmaply (v 1.2.1) and ComplexHeatmap (v 2.7.1.1009). Identified metabolites were annotated using KEGG compound database, subsequently, annotated metabolites and proteins were then mapped to KEGG pathway database. Significantly enriched pathways are identified with a hypergeometric test’s p-value. The map of Differential Abundance Score (DA Score) was a pathway-based method to analyze metabolic changes which could capture the average and overall changes of all metabolites in a pathway. DA score = (a-b)/c, a: number of up-regulated differential metabolites in this pathway, b: number of down-regulated differential metabolites, c: number of all metabolites annotated to this pathway. MSEA enrichment analysis was based on MetaboAnalyst database.

Orthogonal partial least squares-discriminant analysis

Orthogonal Partial Least Squares-Discriminant Analysis (OPLS-DA) was a multivariate statistical analysis, and the specific method was to extract the components of independent variable X and dependent variable Y respectively, and then calculate the correlation between the components.

Correlation analysis of differential metabolites

There were synergistic or mutually exclusive relationships between different metabolites or metabolites and proteins, and correlation analysis could help measure proximities which was conducive to further understanding of the mutual regulatory relationships. The Pearson correlation coefficient was performed by R packages-base package (v 3.5.1) and Hmisc (v 4.4.0). Chord diagram were drawn by R packages-igraph (v 1.2.4.2) and ggraph (v 2.0.2).

Receiver operating characteristic curve analysis of Differential metabolites

Receiver Operating Characteristic (ROC) curve was a quantitative way of making accurate decisions. In this study, the area under ROC curve was used for assessing the discriminative ability of prediction models and finding biomarkers for ICP occurrence.

O2PLS analysis

O2PLS model was used for integrated analysis between proteomics and metabolomics. In this study, we selected all the differential proteins and metabolites to establish O2PLS model. The important variables with high correlation and weight were determined by the load diagram.

Animal models and immunohistochemical staining

The animal models were constructed according to the literature [21]. The antibodies used in immunohistochemistry staining were as follows: PIN1 (10495-1-AP; 1:500; Proteintech, USA), ACAT2 (14755-1-AP; 1:100; Proteintech, USA), PPIL3 (15671-1-AP; 1:100; Proteintech, USA), RABGAP1L (13894-1-AP; 1:100; Proteintech, USA), MYH7 (22280-1-AP; 1:100; Proteintech, USA), BLOC1S1 (19687-1-AP; 1:100; Proteintech, USA), INTS10 (15271-1-AP; 1:50; Proteintech, USA), PPP3CC (55163-1-AP; 1:50; Proteintech, USA), DCD (DF13367; 1:100; Affinity, USA), KRT72 (bs-16831R; 1:100; Bioss, China) and TMEM258 (NBP1-91711; 1:200; Novus, USA).

Safranin O-fast green FCF staining

The kit of Safranin O-Fast green FCF staining was purchased from reagent company (B01115, baiqiandu, China). The sections were observed under a microscope (NIKON ECLIPSE CI, Japan).

Statistics analysis

Results expressed as mean ± standard error (m ± SE) were analyzed by two-tailed Student’s t-test. Fisher’s test was used for clinicopathological parameters and proteins. P < 0.05 using the Prism software (v 9.4.1, GraphPad Software, USA) was considered as significant difference data.

Results

Quantitative proteomics of data-independent acquisition (DIA)

The schematic summary was shown in Fig. 1A. Quantitative proteomic analysis of 10 samples from patients without ICP and 10 ICP samples was performed by Q Exactive plus mass spectrometer coupled with an EASY-nLC 1200 system. The raw MS data were processed using DIA-NN software (v 1.7.16) using library-free method. Briefly, human protein sequences database from SwissProt was used for library prediction by deep learning algorithms. MS resulted in the identification of 77,979 peptides and 5613 proteins. OPLS-DA sensitive to variables with less correlation was used to assess clustering of the different groups. The OPLS-DA results demonstrated that R²X = 0.367, R²Y = 0.993 and Q² = 0.302, and p value of R²Y and Q² was less than 0.05 (Fig. 2A). In total, there were 77 up-regulated proteins and 99 down-regulated proteins between NC and ICP group, and the heatmap displayed top 30 up-regulated and 30 down-regulated proteins with p < 0.05 (Fig. 2B). The top 30 up-regulated differential expression proteins were PSAT1, HBG1, HBM, SPI1, PIP4K2B, HBG2, HBE1, FOXK1, YOD1, KRT72, PKLR, SIGLEC6, SH3BP4, BTN3A3, BIN2, HMBS, CASKIN2, SMIM1, SLC2A3, PTGIS, ERVFRD-1, NOP10, PADI3, MRC2, UBE2O, SLC46A1, TMC4, RABGAP1L, ACAT2 and SERPINB5 orderly (Table 1). MYH7, F7, BLOC1S1, GH2, GPLD1, SP9, TSFM, MBD2, SLC13A3, EGFL7, EVA1A, HSPA13, SMARCB1, RYBP, RIOX2, PPP3CC, SULT2B1, TTC9C, COQ6, SDC1, GPR173, PPTC7, KRT84, COL5A3, LIFR, RASAL2, RIPK2, TAOK1, CAMSAP2 and MLST8 were down-regulated in sequence (Table 2).

Fig. 1

Schematic summary. (A) Schematic outlining sample acquisition, mass spectrometry platforms and subsequent analysis process

Fig. 2

Quantitative Proteomics of Data-independent acquisition (DIA). (A) The OPLS-DA permutations Test diagram of ICP differential expression proteins. (B) The heatmap of top 30 up-regulated proteins and top 30 down-regulated proteins in proteomics

Table 1 Top 30 up-regulated differential expression proteins in NC vs. ICP placentas

Table 2 Top 30 down-regulated differential expression proteins in NC vs. ICP placentas

Quantitative metabolites

The overlap diagram of total ions current (TIC) showed that the retention time and peak intensity of TIC curve were consistent in different quality control (QC) samples which indicated a stable signal (Supplementary file 2: Fig. 1A). Coefficient of Variation (CV) distribution map of samples displayed that the proportion of QC samples with CV value less than 0.3 was higher than 85% (Supplementary file 2: Fig. 1B). The correlation diagram of QC sample verified that the Pearson correlation coefficients between different samples were greater than 0.85 (Supplementary file 2: Fig. 1C). Model validation permutation test diagram of OPLS-DA showed that R2X = 0.25, R2Y = 0.938 and Q2 = 0.578, and p value (R2Y and Q2) < 0.005 in multivariate statistics analysis of metabolomics (Fig. 3A). The volcano plot emphasized 20 up-regulated metabolites and 4 down-regulated metabolites by double screening of Log2FC and VIP values among 762 metabolites (Fig. 3B). The heatmap implied 24 differential metabolites- GCDCA, α-MCA, Gamma-Mercholic Acid, CA, Noncholic acid, β-MCA, glycine deoxycholic acid, Glycohyocholic acid, 3-N-Methyl-L-Histidine, IlE-Leu, Arg-Thr, Cyclo (His-Pro), S-methyl-L-thiocitrulline, Oleamide, (E, Z)-2-Amino-3,14-octadecadien-1-ol, LPC (16:1/0:0), LPC (0:0/16:1), Guanine, B-Nicotinamide Mononucleotide, AG-183, Carnitine C7:0, Pterine, 2-Furanoic Acid and Abietic acid. Among them, GCDCA, α-MCA, Gamma-Mercholic Acid, CA, Noncholic acid, β-MCA, glycine deoxycholic acid and Glycohyocholic acid belonged to bile acids; 3-N-Methyl-L-Histidine, IlE-Leu, Arg-Thr, Cyclo (His-Pro) and S-methyl-L-thiocitrulline were classified to amino acid and its metabolites; Oleamide and (E, Z)-2-Amino-3,14-octadecadien-1-ol belonged to alcohol and amines; LPC (16:1/0:0) and LPC (0:0/16:1) belonged to glutamate pyruvate (GP); Guanine and B-Nicotinamide Mononucleotide were classified to nucleotide and its metabolites; AG-183 was benzene and substituted derivatives; Carnitine C7:0 was fatty acid (FA); Pterine was heterocyclic compounds; 2-Furanoic Acid was classified to organic acid and its derivatives (Fig. 3C). The Z-value diagram of differential metabolites was used to see the distribution of each metabolite in NC and ICP group intuitively, and the content of α-MCA and β-MCA had good intra-group stability (Fig. 3D). The Violin plot of raw values showed that 20 up-regulated metabolites were glycine deoxycholic acid, GCDCA, α-MCA, β-MCA, LPC (0:0/16:1), LPC (16:1/0:0), Arg-Thr, 3-N-Methyl-L-Histidine, IlE-Leu, Cyclo (His-Pro), Noncholic acid, CA, Pterine, Glycohyocholic acid, B-Nicotinamide Mononucleotide, S-methyl-L-thiocitrulline, Gamma-Mercholic Acid, AG-183, Abietic acid and 2-Furanoic Acid in ICP placenta, and Carnitine C7:0, Oleamide, (E, Z)-2-Amino-3,14-octadecadien-1-ol and Guanine were down-regulated (Fig. 3E).

Fig. 3

Quantitative metabolites. (A) The OPLS-DA permutations Test diagram of differential metabolites. (B) The volcano plot. (C) The heatmap of 24 differential metabolites. (D) The Z-value diagram of differential metabolites. (E) The Violin plot of differential metabolites

Enrichment Analysis and correlation of metabolites

The differential metabolites in descending order of FC were Glycine deoxycholic acid, GCDCA, β-MCA, α-MCA, LPC (0:0/16:1), LPC (16:1/0:0), Arg-Thr, 3-N-Methyl-L-Histidine, IlE-Leu, Cyclo (His-Pro), Noncholic acid, CA, Pterine, Glycohyocholic acid, B-Nicotinamide Mononucleotide, S-methyl-L-thiocitrulline, Gamma-Mercholic Acid, AG-183, Abietic acid, 2-Furanoic Acid, Carnitine C7:0, Oleamide, (E, Z)-2-Amino-3,14-octadecadien-1-ol and Guanine (Fig. 4A). The differential metabolites in descending order of VIP score were β-MCA, α-MCA, Glycine deoxycholic acid, GCDCA, LPC (0:0/16:1), LPC (16:1/0:0), IlE-Leu, Arg-Thr, S-methyl-L-thiocitrulline, 3-N-Methyl-L-Histidine, AG-183, Pterine, Abietic acid, Glycohyocholic acid, 2-Furanoic Acid, Noncholic acid, CA, Gamma-Mercholic Acid, Oleamide, (E, Z)-2-Amino-3,14-octadecadien-1-ol, B-Nicotinamide Mononucleotide, Carnitine C7:0, Cyclo (His-Pro) and Guanine (Fig. 4B). The specific values of VIP and FC were shown in Table 3. There were synergistic or mutually exclusive relationships between different metabolites, and correlation analysis could help measure metabolic proximities of significantly different metabolites, which was conductive to further understanding of the mutual regulatory relationships between metabolites in the process of biological state changes [22]. The heatmap and chord map results of differential metabolites demonstrated that there was a strong correlation between B-Nicotinamide Mononucleotide and Arg-Thr, LPC (0:0/16:1) and LPC (16:1/0:0), β-MCA and α-MCA, Glycine deoxycholic acid and GCDCA, Oleamide and (E, Z)-2-Amino-3,14-octadecadien-1-ol, Noncholic acid, CA and Gamma-Mercholic Acid respectively (Fig. 4C and D). The DA Score map included six pathways- Nicotinate and nicotinamide metabolism, Primary bile acid biosynthesis, Cholesterol metabolism, Bile secretion and Purine metabolism. The dots were distributed on the middle or right side of the central axis indicting that the overall expression of pathway was more inclined to up-regulated (Fig. 4E). MSEA enrichment analysis results implied that Retinol metabolism and Riboflavin metabolism were involved in ICP development (Fig. 4F). Small Molecule Pathway Database (SMPD) was designed to support pathway elucidation in metabolomics [23]. Zellweger Syndrome, Familial Hypercholanemia (FHCA), Congenital Bile Acid Synthesis Defect II/III, Cerebrotendinous Xanthomatosis (CTX), Bile Acid Biosynthesis and 27-Hydroxylase Deficiency were enriched by SMPD (Fig. 4G). To further explore the biological significance, we performed ROC curve analysis, and the results showed that the area under ROC curve of β-MCA / α-MCA / Glycine deoxycholic acid / GCDCA was close to 1, followingly, the area under ROC curve (AUC) of Arg-Thr / Glycohyocholic acid / LPC (0:0/16:1) / LPC (16:1/0:0) / IlE-Leu / 3-N-Methyl-L-Histidine was greater than 0.8 among the up-regulated metabolites (Fig. 5A). The AUC of (E, Z)-2-Amino-3,14-octadecadien-1-ol and Oleamide was 0.710 among four down-regulated metabolites (Fig. 5B).

Fig. 4

Enrichment Analysis and Correlation of Metabolites. (A) Fold change distribution diagram of metabolites. (B) The variable importance in projection distribution diagram of differential metabolites. (C) The correlation heatmap of differential metabolites. (D) The chord diagram of differential metabolites. (E) The diagram of differential abundance score. (F) The MSEA enrichment diagram of differential metabolites. (G) The HMDB enrichment diagram of differential metabolites

Table 3 Differential metabolites in NC vs. ICP placentas

Fig. 5

Receiver Operating Characteristic (ROC) curve analysis. (A) ROC curve of 20 up-regulated metabolites. (B) ROC curve of 4 down-regulated metabolites

Multi-omics Integration Analysis of Placental Data

There were 48 pathways enriched by 60 placental specimens and 223 pathways in the protein enrichment pathway. The 20 common regulatory were Steroid hormone biosynthesis, Vitamin B6 metabolism, Drug metabolism-other enzymes, Lysine degradation, Purine metabolism, Butanoate metabolism, Biotin metabolism, Tryptophan metabolism and Pyrimidine metabolism etc. (Fig. 6A). Among the proteomics and metabolomics analysis of 20 placental specimens, there were 4 metabolites and 26 proteins involved in metabolic pathways. Moreover, there were 2 metabolites-GCDCA and Acetaminophen, and 2 proteins encoded by EPHX1 and ADCY4 involved in bile secretion pathway with p value < 0.1 (Fig. 6B). To further elucidate the interaction between proteomics and metabolomics, we selected all differential proteins and metabolites to establish O2PLS model in order to screen out core variables in multi-omic analysis through load diagram. Metabolite loading diagram showed that the top 10 metabolites affecting proteomics were β-MCA, α-MCA, Glycine deoxycholic acid, GCDCA, LPC (0:0/16:1), LPC (16:1/0:0), IlE-Leu, S-methyl-L-thiocitrulline, Pterine and Acetaminophen respectively (Fig. 6C). Protein loading diagram showed that the top 10 proteins affecting metabolomics were encoded by ACAT2, CDC42EP4, DCD, FECH, CRYAB, INTS10, KRT72, PPIL3, RABGAP1L and TMEM258 respectively (Fig. 6D). Further, we made a network diagram of the five metabolites involved in the four common pathways and their related regulatory proteins (Fig. 6E). ICP was a disease related to bile acid metabolism disorder, and our analysis showed that bile secretion was the only pathway with p value < 0.1. To explore the mechanism of ICP occurrence, we illustrated the relationship between core metabolites (bile acids - GCDCA, α-MCA, Gamma-Mercholic Acid, CA, Noncholic acid, β-MCA, glycine deoxycholic acid and Glycohyocholic acid) and their related proteins. The results showed that glycohyocholic acid was positively related to GLRX3, PGGT1B, SP3, LACTB2, C2CD5 and IPO9, negatively related to MYL1, MYH7, ACTG1, APBB2 and MYH2 expression. BLOC1S1, PPP3CC, DNAJC7, RASAL2, PIN1 and ATXN3 were associated with multiple bile acids (Fig. 6F). The above results emphasized that proteins in bile secretion pathway, protein loading diagram or bile acids-related proteins might play an essential role in metabolic disturbance of ICP development.

Fig. 6

Multi-omics Integration Analysis of Placental Data. (A) Venn diagram summarizing the overlapping pathways. (B) The bar graph of KEGG enrichment analysis. The red and green bars represent metabolome and proteome respectively. (C) The loadings diagram of metabolites. (D) The loadings diagram of proteins. (E) The network diagram of metabolite-protein. The solid line represents a positive correlation and the dashed line represents a negative correlation. (F) The expression correlation analysis between metabolite and protein

Core Proteins Analysis in ICP Development

We further analyzed the expression levels of all bile acids-related proteins in proteomics, and the results verified that the expression of PIN1, TMEM258, KRT72, PPIL3, BLOC1S1, CDC42EP4, INTS10, RABGAP1L, PPP3CC, ACAT2, DCD, MYH7 and MYH2 were significantly different between NC and ICP group (Fig. 7A). Correlation results between ICP and clinicopathologic parameters showed that ICP development was closely to preterm birth (p < 0.0001), low fetal weight (p = 0.0098), gestational diabetes mellitus (p = 0.0303), the expression of PIN1, ACAT2, PPIL3, INTS10, DCD and TMEM258 (p = 0.023), the expression of KRT72 (p < 0.0001), the expression of PPP3CC and MYH7(p < 0.0230), and expression of RABGAP1L and BLOC1S1 (p = 0.0011), moreover, CDC42EP4, CRYAB and MYH2 weren’t associated with ICP occurrence in our study (Table 4).

Fig. 7

The expression of core proteins. (A) Expression of core proteins. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001

Table 4 Correlation between intrahepatic cholestasis of pregnancy and clinicopathologic parameters

Validation of animal models

We constructed animal models of ICP, and the results showed that fetal rats were smaller in terms of morphology in ICP group (Fig. 8A). HE staining of fetal long bone showed that there was reduced trabeculae in ICP group (Fig. 8B). Additionally, Safranin O-Fast green FCF staining demonstrated that fetal rats in ICP group had poor osteogenic and chondrogenic differentiation (Fig. 8C). ACAT2 was involved in multiple metabolic processes, such as fatty acid metabolism and pyruvate metabolism (Supplementary file 2: Fig. 2). We further explore the expression of ACAT2 in liver, and the results verified that ACAT2 was up-regulated in maternal and fetal liver in ICP group (Fig. 8D). PIN1, ACAT2, PPIL3, INTS10, DCD, TMEM258, KRT72, PPP3CC, MYH7, RABGAP1L and BLOC1S1 were the core proteins in ICP development validated by women placental specimen. Immunohistochemical staining of rat placenta also demonstrated that BLOC1S1, PPP3CC, CRYAB, MYH7 and MYH2 were down-regulated in ICP group (Fig. 8E). PIN1, TMEM258, KRT72, PPIL3, INTS10, RABGAP1L, ACAT2 and DCD were up-regulated in ICP group (Fig. 8F).

Fig. 8

Validation of Animal Models. (A) The general view of fetal rats. (B) The HE staining of long bone. (C) The diagram of Safranin O-Fast green FCF staining. (D) The immunohistochemical staining of liver. (E) The immunohistochemical staining of down-regulated proteins in ICP placenta. (F) The immunohistochemical staining of up-regulated proteins in ICP placenta

Discussion

In our study, bile secretion pathway played a vital role in ICP development, and it was confirmed by metabolomics that the metabolism disorder of GCDCA, α-MCA, Gamma-Mercholic Acid, CA, Noncholic acid, β-MCA, glycine deoxycholic acid and Glycohyocholic acid classified into bile acids was the main cause of ICP occurrence. In multi-omics joint analysis, core proteins-PIN1, ACAT2, PPIL3, INTS10, DCD, TMEM258, KRT72, PPP3CC, MYH7, RABGAP1L and BLOC1S1 were related to bile acids metabolism. Moreover, we also elaborated that ICP development was closely to preterm birth (p < 0.0001), low fetal weight (p = 0.0098), gestational diabetes mellitus (p = 0.0303). Furthermore, animal models validated that fetal rat smaller in size in ICP group had poor osteogenic and chondrogenic differentiation. In addition, the expression trend of core proteins in rat placenta were consistent with proteomics results. Collectively, our results provided a new insight for ICP liquid biopsy and therapy.

Biological processes were complex and integrated. Single omics data couldn’t analyze the macro development process of biological systems and it was difficult to explain complex biological network regulation. Combined analysis of proteomics and metabolomics was used in our research. On the one hand, integrating multiple omics data could make up for data problems caused by data missing and noise during single omics data analysis; On the other hand, the combined analysis provided us with a panoramic window to study phenotype and regulation mechanism of biological processes systematically and comprehensively [24, 25]. ICP as a metabolism disorder disease was proved to be associated with gestational diabetes mellitus. Liu et al. also raised that women with ICP were more likely to have gestational diabetes mellitus by a retrospective cohort study of 95,728 singleton births and the adjusted odds ratio was 1.406 [26]. Mechanically speaking, FXR missing in ICP might impress normal glucose homeostasis [27].

Bile acid metabolites - GCDCA, α-MCA, Gamma-Mercholic Acid, CA, Noncholic acid, β-MCA, glycine deoxycholic acid and Glycohyocholic acid with differential expression were included in our study. Bile acid synthesis homeostasis was regulated by FXRs expressed in the liver and gut, and CA-FXR agonistic bile acid, α-MCA and β-MCA - FXR antagonistic bile acids were reported in a range of disease [28]. Frank J Gonzalez and his colleagues summarized that activation and inhibition of FXR was associated with obesity, insulin resistance and non-alcoholic fatty liver disease [29]. Maria J Perez demonstrated that GCDCA as the oxidative insult could lead to mitochondrial function impairment and cell apoptosis of fetal and maternal hepatocytes by ascending Bax-alpha/Bcl-2 level and caspase-3 activity [30].

Core proteins-PIN1, ACAT2, PPIL3, INTS10, DCD, TMEM258, KRT72, PPP3CC, MYH7, RABGAP1L and BLOC1S1 were associated with ICP development. Moreover, these core proteins also regulated the progression of other pregnancy-related disease. For instance, there was a good communication between PIN1 and craniosynostosis [31], ACAT2 and birth weight and neonatal lean mass [32], DCD and placental function [33]. Some core proteins were only reported in different tumor, such as PPIL3 in breast cancer [34], RABGAP1L in estrogen receptor-positive breast cancer [35], INTS10 in hepatocellular carcinoma [36], PPP3CC in epithelial ovarian cancer [37]. TMEM258 could regulate endoplasmic reticulum homeostasis and was reported to be associated with inflammatory bowel disease risk and impaired secretion of glycoproteins [38]. KRT72 was related to a wide range of epithelial disorders [39] and BLOC1S1 was a critical mediator of autophagy [40, 41]. Nonetheless, the study of these core proteins in ICP development was still blank at present, and our study provided a new direction for ICP progress. But we used a small sample size, which limits the generalizability of our results and we will enlarge the sample size to do more in-depth research on the mechanism of ICP in the future.

Conclusion

Summarily, our study was the first to illustrate the potential molecular mechanism of ICP development by virtue of proteomics and metabolomics analysis. The results of clinicopathologic parameters showed that ICP group had early delivery gestational age, low fetal weight and susceptibility to gestational diabetes mellitus compare with NC group. The enrichment analysis verified up-regulated bile secretion in ICP placenta. Bile acid metabolites - GCDCA, α-MCA, Gamma-Mercholic Acid, CA, Noncholic acid, β-MCA, glycine deoxycholic acid, Glycohyocholic acid and bile acid metabolism related protein- PIN1, ACAT2, PPIL3, INTS10, DCD, TMEM258, KRT72, PPP3CC, MYH7, RABGAP1L and BLOC1S1 might be potential drug therapeutic targets for ICP to improve poor perinatal outcomes.