Analysis of the relationship between bile duct and duodenal microbiota reveals that dysbacteriosis is the main cause of primary common bile duct stones

Background Bacteria play an important role in the formation of primary CBD stones. While the composition and function of the microbiota of bile duct in patients with primary CBD stones could not be explained in the whole. Intestinal dysbacteriosis is related to the formation of many diseases, but there are relatively few reports on the effects of intestinal microbial community on the formation of primary CBD stones, but none of these studies have a normal population of biliary tract and gut microbiota as control. Result We use the 16S rRNA gene high-throughput sequencing technology to analyze the microbial diversity and community composition of biliary and duodenal microbiota in 15 patients with primary CBD stones (EG) and 4 patients without biliary tract disease (CK), and we divided above samples into four subgroups: B (EG bile), D (EG duodenal juice), BCK (CK bile)and DCK (CK duodenal juice). Alpha diversity analysis showed that only the abundance and evenness in B were lower than that in DCK (P<0.05); Beta diversity analysis showed that there were some differences in the composition between biliary microbiota and the duodenal microbiota, but the abundance of the main groups showed similarities. Proteobacteria and Firmicutes were the dominant bacteria at phylum level accounted for at least 75% of the total reads in each subgroup. It is noteworthy that Clostridium sensu_stricto (VIP>1), Lachnospiraceae _UCG-008 (VIP>1), Butyrivibrio (VIP>1) and Roseburia (VIP> 1) which could produce short chain fatty acids (SCFA), were signicantly decreased in B (p<0.05) than that in BCK. Conclusion According to the results, dysbacteriosis of biliary and duodenal microbiota plays an important role in the formation of primary CBD stones. The in abundance of certain major acid-producing bacteria affects the health biliary leads to the formation of stones.

The previous studies could not give an overall explanation of the composition and function of the biliary microbiota, especially the normal biliary microbiota, due to technical reasons. What's the relationship between biliary and gut microbiota and how do they affect the formation of primary CBD stones have rarely been reported. With the development of science and technology, high-throughput sequencing based on 16S rRNA gene PCR ampli cation has been widely used in microbial research because of its low cost, high throughput, high sensitivity, and rapid sequencing [9].
Based on the 16S rRNA gene high throughput sequencing technique, the bile and duodenal juice of the patients with primary CBD stones and CK in our hospital are sequenced using the Illumina Miseq sequencing platform, the composition and diversity of the microbiota were analyzed, and the differences between the bile duct bacteria and intestinal bacteria were compared, to explore its signi cance for the formation of primary CBD stones.

