Whole-exome sequencing identifies novel mutations in ABC transporter genes associated with intrahepatic cholestasis of pregnancy disease: a case-control study

Background Intrahepatic cholestasis of pregnancy (ICP) can cause premature delivery and stillbirth. Previous studies have reported that mutations in ABC transporter genes strongly influence the transport of bile salts. However, to date, their effects are still largely elusive. Methods A whole-exome sequencing (WES) approach was used to detect novel variants. Rare novel exonic variants (minor allele frequencies: MAF < 1%) were analyzed. Three web-available tools, namely, SIFT, Mutation Taster and FATHMM, were used to predict protein damage. Protein structure modeling and comparisons between reference and modified protein structures were performed by SWISS-MODEL and Chimera 1.14rc, respectively. Results We detected a total of 2953 mutations in 44 ABC family transporter genes. When the MAF of loci was controlled in all databases at less than 0.01, 320 mutations were reserved for further analysis. Among these mutations, 42 were novel. We classified these loci into four groups (the damaging, probably damaging, possibly damaging, and neutral groups) according to the prediction results, of which 7 novel possible pathogenic mutations were identified that were located in known functional genes, including ABCB4 (Trp708Ter, Gly527Glu and Lys386Glu), ABCB11 (Gln1194Ter, Gln605Pro and Leu589Met) and ABCC2 (Ser1342Tyr), in the damaging group. New mutations in the first two genes were reported in our recent article. In addition, compared to the wild-type protein structure, the ABCC2 Ser1342Tyr-modified protein structure showed a slight change in the chemical bond lengths of ATP ligand-binding amino acid side chains. In placental tissue, the expression level of the ABCC2 gene in patients with ICP was significantly higher (P < 0.05) than that in healthy pregnant women. In particular, the patients with two mutations in ABC family genes had higher average values of total bile acids (TBA), aspartate transaminase (AST), direct bilirubin (DBIL), total cholesterol (CHOL), triglycerides (TG) and high-density lipoprotein (HDL) than the patients who had one mutation, no mutation in ABC genes and local controls. Conclusions Our present study provide new insight into the genetic architecture of ICP and will benefit the final identification of the underlying mutations. Supplementary Information The online version contains supplementary material available at 10.1186/s12884-021-03595-x.


Background
Intrahepatic cholestasis of pregnancy (ICP) is a reversible pregnancy-specific liver disease characterized by pruritus and abnormal liver function, such as elevated liver enzymes and increased serum total bile acids (TBA) (≥ 10 μmol/L), that appears in the second and third trimesters of pregnancy and resolves completely after delivery in the early postpartum period [1]. The incidence of ICP disease is reported to be between 0.1 and 15.6% depending on geographical differences [1,2]. ICP has been associated with adverse fetal outcomes, including spontaneous preterm birth, respiratory distress, low Apgar scores, fetal distress and fetal death [3][4][5][6][7]. For example, Li Li et al. reported that birth size was reduced, whereas the rate of small-for-gestational-age infants was increased across increasing categories of serum total TBA [8]. It has been noted that approximately 2-4% of ICP pregnancies are affected by fetal mortality [9,10]. The level of TBA increases the risk of fetal morbidity and stillbirth [11][12][13]. Therefore, untangling the mechanisms of ICP and its association with fetal complications is very important.
The abnormal synthesis, metabolism transport, secretion and excretion of bile acids may lead to ICP disease [14]. Therefore, the etiology contributing to the development of ICP disease is complex and depends on multiple factors, including hormonal, genetic and environmental backgrounds [15]. Familial clustering analysis in pedigree studies showed a high incidence in mothers and sisters of patients with ICP, which might indicate a genetic predisposition for the condition [16][17][18][19]. Among these, gene mutations in the hepatocellular transporters of bile salts play a pivotal role in the pathogenesis of ICP [20].
Bile salt transport is the key physiological function of ATP-binding cassette (ABC) membrane proteins, covering seven distinct members: ABCA, ABCB, ABCC, ABCD, ABCE, ABCF and ABCG. Of these genes, ABCB4, ABCB11 and ABCC2 are functionally known genes that affect the development of ICP [21][22][23]. Except for the ABCB4, ABCB11, ABCC2, ATP8B1, TJP2 and ANO8 genes [23,24], the role of other ABC transporter genes seems to be less studied. By taking advantage of the high-throughput genotyping technologies in a larger-scale population, the WES approach that combined genotype data for all patients has proven to be far more efficient in anchoring exonic mutations for the target gene at once. In particular, the method can greatly accelerate screening for new potential pathogenicity sites of all mutations. Therefore, examining exonic variants across ICP disease groups likely augments the excavation of novel loci. However, to the best of our knowledge, there have been no reports of the use of WES to identify the genetic variants in the ABC gene series of bile acid transporters for ICP disease.
ICP is a complex disease that is influenced by multiple genes. Previously, we reported the identification of the novel gene ANO8 as a genetic risk factor for intrahepatic cholestasis of pregnancy in 151 ICP samples [24]. In the present study, we hypothesized that genetic variation in ABC transporters confers susceptibility to ICP. Therefore, we performed WES to extensively investigate the presence of mutations, especially for the discovery of new functional variants, of the ABC gene series involved in bile acid transport in 151 patients with ICP disease and related them to the clinical data and pregnancy outcomes.

