Alterations of the Human Gut Microbiota in Intrahepatic Cholestasis of Pregnancy

Background and Aims Women with severe intrahepatic cholestasis of pregnancy (ICP) are at higher risks of fetal complications and without effective treatments. Changes in gut microbiota in pregnancy were found to be related to the altered intestinal bile acid composition, so we aimed to explore the alterations of microbiota in the gut of ICP patients. Methods A total of 90 women were recruited, including 45 ICP patients and 45 healthy controls. The gut microbiota communities of ICP group were compared to control group through 16S ribosomal RNA gene sequencing. The results were then confirmed by real-time polymerase chain reaction (PCR) and generalized linear model (GLM). Furthermore, we analyzed the relationships between microbiota and the severity of ICP. Results A total of seven genera and nine taxa with differential abundances between the ICP patients and the controls were identified. All of the seven genera were verified through real-time PCR, and three key genera Parabacteroides, Flavonifractor, and Megamonas were confirmed by using the GLM model. Further analysis found that the genera Escherichia_Shigella, Olsenella, and Turicibacter were enriched in the severe ICP group, the microbial gene function related to biosynthesis of unsaturated fatty acids and propanoate metabolism were also increased in them. Conclusions Overall, our study was the first in Asia to demonstrate an association between gut microbiota and ICP. Our findings would contribute to a better understanding of the occurrence of ICP.


Illumina MiSeq PE250 sequencing
Amplicons were extracted from 2% agarose gels and purified using the AxyPrep DNA

Process of sequencing data
PANDAseq(V2.9) [2] is employed to assemble pairs of reads into complete sequences Tags, trimmed of barcodes and primers, were further checked on their rest lengths and average base quality. 16S tags were restricted between 250 and 500 nt such that the average Phred score of bases was no worse than 20 (Q20) and no more than 3 ambiguous N. The copy number of tags was enumetated and redundancy of repeated tags was removed. Only the tags with frequency more than 1, which tend to be more reliable, were clustered into OTUs, each of which had a representative tag.
The α-diversity and β-diversity indices were calculated based on the rarefied OTU counts using the Qiime program. We randomly selected 23,410 reads for each sample. The α-diversity represents an analysis of diversity in a single sample reflected by parameters including good coverage, Chao 1, whole tree, Shannon index and Simpson index using Qiime. Wilcoxon rank sum test was used to compare each α-diversity index. β-diversity is used as a measure of the microbiota structure between the two groups. Both the weighted and unweighted Unifrac distance matrices were plotted in the principal coordinate analysis (PCoA), and analyses of similarities (ANOSIMs) were performed using the R "vegan" package.
For taxa with a prevalence≥10%, differential abundance analysis was performed using the Wilcoxon rank-sum test at the phylum, class, order, family, and genus levels.
For multiple comparisons of bacterial counts, the false discovery rate was calculated using the Benjamini and Hochberg method. Microorganism features used to distinguish the microbiotas specific to ICP were identified using the linear discriminant analysis (LDA) effect size (LEfSe) [5] method with an alpha cutoff of 0.05 and an effect size cutoff of 2.0.

Phylogenetic Investigation of Communities by Reconstruction of Unobserved
States (PICRUSt) [6] was used to predict the abundances of functional categories in the Kyoto Encyclopedia of Genes and Genomes (KEGG) orthologs (KO). The graph of KEGG pathways in level 2 (41pathways) and level 3 (328pathways) 2 was performed with STAMP [7] , and p values were calculated with White's non-parametric t-test.
For the PICRUSt analysis, the OTU clustering method we used is de novo OTU clustering. The analysis software used is Usearch (version 7.0), and the database used is RDP database (http://rdp.cme.msu.edu). Firstly, the singletons in the spliced long reads were filtered out, because the singletons may be caused by sequencing errors.
This part of the sequence was removed without cluster analysis. Usearch was used to cluster under the similarity of 0.97, and the clustered sequences were filtered by chimerism to obtain OTUs for species classification.

Genus-specific quantification by real-time PCR
The comparison and analysis of the 16S region sequences with those of differential genera reference strain were carried out by the multiple alignment program DNAman software. The specific primer target sites for the quantitative analysis were performed using Primer 5. Universal 16S rRNA gene was used as the internal reference and the abundances of the genera were expressed as relative levels to 16S rRNA. The PCR reaction and condition were the same with 16S rRNA gene quantification(95 ℃ for 3 min, followed by 30 cycles at 98℃ for 20 s, 58 ℃ for 15s, and 72 ℃ for 20 s and a final extension at 72 ℃ for 5 min). The genus-specific primer sequences used were listed in TableS1. Data (three biological repeats and three technical repeats at each time point) were analyzed by the comparative threshold cycle (CT) method and the standard formula. Each reaction was performed in triplicate.