Results
Through the metagenomic DNA extraction, 16S rRNA gene ampli cation and amplicon sequencing of the 38 biliary bile and duodenal juice samples from 19 persons including 15 patients and 4 controls, a total of 1,811,427good reads were obtained and 2,801 OTUs were annotated. Sequencing data statistics and diversity analysis each sample yielded 47669.13 sequences (range = 33338 ~ 76885, SD = 10186.88) on average.
The good reads coverage of each group was above 99%, and there was no signi cant difference in the amount of sequence and OTUs at phylum, genus and species among the four subgroups, namely EG group bile (B subgroup) and duodenal juice (D subgroup) and CK group bile (BCK subgroup) and duodenal juice (DCK subgroup).
The Species accumulation curves [10]for the total number of OTUs corresponding to each sample were draw. It showed when the sample size increases, the slope of the curve gradually decreases and tends to atten. At which point the total number of OTUs in the community will no longer increase signi cantly with the addition of new samples. It indicates that the current sample size was enough to re ect the richness of the community.
We applied the diversity index (Chao1, ACE, Simpson and Shannon) to re ect its alpha diversity. The Shannon curve in each sample reached plateau phase, which con rmed the sequence numbers is large enough to cover most of the microbial information in the sample, and abundance and evenness of the microbial species in the current study meet the requirements for sequencing and analysis. Although the Simpson, Chao 1 and ACE index demonstrated there was no striking difference in the microbial community abundance between groups, the Shannon index showed signi cant difference between B and DCK group (p < 0.05). Shannon index is more sensitive to microbiota richness and rare OTUs. Thus, it indicated that the microbiota richness in B group is lower than that in DCK, and the difference is mainly in the content of the rare species (Fig. 1).
Beta diversity analyses were used by employing Principal component analysis (PCA) and Multidimensional scaling (MDS) of the similarity of microbial community structure in different groups. The PCA [11] result was shown in Fig. 2, partial samples of each group gathered on the left side of the longitudinal axis (dotted line).The BCK group and DCK group were distributed on the horizontal axis (dotted line) and gathered together. While the samples of the B subgroup and D subgroup showed more discrete and distributed in different regions. Mapping data were obtained from UniFrac PCoA master coordinate analysis via QIIME [12,13].Based on the Unweighted UniFrac distance, almost all bile samples were on the left side of the vertical axis and the duodenal uid samples were on the right side ( Fig. 3.a), and Non-metric Multidimensional scaling (NMDS) analysis based on Unweighted UniFrac distance obtained similar mapping data (Fig. 3.c). On the contrary, based on the Weighted UniFrac distance, the distribution of all samples is more concentrated and some of them are clustered and overlapped ( Fig. 3.b), and NMDS analysis based on Weighted UniFrac distance also showed similar results (Fig. 3.d).
We presented the proportion of shared and unique OTUs in each subgroup by using the Venn diagram visually (Fig. 4).There were 1205 shared OTUs among the four subgroups, occupying 43% of all OTUs (total 2801 OTUs). More than 70% OTUs were shared in B subgroup and D subgroup, but only about 60% shared in BCK subgroup and DCK subgroup. The most dominant shared phylum and genus were Proteobacteria(95%) and Pseudomonas(83%), respectively. In the Venn diagram, we could also observe the unique OTUs in each subgroup. The unique OTUs of genus in the B subgroup were mostly Pseudomonas, which corresponding species were unidenti ed. Only 6 unique OTUs of genus were observed in BCK subgroup: Dechlorobacter, Halomonas, Pseudomonas and other unidenti ed genera.
Analyses the microbial community structure in different groups According to the OTU division and classi cation status, the speci c microbial composition of samples at each classi cation level can be obtained. A total of 23 phyla were identi ed in this work. At each sample or group, Proteobacteria and Firmicutes accounted for at least 75% of all phyla (Fig. 5). The number of bacteria in B subgroup was less than that of the D subgroup and DCK subgroup at the phylum, genus and species level (Table 1). A total of 291 genera were identi ed in the two groups, and the dominant genera were Pseudomonas, Escherichia-Shigella and Streptococcusin sequence. 14 genera were signi cant different with p < 0.05 between two groups. The abundance of Escherichia-Shigella, Geobacillus, Stomatobaculum, Campylobacter and Solibacillus in D subgroup was signi cantly increased, the remaining 9 genera such as Ochrobactrum, Moraxella, etc. increased signi cantly in DCK subgroup (Fig. 6.c).
3 Comparison of the microbial community structure in B subgroup and D subgroup Similarly, Proteobacteria and Firmicuteswere still the most abundant phyla. The bacterial phyla in bile were contained in the duodenal uid, however the duodenal uid additionally contained Nitrospirae and Gracilibacteria. Eight phyla were signi cant different with p < 0.05 between B subgroup and D subgroup. In addition to the abundance of Proteobacteria and Cyanobacteria increased, the others decreased signi cantly in bile ( Fig. 6.d).
There were 295 genera were identi ed in the two groups, of which Pseudomonas, Escherichia-Shigella and Clostridiumsensu_stricto_1were the dominant genera in sequence. There are 15 unique genera in the bile, but there are up to 40 unique genera in the duodenal uid. Even more striking is that 69 genera were signi cant different with p < 0.05 between two groups ( Fig. 6.e). Among 69 genera, the abundance of Pseudomonas was signi cantly increased in the bile, others such as Streptococcus, Veillonella, Lactobacillus, Prevotella_7\Prevotella_6, etc. increased signi cantly in duodenal uid.
4 Comparison of the microbial community structure in BCK subgroup and DCK subgroup A total of 18 phyla were identi ed in the two groups, of which Proteobacteria, Firmicutes and Actinobacteria were the dominant phyla in sequence. Only phylum Bacteroidetes was found obvious difference between bile and duodenal uid in CK, of which the abundance was signi cantly increased in the duodenal uid (p < 0.05) ( Fig. 6.f).
At the genus level, we identi ed 234 genera from the two groups. Similarly, Pseudomonas and Escherichia-Shigellain the phylum Proteobacteria are in turn the two most dominant. 17 genera were found obvious signi cant difference between two groups ( Fig. 6.g), of which only the abundance of Aquabacterium and Enhydrobacter was signi cantly increased in the duodenal uid.
Key species screening and association analysis 1 Partial Least Squares Discriminant Analysis (PLS-DA) The PLS-DA discriminant model was constructed based on the species abundance matrix and sample grouping data. Just like Fig. 7 showed, both bile samples and duodenal uid samples received relatively better grouping models. At the same time, at the Genus level, we calculated the VIP (Variable importance in projection) coe cient for each genus (VIP > 1, the greater the value, the greater the contribution of the species to differences between groups). In bile samples, there were 75 genera that contributed to the bile grouping. Correspondingly, 87contributing genera were discriminated for duodenal uid grouping in the duodenal uid.