Patients and clinical data
A total of 151 pregnant women, who were also used for our previous study [24], were included who did not have other liver diseases and were diagnosed as ICP on the basis of skin pruritus in combination with abnormal liver biochemistry indexes, e.g., TBA, ALT and AST. Among them, the cut-off level of TBA was 10 μmol/L. In addition, 26 clinical features covering six main indicators of maternal and neonatal data, including basic patient features (age, gestational age and body mass index: BMI), ion concentrations (K, Na, Cl, Ca, Mg, and P), routine blood tests (white blood cell count: WBC, red blood cell count: RBC, platelet count: PLT, and red blood cell distribution width.SD: RDW.SD), liver function indexes (TBA, ALT, AST, total bilirubin: TBIL, DBIL, and indirect bilirubin: IDBIL), lipid indexes (CHOL, TG, HDL, low-density lipoprotein: LDL, and uric acid) and outcomes of pregnant women and newborns (birth weight and bleeding amount) were recorded. The concrete methods for these biochemical variables can be found in our previous study [24]. Briefly, the ion concentration, liver function and lipid indexes were determined by an AU5800 automatic biochemical analyzer (Beckman Coulter). Routine blood tests were determined by a Sysmex-xn-2000 automatic blood cell analyzer. We collected all the samples from June 2018 through July 2020. Blood samples were collected from 17 to 41 + 3 weeks of gestation. The delivery gestational age of pregnant women ranged from 28 to 41 + 3 weeks (median: 38 weeks). In addition, we also recruited 1029 pregnant women without ICP disease as negative controls in the same periods, which were also used in our previous study. The controls were in the same gestational age range as the patients with ICP.

Whole-exome sequencing
A total of 151 genomic DNA samples were extracted from peripheral blood by an Axy Prep Blood Genomic DNA Mini prep Kit (Item No. 05119KC3). To fulfill the requirements for genotyping, the concentration and integrity of all the extracted DNA samples were determined by a Nanodrop-1000 spectrophotometer (Thermo Fisher, USA) and gel electrophoresis, respectively. A total of 151 ICP patient samples were subjected to exome capture with a BGI Exon Kit according to the manufacturer's protocols. DNA libraries were constructed using the combinatorial probe anchor ligation (cPAL) method. Each resulting qualified captured library was then loaded on BGISEQ-500 sequencing platforms. In addition, whole-exome sequencing was also performed in 1029 controls.

Variant annotations, filtering and prioritizing
The bioinformatics analysis began with the sequencing data. We carried out the analysis mainly according to Huusko et al. [25]. Briefly, first, joint contamination and low-quality variants (read depth < 15 and genotype quality score < 20) were removed. The reads were mapped to the reference human genome (UCSC GRCh37/hg19) using BWA (Burrows-Wheeler Aligner) software [26]. After that, variant calls were conducted by GATK (v3.7) [27]. Finally, the ANNOVAR tool was applied to perform a series of annotations [28], such as the frequency Fig. 1 Overview flow chart of the whole-exome sequencing study for intrahepatic cholestasis of pregnancy in the databases, conservative prediction and pathogenicity prediction with available website tools for variants.
The quality control workflow was carried out as shown in Fig. 1. First, we removed variants with an MAF greater than 0.01 in the 1000 Genomes Project (http:// www.internationalgenome.org/), Exome Aggregation Consortium (ExAC) (http://exac.broadinstitute.org/), and dbSNP (https://www.ncbi.nlm.nih.gov/snp/) databases. Second, missense, damaging or loss-of-function variants of ABC genes were included in the subsequent analyses. Third, we concentrated on novel variants that were filtered by NCBI and Ensembl. In addition, we prioritized paying attention to novel variants that would likely have functional effects. For instance, a variant was highlighted when it was predicted to be simultaneously deleterious by SIFT, Mutation Taster and FATHMM software.

Statistical analysis
We applied the sapply function to calculate descriptive statistics for clinical data. The Pearson correlation coefficient estimation between clinical data was evaluated by using the cor function [29]. The cor.test function was used to examine the significance of the correlation coefficients. Fisher's test was used to test the significance of the frequency between 151 patients with ICP and 1029 controls in 3 databases [30]. Moreover, the pie function was used to draw the percentage of all ABC mutations. The comparison of average value differences in six clinical datasets (TBA, AST, DBIL, CHOL, TG and HDL) among the four groups was performed by one-way ANOVA [31]. All the analyses were preformed using R software. In addition, we used three web-available tools, namely, SIFT, Mutation Taster and FATHMM, to predict protein damage. Predictions were defined as damaging, probably damaging, possibly damaging, or neutral when all three prediction tools, two out of three, one out of three and none reached the damaging level, respectively.

Verification of novel candidate sites by sanger sequencing
Since the WES approach is associated with false positive results, the presence of 13 interesting novel mutations from WES analyses was confirmed using Sanger sequencing. The thirteen loci were selected for sequencing mainly based on functional prediction results for mutations (the damaging, probably damaging, possibly damaging, and neutral groups) and classification of members of the ABC transporter gene superfamily (ABCA, ABCB, ABCC, ABCD, ABCE, ABCF and ABCG). The primers were designed by Primer Premier 5 software. The details of the PCR primers and their optimum annealing temperatures are shown in Table 1.

Evolutionary conservation analysis
We also selected the same twelve loci as Sanger sequencing for evolutionary conservation analysis for the same reasons. Evolutionary conservation analysis [32,33] of the amino acids encoded by the new pathogenicity sites in ABC genes was performed among vertebrates, including Gibbon, Macaque, Olive Baboon, Gelada, Marmoset, Mouse, Rat, Cow, Sheep, Pig, Chicken, Zebra Finch, and Zebrafish, using genomic alignments of the Ensembl Genome Browser.

Protein structure modeling
Protein structure modeling involves in two steps. First, the reference and modified (ABCC2 Ser1342Tyr) protein sequences were submitted to SWISS-MODEL (http:// swissmodel.expasy.org/) software for structure modeling. Then, these two protein models were compared simultaneously using the Chimera 1.14rc package.

Results
Whole-exome sequencing results of variants of ABC transporter genes in 151 ICP samples In total, ABC transporters covering 44 genes underwent targeted WES to identify variants in a cohort of 151 patients with ICP disease. Overall, we identified 2953 genetic variants, including 1254 mutations in 12 ABCA genes that consisted of ABCA1-ABCA10 and ABCA12-ABCA13, 479 variants in ABCB1 and ABCB4-ABCB11, 812 genetic variants in eleven genes covering ABCC1-ABCC6 and ABCC8-ABCC12, and 408 variants in the ABCD-ABCG gene series. These types of variants covered 2057 introns, 297 synonymous, 76 splice, 13 stop gained/start lost, 36 3′ primer UTR, 21 5′ primer UTR, 10 downstream gene, 32 upstream gene, 3 structural interaction and 408 missense. The percentage of these types of variants is shown in Fig. 2a. After quality control (MAF < 0.01 in all databases), 42 out of 320 variants were detected that were first reported, including 17 out of 146 variants in ABCA genes, 10 out of 59 in ABCB, 9 out of 71 in ABCC and 6 out of 44 in ABCD-ABCG genes (Fig. 2b). For these 42 novel mutations, we classified these mutations as four echelons (damaging, probably damaging, possibly damaging and neutral) according to the prediction results ( Table 2). The damage group had fifteen mutations, which contained six mutations in three known functional genes associated with ICP disease, such as ABCB4 Trp708Ter, Gly527Glu and Lys386Glu and ABCB11 Gln1194Ter, Gln605Pro and Leu589Met. These six novel mutations have been shown in our recently published article [24]. In addition, 9      10E-14). Twenty-one mutations were assigned to the second tier (probably damaging), including 9 in ABCA genes (ABCA2 Asp1108Glu and Ala583Val, ABCA5 Val997Ala, ABCA7 Pro449His, ABCA8 Val8Ala, ABCA10 Leu260Ser, ABCA12 Ile1022Thr, ABCA13 Leu4624Val and Thr4912Ala), Glu386Lys and Ile339Val in all ABCB9 genes, 6 in the ABCC gene series (ABCC1 Ser497Gly, ABCC3 Leu137Phe, ABCC5 Gln851Pro, ABCC6 His1043Gln, ABCC9 Glu1034Val and ABCC12 Pro37Ser), ABCD4 Ser395Phe and Thr385Ala, ABCF2 Gly47Ser and ABCG1 Ile273Val in the ABCD-ABCG gene series. These mutations were also absent in controls. The MAF in these mutations between cases and controls showed a significant difference (P = 2.20E-16).
We identified five possible potential pathogenic mutations associated with ICP disease, including ABCA3 Ser593Gly, ABCA10 Leu823Val, ABCA12 Asn2492Ser, ABCA13 Ser3111Asn and ABCG1 Thr378Ile. Only one mutation was placed in the neutral group. All 42 novel mutations were absent in the databases, including the 1000 Genome Project, ExAC, dbSNP, ChinaMAP and 1029 local controls.

Confirmation of the novel variants by sanger sequencing
We used Sanger sequencing to confirm 13 possible candidate novel pathogenicity loci in the ABC family genes from four lines. The results (Fig. 3) were all consistent with WES.

Evolutionary conservation analysis
We performed evolutionary conservation analysis of twelve loci, which were also sequenced by Sanger sequencing. The results suggested that these mutations were highly conserved among vertebrate species, including rats, sheep, cows, pigs, dogs, horses, chickens and so on (Fig. 4).
Comparison of the protein structural model of the ABCC2 Ser1342Tyr mutation ABCC2 consists of 32 exons involved in bile formation that mediate hepatobiliary excretion of numerous organic anions and conjugated organic anions such as methotrexate and transport sulfated bile salts such as taurolithocholate sulfate [34,35]. This gene has two main nucleotide binding domains located at positions 671-678 and 1334-1341 (Fig. 5a). We only identified the ABCC2 Ser1342Tyr mutation when the quality control of the minor allele frequency was below 0.01. The location of this variant (Ser1342Tyr) is close to the ATP binding functional domain (1334-1341). In addition, the mutation was found to be harmful to the function of the protein according to the results of three prediction software programs. To further investigate the possible effects of the missense variant on protein structure, the reference and modified protein structures of the ABCC2 gene were compared simultaneously using UCSF Chimera 1.14rc. The results showed that, compared with the reference molecular structure, the 3D model of mutation had a slight change in the chemical bond lengths of ATP-ligand binding amino acid side chains at positions Ser1342, Ser678 and Gln706 (Fig. 5b). The change in the amino side chains could affect the binding efficiency of the ATP molecule.
In this study, we did not collect placental tissue from these 151 patients with ICP. Therefore, to further analyze the genetic basis of ABCC2, we analyzed the mRNA expression level of the ABCC2 gene in placental tissue between 2 healthy participants and four patients with ICP using GEO datasets derived from NCBI (GEO accession: GSE46157) from Du Q et al.'s report [36]. A significant (P < 0.05) difference in gene expression was observed between the two groups (Fig. 5c). The expression of ABCC2 was upregulated in the ICP group. In addition, we also detected that the expression of three other genes, namely, ABCC6, ABCE1 and ABCG5, changed in placental tissue (Fig. 6).

Correlations among clinical data
We examined the correlation coefficients (Supplementary file 1) among these 26 clinical features and found that TBA was significantly negatively correlated with gestation days (r = − 0.34), birth weight (r = − 0.30), and RDW.SD (r = − 0.16) and significantly positively correlated with TBIL (r = 0.52), DBIL (r = 0.56), IDBIL (r = 0.18), ALT (r = 0.17), and AST (r = 0.18). There was also a close relationship among liver function indexes. For example, ALT was highly positively correlated with AST (r = 0.93). ALT and AST were positively correlated with TBIL and DBIL. TBIL positively correlated with DBIL (r = 0.89). CHOL, TG, LDL and uric acid were also positively correlated with ALT and AST. In addition, the concentrations of Ca and Mg ions were positively correlated with AST. The correlations between the abovementioned clinical data were significant (P < 0.05).

Biochemical and clinical features of ICP cases with variants
Descriptive statistics of 26 clinical features for patients with ICP with 42 new mutations are shown in Table 3.  For all the clinical data, the levels of ALT, AST, TBA, DBIL, CHOL and TG in 151 patients with ICP disease were higher than the reference levels. Notably, individuals with the novel ABC mutations had a fourfold higher TBA and twofold higher ALT, AST, and TG average values than the reference values, which confirms that ICP disease presents with abnormal liver function and elevated bile acids associated with abnormal lipid metabolism. Compared to the ICP group without the 42 novel mutations, the group with the 42 novel mutations had a higher age, BMI, and levels of K, Na, Cl, RBC, TBIL, IDBIL, TG, and HDL. Additionally, we found that six patients with ICP containing two mutations exhibited higher TBA, AST, DBIL, CHOL, TG and HDL than 31 patients with one mutation, 114 patients with ICP with no mutations and 414 local healthy controls (Fig. 7). In particular, for TBA, as a clinical characteristic of ICP, the trends for the average values measured in patients with ICP with mutations of ABC transporter genes and local healthy controls were ranked as follows: ICP with two mutations > ICP with one mutation > ICP with no mutations > healthy controls.

Discussion
This study is the first to use whole-exome sequencing technology to uncover potential novel pathogenic mutations in ABC family genes involved in bile acid transport. Altogether, we identified 42 new loci covering 44 ABC genes in 151 patients who were diagnosed with ICP disease.
The present study has 4 major strengths. First, following the advent of WES technology, it has proven to be efficient in unearthing sequence variation across the MAF spectrum in obstetric and gynecological diseases [25,37]. Using this method, we successfully identified novel candidate pathogenicity loci with known functional genes ABCB4, ABCB11 [24] and ABCC2. In addition to these three genes, we also identified some novel variations in other genes in the ABC series. Our results revealed the genetic mutations of the first WESidentified ABC family genes associated with ICP disease. Second, to date, no studies have unraveled the genetic mutations in ABC genes of hepatic disease among pregnant women from a relatively large nationally representative sample (n = 151) in China. Using this same population, we previously identified a novel gene, ANO8, associated with ICP. Moreover, the clinical data of these patients are relatively complete, providing data support for association analysis between mutations and clinical data and subsequent functional verification. Finally, combining WES and clinical data favors deciphering the molecular mechanism of ICP disease.
Certainly, our study also has limitations. First, the WES approach needs a larger sample size to target lowfrequency and rare variants. Otherwise, it may miss some valuable variants. Moreover, it also results in inaccurate MAFs for rare variants. However, the strict exclusion conditions guaranteed the selection of a defined cohort, such as employing the MAF in the databases and identifying 151 ICP cases and using 1029 local controls, combined with prediction tools. Second, the samples in (See figure on previous page.) Fig. 5 Effects of the ABCC2 Ser1342Tyr variant on the protein sequence and structure. a Distribution of the ABCC2 variants. ABCC2 is a 1545-amino acid protein containing two major nucleotide-binding functional domains (ATP binding cassette domains). The locations of two ABCC2 variants from whole-exome sequencing are shown in the protein sequence. b In silico comparison of the higher assembly of the reference and modified (Ser1342Tyr) ABCC2 protein models. The reference and Ser1342Tyr molecules are presented as gold and blue rounded ribbon structures, respectively. The enlarged image shows that the functional domains have small changes in the chemical bond lengths. DBD: DNA-binding domain. c Comparison of the expression level of the ABCC2 gene between two healthy controls and four patients with ICP. Different letters above the expression column denote significant differences (P < 0.05) Fig. 6 Comparison of the expression levels of the ABCC6, ABCE1 and ABCG5 genes in placental tissue between two healthy controls and four patients with ICP this study were all derived from Jiangxi Province, which lacks geographic diversity and may limit the applicability of the generality of these results. However, this study is still valuable due to the high incidence (~5%) of ICP in Jiangxi. Next, the frequency of damaging mutations in 151 ICP samples was 10.60%, accounting for only a small portion of the genomic mutations of ICP disease. Therefore, we need to identify more genes/mutations associated with ICP. Finally, our results provided some possible interesting causative mutations; however, the causality between these potential pathogenic candidate loci and ICP disease needs to be verified by validation functional experiments.
Previous studies have detected genetic mutations in ICP primarily by Sanger sequencing or offspring studies based on a limited number of individuals. The rare  [25]. To our knowledge, this study is the largest analysis to date of the role of mutations in the ABC series of genes in 151 ICP-susceptible patients. Intriguingly, seven novel pathogenic loci in the ABCB4, ABCB11 and ABCC2 genes were simultaneously assigned to the damage group, suggesting the accuracy of our results. Several studies have shown that heterozygous missense mutations in ABCB4, especially lowfrequency and rare variations, are commonly responsible for the occurrence of ICP disease, which is consistent with our results [38][39][40]. After filtering the frequency of the data, we identified a total of eight mutations with a MAF < 0.01, seven of which were heterozygous missense mutations and one of which was a nonsense mutation in ABCB4. The role of ABCB11 in ICP has also been clearly identified, although its contribution seems less than that of ABCB4. A comprehensive analysis of multiple previous studies concluded that up to 5% of ICP cases harbor monoallelic mutations in ABCB11 [41]. Our study confirmed the role of ABCB11 and further expanded the role of the ABCB11 gene.
The function of ABCC2 is related to ICP disease. Therefore, genetic variants of ABCC2 may cause bile acid metabolism disorder. Our results demonstrated a small number of variations, in particular one novel pathogenic heterozygous (ABCC2 Ser1342Tyr) variant. The study of Corpechot C et al. showed that genetic variants of ABCC2 contribute to inherited cholestatic disorders [42]. Moreover, More et al. also found that the expression of ABCC2 was more closely related to the livers from an alcohol cirrhosis cohort, indicating that ABCC2 expression changed in liver-related disease [43], which is consistent with the fact that ICP is a liver disease. For the other three genes that were expressed differently, we found ABCC6 to have a novel mutation, His1043Gln, which was located in the probably damaging group. Therefore, the ABCC6 gene and the His1043Gln locus are of interest in the development of ICP disease. In our study, we failed to detect novel mutations in the ABCE1 and ABCG5 genes; therefore, we conjectured that this is probably because of the low frequency of novel mutations. It might also be that the known mutations in these two genes are associated with ICP, or other patterns of mutations lead to changes in the mRNA expression levels of the two genes. Moreover, Fig. 7 Effects of ABC gene mutations on biochemical indexes. The difference in the average TBA, AST, DBIL, CHOL, TG and HDL in the following four groups: a having both mutations, b one mutation, c no mutation in ABC series genes with ICP, and d healthy controls without ICP disease. Means of different groups in rows with different superscript letters are significantly (P < 0.001) different from each other compared to healthy samples, the levels of ABCE1 and ABCG5 mRNA expression were upregulated and downregulated in ICP, respectively. A reasonable explanation for this condition is that the unique functions of ABCE1 and ABCG5 lead to different mechanisms of ICP disease.
Because there have been relatively few studies on the genetic basis of ICP disease, many functional genes have not been mined. However, based on previous studies, several researchers found some evidence that certain genes may predispose patients to liver or ICP disease, such as ABCB9 associated with hepatocellular carcinoma and ABCG2 associated with ICP [44,45]. Therefore, we hypothesized that genetic mutations in these unknown functional genes, excluding ABCB4, ABCB11 and ABCC2, confer susceptibility to ICP, especially in the damaging group. These findings extended the knowledge on the role of mutations and strengthened the understanding of the genetic basis of ICP disease, including many additional ABC family genes.
Bacq et al. suggested that elevations in serum TBA, ALT and AST activity were identified among patients with ICP who have mutations in ABC transporter genes [46]. In agreement with this, our results confirmed that the mutation group containing 42 novel variants was associated with higher TBA, AST, DBIL, CHOL, TG and HDL, especially among the patients with ICP with 15 damaging novel mutations (Table 3, Fig. 7). Furthermore, Piatek K et al. suggested that the analysis of genotype coexistence pointed to the possibility of mutated variants of polymorphisms of genes having a summation effect on the development of ICP disease [47]. Consistent with this result, in our study, we found that six patients with ICP had two mutations at once. Notably, ICP harboring both mutations increased the average TBA level by 17.33, 20.22 and 50.32 μmol/L based on the difference in average TBA levels in ICP with one mutation, no mutation in ABC transporter genes and no mutation in healthy controls, respectively. As expected, AST, DBIL, CHOL, TG and HDL were also increased. A reasonable explanation for this situation is that double mutation has an additive effect on TBA levels, and this accumulation of TBA levels results in ICP exacerbation.
Moreover, significant differences between the wild type and mutant type were seen for the biochemical indexes, which suggested that the accumulation of TBA could lead to lipid abnormalities. Consistent with this result, researchers found that bile acids (BAs), which are well known for their amphipathic nature, are essential to lipid absorption and energy balance in humans [48][49][50]. Thus, taken together, our study provides additional support for the effect of mutations in ABC family genes on ICP disease and lipid metabolism.
Not surprisingly, to date, there are only a handful of reports of genetic analysis studies for ICP disease. The present findings not only enrich the molecular basis of the known functional genes (such as ABCB4, ABCB11 and ABCC2) but also expand the new candidate genes (ABCA, ABCB, ABCC, and ABCD-ABCG) associated with ICP. Our results support the fact that mutations in ABC transporter genes lead to ICP disease. Therefore, further work on the genetic mutations involved in ICP pathogenesis has the potential to inspire novel therapies for patients with ICP.
Recently, there has been much interest in whether variations in bile salt concentrations and mutations in ABCB4/ABCB11/ABCC2/other new candidate genes could be biomarkers for various forms of drug-induced liver injury [51]. Furthermore, considering the negative effect of ABC mutations on maternal and fetal outcomes, it is important to genotype these mutations to genetically diagnose them in a timely manner and provide immediate attention and treatment for all ICPsusceptible people.

Conclusion
To the best of our knowledge, this is the first study to conduct WES to reveal the genetic variants in ABC family genes associated with ICP disease. We detected 42 novel potential pathogenic mutations in 44 ABC family genes. Among them, seven loci were identified in ABCB4, ABCB11 and ABCC2, and the remaining 35 loci were in other genes. In particular, the 15 novel mutations that were classified as the damaging group were the main focus of further study. Their functional validation and experimental verification need to be further investigated. Moreover, we detected genetic variants that were significantly associated with six biochemical indexes, including TBA, ALT, AST, DBIL, CHOL and TG (P < 0.05). Our findings provide new valuable insights into the genetic basis of ICP disease and suggest potential candidate variants for clinical diagnosis.
Additional file 1. Correlation coefficients among twenty-six clinical features. The dots indicate the significant (P < 0.05) correlation coefficients between each pair of features. The size and colors separately represent the degree and direction of correlation coefficients.