Spearman Association Network Analysis of Dominant genera
We calculated the Spearman rank correlation coe cient between the dominant Genera, whose abundance was located in the top 50, and constructed the associated network for the related advantage where rho > 0.6 and P value < 0.01 [15].The results showed ( Fig. 8) that there were abundant correlated networks between EG and CK. Most of the bacterial genera existed in the same environment with collaborative manners. However the associations among bacterial genera in EG were more complicated. In CK, the relationship between the genera Ralstonia and Ochrobactrum, Streptococcus and Granulicatella, Acinetobacter and Alloprevotella were antagonistic to each other( Fig. 8.a), but they showed cooperative relationships in EG( Fig. 8.b).

Microbial function prediction
In order to further understand the differences in the functions of biliary and intestinal microbiota, we performed functional prediction of the 16S rRNA gene sequence by PICRUSt in the KEGG database [16].At KEGG level 3, we obtained a total of 304 pathways. 14 pathways were discovered to have signi cant difference in B subgroup and BCK subgroup (p < 0.05), and they were down-regulated in B subgroup than those in BCK subgroup. 8pathwayswere signi cant different with p < 0.05 in duodenal uid samples between D subgroup and DCK subgroup, including 6 up-regulated pathways and 2 down-regulated pathways.
The 6 major metabolic pathways, Metabolism, Genetic Information Processing, Environmental Information Processing, Cellular Processes, Organismal Systems and Human Diseases, all showed differences in some pathways at KEGG level 3 in bile samples between EG and CK. Microbiota in B subgroup had lower abundance of functions involved in Metabolism pathways such as Cysteine and methionine metabolism, Glycerophospholipid metabolism and Tetracycline biosynthesis than those in BCK subgroup. In addition, the abundance of cellular processes (Apoptosisand Endocytosis), environmental information processing (Two-component system and CAM ligands) and the immune system (RIG-I-like receptor signaling pathway) decreased signi cantly in B subgroup ( Fig. 9.a).
In duodenal uid samples ( Fig. 9.b), microbiota in D subgroup had higher abundance of functions involved in Metabolism pathways (Purine metabolism, Metabolism of xenobiotics by cytochrome P450 and Naphthalene degradation) and human diseases (Bacterial invasion of epithelial cells, Measles).But the abundance of Digestive system (Pancreatic secretion) and nervous system (Glutamatergic synapse) decreased signi cantly in D subgroup.

Discussion
Due to differences in lifestyle, dietary structure, and living environment, primary CBD stones maintain a high incidence in China, which has brought about property damage and health hazards for our country and people. As an important condition for the occurrence of the disease, bacteria have been studied by scholars. Most previous experiments were based on Cholecystolithiasis or in vitro model to infer the relationship between the two. In addition, the intestinal microbiota plays a role in human diseases with its rich functions. There are a number of diseases such as Crohn's disease, colon cancer, and metabolic diseases that have been shown to be related to the intestinal ora [5][6][7].The duodenal microbiota, as the nearest neighbor to the biliary tract, was considered to be the main source of possible biliary bacterial infections [17]. In this study, bile and duodenal uid were directly obtained through ERCP techniques and diversity analysis was performed. Through the rigorous screening of the subjects, we obtained a total of 19 samples of patients. It is worth noting that 4 of them were pancreatic-related diseases patients without hepatobiliary system disease. For the rst time, the project used normal normal bile duct bile was used as CK to investigate the relationship between CBD and microbiota. The study was ethical and informed consents were obtained from the patients. To ensure that samples were not contaminated with each other was a di cult point in the sampling process. Thus, the duodenal uid was rst taken and then bile was extracted under the premise of strict aseptic processing to ensure the delity of this study.
In the eld of microbial research, the widespread application of high-throughput sequencing technology has broken the limitations of traditional culture techniques. In our study, the 16S rRNA gene highthroughput sequencing technology was used for biliary and intestinal (duodenal) microbiota in patients with CBD stones and in patients without hepatobiliary system diseases. Through Illumina platform sequencing, a total of 1,811,427 effective sequencing sequences were obtained, and a total of 2801 bacterial OTUs were identi ed, The Goods coverage (%) of each group was above 99%, that is, the amount of sequencing of each sample is nearly saturated. In addition, the Species accumulation curves indicate that the current sample size was enough to re ect the richness of the community. The above two fully illustrate that this study has a very high reliability and reference value. 748,340 high-quality gene sequences were obtained in B subgroup, and 2062 bacterial OTUs were identi ed, which is similar to previous studies by Shen et al [18]. Interestingly, at the same time, we obtained a large number of effective high-quality sequences of 172,025 and 1597 OTUs in BCK subgroup (without hepatobiliary system disease). In the past, based on bacterial culture, it was found that there was no bacteria in normal bile which without biliary tract infections, because of its potent antimicrobial properties [19]. This study has broken through the limitations of this conclusion. Some scholars have also demonstrated the existence of bacteria in bile without common bile duct stones. For example, Pereira P et al. found the presence of bile bacterial communities in PSC patients [20]. Liu et al. found bacterial structures in acute liver failure and normal mice [21]. However, none of the patients enrolled in this study had the above-mentioned liver and biliary tract diseases, which provided valuable resources for the study of normal biliary microbiota.
The sequence quantity and the number of OTUs identi ed indicated that the bile duct bile and duodenal juice are rich in bacteria, regardless of whether the patients have hepatobiliary disease. Our result showed that the Alpha diversity of bile microbiota in the experiment group was lower than that in CK. But there was also study showing that the Alpha diversity of biliary microbiotawas higher than that of normal intestinal ora and had a higher richness [8]. The factors for this difference may be related to the sites of the bile samples (gallbladder or common bile duct) and the sites of the intestinal microbiota samples (fecal or duodenal uid).At the same time, we observed that in the same group, the Alpha diversity of bile and duodenal microbiota was similar, at least for their respective dominant bacterium.
We found that the most dominant phylum in bile and duodenal uid was Proteobacteria in obtained 23 phyla. Other phyla with higher abundances included Firmicutes, Actinobacteria, Bacteroidetes, etc., this is similar to previous studies [4,22].The difference is that Synergistetes (0.63%) with a high abundance in B subgroup were not found in BCK subgroup. Synergistetes, the main pathogens of oral diseases such as periodontitis [23], may also participate in the formation of CBD stones. The analysis of the composition of the biliary and duodenal microbiota showed that phyla contained in the bile could be found in the duodenal uid almost. Similar result was also found in previous work of Ye et al. [4], their work showed the biliary microbiota had a comparatively higher similarity with the duodenal microbiota. However, the abundance of Proteobacteria and Cyanobacteria in the bile of the experiment group was signi cantly higher than that in the duodenal uid. It shows that although the composition of biliary and duodenal bacterial community is similar, but the composition ratio is different.
At the genus level, the most dominant genus was Pseudomonasof the phylum Proteobacteria, followed by Escherichia-Shigella. This is different from previous research results, Ye [4,24].According to Metastatsand PLS-DA, Pseudomonas (VIP < 1) and Escherichia-Shigella (VIP < 1) were discovered to have signi cant difference in bile between EG and CK, and were not critical to bile grouping. In our view, the two most abundant genera may only exist as normal ora in bile, not as a risk factor for CBD stones. At the same time, Escherichia-Shigella (VIP > 1) was found to be more abundant in the duodenal uid of EG than that of CK. It shows that Escherichia-Shigellaplays an important role in the duodenum of patients with CBD stones. The bile or duodenal uid in EG contained more genera than CK. The results showed that although the oras of the samples were relatively conservative at phylum level, there was dysbacteriosis and the imbalance of the composition ratio in the bile or duodenal uid on genus level in EG. Meanwhile, the associations among bacterial genera in EG were more complicated by Spearman Association Network Analysis. It further proves that dysbacteriosis occurs in both the biliary tract and the duodenum in the relationship among bacterial genera.
Studies have shown that changes in the bacterial community of the duodenum were associated with various diseases such as small intestinal bacterial overgrowth (SIBO), irritable bowel syndrome (IBS) and celiac disease (CD) [25]. It was found that Clostridium sensu_stricto (VIP > 1), Lachnospiraceae_UCG-008 (VIP > 1), Butyrivibrio (VIP > 1), Roseburia (VIP > 1), which can produce short chain fatty acids (SCFA), were signi cantly reduced in bile in EG compared with CK. SCFA has the functions of sterilization, bacteriostasis, and provide energy for host antibodies to improve immunity [26,27]. SCFA not only affects the intestinal epithelium, but has a certain effect on airway epithelium, enteric nervous system (ENS), central nervous system (CNS), and peripheral neurons [28]. The occurrence of multiple diseases is associated with a decrease in the proportion of "good bacteria" that can produce SCFA in the intestine, such as in ammatory bowel disease [29], type 2 diabetes [30], Asthma [31] and so on. So we can conclude that the low expression of "good bacteria" in bile has a negative impact on biliary health in patients with common bile duct stones, which ultimately leads to the development of stones. In addition, Klebsiella (VIP > 1) and Helicobacter (VIP > 1) have been reported in previous studies to be associated with CBD stone formation. For example, Klebsiellapneumoniae and Helicobacter pylori have been cultured by bacteria or identi ed by PCR technique in bile from patients with CBD stones [18,32].Similar results were obtained in this study. But the two genera were more abundant in BCK subgroup. This phenomenon suggests that the disorder of the biliary bacteria community may be more important than the number of related pathogenic bacteria in the formation of CBD stones.
The comparison of biliary microbiota and duodenum microbiota indicated that gut microbiota showed a higher diversity of genera, both in EG and in CK. At the same time, the speci city and complexity of the two bacterial communities could be seen from as many as 60 differential genera between biliary microbiota and duodenum microbiota in EG. According to the Spearman Association Network Analysis of Dominant genera in EG and the normal group, we found that the relationship between Ralstonia and Ochrobactrum, Streptococcus and Granulicatella, Acinetobacter and Alloprevotella were antagonistic to each other in CK, but they showed cooperative relationships in EG. This shift in bacterial relationship may be the biological basis for the occurrence and development of biliary tract diseases (e g. primary CBD stones). The reasons for this shift are not referenced in the relevant literature, may be related to the elimination or reduction of speci c genus, or may be affected by other genera. These speculations need further research to explore.
The relationship between biliary and intestinal microbial communities has been a subject of continuous research and discussion. Bile bacterial infection was an important factor in the formation of CBD stones. Based on human anatomy and related studies, biliary microbiota in patients with CBD stones originated from the intestinal tract [33,34]. In recent years, research has also tried to con rm this view. Ye et al. provided evidence that the biliary ora in patients with CBD stones originated in the intestine and mainly from the upper gastrointestinal tract by the 16S rRNA gene amplicon sequence [4].In this study, the beta diversity analysis of the sample community structure showed that the biliary and the duodenum microbiota showed some similarities in both EG and CK, especially in terms of bacterial population abundance (slightly different in bacterial composition), moreover this similarity was more evident in the EG. At the same time, a Venn diagram showed "core microbiome" (1205 shared OTUs), which was rst proposed by Turnbaugh et al [35].Thus we can conclude that the most of biliary microbiota and the duodenal microbiota have a common source. More importantly, this study con rmed the existence of self-microbiota in the normal biliary tract, which constituted a microenvironment of bile together with bile salts, bile acids, and so on. In addition, though the abundance of unique OTUs was low in four subgroups, we still have the possibility to nd biomarkers related CBD stones among them.
In order to gain more understanding of the biliary and duodenal microbiota, we performed microbial function prediction by PICTUSt. Both the duodenal and the biliary microbiota have gained rich function pathways. There are differences in some pathways at KEGG level 3 in bile samples between EG and CK, such as cellular processes (Apoptosis, Endocytosis), environmental information processing (Twocomponent system, CAM ligands) and the immune system (RIG-I-like receptor signaling pathway), of which the abundance in EG decreased. We believe that some genera appeared to have a decline of selfrepair ability and environmental adaptability in bile in EG and were in a "pathological state." In duodenal uid samples, microbiota in EG hadhigher abundance of functions involved in Naphthalene degradation and Bacterial invasion of epithelial cells. Studies have shown that Pseudomonas sp. JLR11 and E. coli were related to pathway of Naphthalene degradation [36,37].Some bacteria such as Shigella, Streptococcus and Listeria could enter the epithelial cells through pathway of Bacterial invasion of epithelial cells [38][39][40].It is indicated that the increase of bacteria with invasive ability of duodenal microbiota may indirectly lead to the formation of the stone, which may also explain the increase in the abundance of Escherichia-Shigella in the duodenal uid of EG. However, in the Infectious Diseases pathway, Shigellosis pathways and pathogenic Escherichia coli infection pathways had no difference in duodenal microbiota between two groups. This is also a good illustration of why the high abundance of Escherichia-Shigellain duodenal ora does not cause related diseases.
This study further explored the composition of the biliary and intestinal microbiota in patients with primary CBD stones and the differences in key genera. And we had a completely new understanding of the biological communities in the normal biliary tract. However, this experiment still has its limitations.
The small sample size may lead to certain deviations. As far as research methods are concerned, in the future, metagenomic sequencing, animal experiments and other methods may be needed to study and validate the biliary ora, so as to further clarify the role of bacteria on the formation of primary CBD stones.

Conclusion
In conclusion, our results suggest the existence of the intrinsic ora in normal bile and give a more comprehensive disclosure of its composition. The disorder of the biliary ora is of great importance in the formation of CBD stones, which is mainly re ected in the reduction of some dominant bacterium and the change of the association network of dominant genera. These are new understandings of the biliary microbiota and pathogenic mechanisms. Equally important is the discovery of reduction of bacteria which could produce short chain fatty acids (SCFA), the decrease may affect the health of the biliary tract and lead to the formation of stones. Maybe this discovery will guide us in the prevention and treatment of primary CBD stones in the future.

Study Design and Sample Collection
We The bile (B group and BCK group) and duodenal juice (D group and DCK group) were respectively extracted from two groups of patients under ERCP.
The entire ERCP process was strictly aseptic and all experimental materials were sterilized according to the procedure. During the entry of the duodenoscope (TJF240/JF260V; Olympus), any inhalation was strictly prohibited to avoid contamination. Before duodenal uid was drawn, the duodenoscope channel was washed with 10 ml acidi ed water (PH = 2.5±0.2) and 20 ml physiological saline at the level of the duodenum, and used balloon extraction catheter to draw at least 2 ml of intestinal uid within 3 cm of the duodenal papilla. Then, the channels and papilla were ushed again with 20 ml of normal saline, and 5-10 ml of bile was extracted from the common bile duct using a Tripe Lumen Sphincterotome. The samples were placed in sterile sputum cup and immediately stored in the − 80℃ refrigerator for later bacterial DNA extraction and high throughput sequencing analysis.
All patients in this study expressed their understanding and signed informed consent and obtained approval from the Ethics Committee of the Second Hospital of Hebei Medical University.

16S rRNA gene ampli cation and sequencing of amplicons
The genomic DNA of bile and duodenal uid samples were extracted by means of Soil DNA Kit (Omega).
The obtained DNA was subjected to 0.8% agarose gel electrophoresis for molecular size determination, and the DNA was quanti ed by ultraviolet spectrophotometer. The DNA is stored in the − 20℃ refrigerator for subsequent analysis after passing the test. In order to ensure the quality of sequencing, the insertion fragment range of the best sequencing was 200-450 bp by Miseq sequence reading length. In this experiment, the V3 ~ V4 region of 16SrRNA gene was selected for sequencing. The PCR primers used were 515F (5'-barcode-GTGCCAGCMGCCGCGG-3') and 907R (5'-CCGTCAATTCMTTTRAGTTT-3'). Conventional PCR thermal cycling program: Initial denaturation at 98°C for 30 s, Denaturation at 98°C for 15 s, Annealing at 50 °C for 30 s, Extension at 72°C for 30 s, The above three processes were carried out for 25 cycles and lastly extended for 5 min at 72°C, nally stored at 4-10°C conventionally. The ampli cation results were performed on 2% agarose gel electrophoresis. After passing the test, the target fragment was excised and the target fragment was recovered using the Axygen Gel Recovery Kit. The PCR products were quanti ed on a Microplate reader (BioTek, FLx800) using the Quant-iTPicoGreen dsDNA Assay Kit, and then mixed according to the amount of data required for each sample. The obtained PCR products above were used to construct a sequencing library by using Illumina trusEq. The Quali ed sequencing library (index sequence cannot be repeated) was diluted according to the gradient, then mixed according to the required amount of sequencing, and denatured by NaOH into single strands for on-machine sequencing. Barcoded V3-V4 amplicons were sequenced using the pair-end method by Illumina Miseq [41]. The optimal sequencing length was chosen to be 200-450 bp. Using FLASH software [42] to pair pair-end sequences based on the initial screening according to the overlapping bases: the overlapped bases of both Read 1 and Read 2 sequences are required to be≥10 bp in length, and base mismatches are not allowed. Afterwards, the valid sequences were obtained according to the Index information corresponding to each sample. Finally, use QIIME software (Quantitative Insights Into Microbial Ecology, v1.8.0) [43] to identify the error sequence and remove it. Ampli cation and sequencing of 16S rRNA gene was completed by Personal Biotechnology Co., Ltd. (Shanghai, China).

OTU clustering and statistical analysis
Operational Taxonomic Unit(OTU) [44]usually refers to a set of human-de ned sequence similarity thresholds. The sequences from one or more samples are merged, and the sequences whose similarity is higher than the threshold will be merged into one OTU. Using the UCLUST sequence comparison tool in QIIME software [45], the sequences obtained above were merged and OTU-divided with a sequence similarity of 97%, and selected the highest abundance sequence in each OTU as the representative sequence of the OTU. The OTU with an abundance value less than 0.001% of the all sample total sequencing was removed [46], then used Silva database as the template sequence for the identi cation of the OTU taxonomic status.
Statistical analyses were performed using the SPSS 23.0 software. The Fisher's exact test was used to enumeration data, and the Kruskal-wall is test to measurement data. The OTU was classi ed by QIIME software using the Silva database and the bacterial annotation was carried out. Differences within or between groups were compared by Alpha diversity, Beta diversity, Metastats analysis, PLS-DA analysis,