Gallbladder microbiota composition is associated with pancreaticobiliary and gallbladder cancer prognosis

The microbial population of the intestinal tract and its relationship to specific diseases has been extensively studied during the past decade. However, reports characterizing the bile microbiota are rare. This study aims to investigate the microbiota composition in patients with pancreaticobiliary cancers and benign diseases by 16S rRNA gene amplicon sequencing and to evaluate its potential value as a biomarker for the cancer of the bile duct, pancreas, and gallbladder. We enrolled patients who were diagnosed with cancer, cystic lesions, and inflammation of the pancreaticobiliary tract. The study cohort comprised 244 patients. We extracted microbiome-derived DNA from the bile juice in surgically resected gallbladders. The microbiome composition was not significantly different according to lesion position and cancer type in terms of alpha and beta diversity. We found a significant difference in the relative abundance of Campylobacter, Citrobacter, Leptotrichia, Enterobacter, Hungatella, Mycolicibacterium, Phyllobacterium and Sphingomonas between patients without and with lymph node metastasis. There was a significant association between the relative abundance of certain microbes and overall survival prognosis. These microbes showed association with good prognosis in cholangiocarcinoma, but with poor prognosis in pancreatic adenocarcinoma, and vice versa. Our findings suggest that pancreaticobiliary tract cancer patients have an altered microbiome composition, which might be a biomarker for distinguishing malignancy.

stage due to the lack of early signs and clinical symptoms, as is pancreaticobiliary tract cancer, which limits the selection of therapy and undermines a better prognosis. Therefore, it is crucial to identify an effective biomarker to enable early diagnosis and predict the prognosis of gallbladder and pancreaticobiliary cancer.
The human body is colonized by over 100 trillion symbiotic microorganisms (almost equivalent to the number of cells in a human) and collectively referred to as the microbiota [9,10]. Due to environmental differences, each site in the body is home to distinct microbial ecosystems [10]. Of them, the most diverse bacterial populations occur in the intestinal tract [11,12]. The human gut microbiota contributes to host physiologic development and maintenance, including education of the host immune system, nutrient digestion, and defense against colonization by pathogenic microorganisms [13,14]. The gut microbiota is increasingly considered an important factor associated with both tumor development and the efficacy of anticancer therapies [15,16].
Bile juice was considered sterile due to the difficulties in accessing biological samples, but recent reports indicate the existence of a microbial ecosystem in people with and without hepatobiliary disorders [17,18]. Furthermore, a clinical study using metagenomic analysis showed the association between carcinogenesis with liver flukes and microbiota in biliary tract cancers [19]. Other studies have shown that intrapancreatic microbiota may mediate tumor resistance to gemcitabine [20]. Hence, a better understanding of the roles of microbes in the development of hepatobiliary-pancreatic tumors may reveal opportunities to develop new prevention and treatment strategies for patients with hepatobiliary-pancreatic cancers by targeting microbes and the microbiota.
In this study, we performed 16S rRNA gene amplicon sequencing analysis of bile juice collected from resected gallbladders in cases of pancreaticobiliary tract and gallbladder cancers to investigate whether alterations in microbiota composition in the gallbladder affect the patient's prognosis after surgery. We anticipated that the composition of individual microbiomes might be a novel biomarker to predict the prognosis of pancreaticobiliary tract and gallbladder cancers.

Differences of microbiome composition in the gallbladder
We isolated the bacterial DNA derived from the bile juice in resected gallbladder samples with pancreaticobiliary tract cancers, gallbladder cancers, pancreas cystic lesions, other cancers, and benign inflammatory lesions. Then, the variable regions (V3-V4) of the 16S rRNA genes were amplified by polymerase chain reaction (PCR). The number of 16S rRNA sequences per bile sample ranged from 10,254 to 342,362. We identified 11,358 ASV (Amplicon Sequence Variants) by subsequent DNA sequencing analysis. Of them, we assigned 19 ASV at the phylum level, 28 ASV at the class level, 61 ASV at the order level, 122 ASV at the family level, and 262 ASV at the genus level (Fig. 1A). Then we assigned ASV at the genus level using BLAST searches (Supplementary Table 1) [21]. There were no differences in the alpha diversity among the lesion locations (Pielou evenness index: p = 0.431; Faith PD: p = 0.703 and Chao1: p = 0.403) without Shannon index (p = 0.024) or lesion type (Pielou evenness index: p = 0.902; Faith PD: p = 0.853; Chao1: p = 0.403 and Shannon index: p = 0.131) (Fig. 1B, C). Beta diversity of the biliary microbiome was also compared. Non-metric multidimensional scaling (NMDS) of the centered log-ratio-transformed data did not show any distinct clustering, indicating the absence of overall microbiome differences among the types and locations of lesions ( Fig. 2A). Similarly, principal coordinate analysis (PCoA) did not show distinct clustering, indicating the absence of overall microbiome differences among the types and locations of lesions (Fig. 2B). However, there was a significant association between trends in these variances and the N-score of the TNM staging system (Supplementary Table 2). Unsupervised hierarchical clustering analysis was performed for relative abundance of microbiota data sets including bile duct lesion, pancreas lesion and other lesion (Fig. 3A). The samples were divided into two clusters according to the clustering results. In bile duct lesion, Cluster 1 showed a significantly better prognosis than Cluster 2 (HR = 0.195, p = 0.015; Fig. 3B). However, in pancreas lesion, it has no significant difference prognosis between Cluster 1 and Cluster 2 (HR = 1.032, p = 0.921; Fig. 3C). The relative abundance of Bradyrhizobium, Carnobacterium, Cutibacterium, Enterococcus, Fusobacterium, Methylobacterium, Phyllobacterium, Pseudomonas, Serratia and Streptococcus were significant difference between the cluster 1 and cluster 2 (Fig. 4).

Association between microbial abundance and clinical features
The statistical analysis for the association between microbial individual relative abundance and clinical features is summarized in Table 2. Females showed higher Escherichia and Streptococcus abundance than males (p = 0.007 and p = 0.030, respectively). There were significant differences in Sphingomonas and Fusobacterium between those participants aged over and under 70 years (p = 0.011 and p = 0.044, respectively). In the cholangiocarcinoma, Campylobacter, Citrobacter and Leptotrichia abundance showed significant increase between Clinical stage with and without lymphnode metastasis (p = 0.025, p = 0.001 and p = 0.007, respectively). In the pancreatic ductal adenocarcinoma, Enterobacter, Hungatella, Mycolicibacterium, Phyllobacterium and Sphingomonas showed significant difference between with and without lymphonude metastasis (p = 0.004, p = 0.007, p = 0.018, p = 0.023 and p = 0.058, respectively). Additionally, there were significant differences in the relative abundance of Schaalia, Alloprevotella, Bilophila, Dialister, Eggerthella, Selenomonas and Streptococcus between Intraductal papillary mucinous carcinoma (IPMC) and intraductal papillary mucinous neoplasm (IPMN) (p = 0.038, p = 0.047, p = 0.047, p = 0.047, p = 0.047, p = 0.047 and p = 0.009, respectively). It showed no significant differences in abundance of microbiota to treatment of chemotherapy and radiotherapy before surgery (data not shown).

Association between microbial abundance and prognosis
To investigate whether the individual relative abundance of microbiota was associated with prognosis, we performed univariate and/or multivariate Cox regression analysis for overall survival ( Table 3). The relative abundance of Abiotrophia, Amaricoccus, Blastococcus, Bosea, Delftia, Dokdonella, Flavobacterium, Gemella, Granulicatella, Haemophilus, Leucobacter, Pelagibacterium, Sphingopyxis, Streptococcus and Williamsia were significantly correlated with prognosis on univariate and multivariate Cox regression analysis after adjustment for clinicopathologic variables, such as ASA score, age, sex, and preoperative chemotherapy (Table 3-1). In the bile duct lesions, the relative abundance of Delftia, Dermacoccus, Haemophilus, Leucobacter, Methylocapsa and Staphylococcus was significantly associated with prognosis on univariate and multivariate Cox regression analysis after adjustment for clinicopathologic variables, such as ASA score and the presence of lymph node metastasis (Table 3-2). In the pancreas lesion, the relative abundance of Abiotrophia, Aureimonas, Flavobacterium, Gemella, Howardella, Klebsiella, Proteus, Pelagibacterium, Sphingopyxis and Williamsia showed a significant correlation with prognosis on univariate and multivariate Cox regression analysis after adjustment for clinicopathologic variables, such as ASA Score, Age, Sex, the presence of lymph node metastasis, and preoperative chemotherapy (Table 3-3). There is no common microbe correlate with prognosis between the bile duct lesion and the pancreatic lesion.

Evaluation of threshold value of individual microbiota relative abundance for prognosis
We evaluated the threshold value of individual microbiota relative abundance for predicting prognosis using Cox proportional hazards model analysis (summarized in Table 4). In total samples, these groups with a high relative abundance of Enterococcus, Eggerthella, Klebsiella, Corynebacterium, Moraxella, Hungatella, Paracoccus, Dermacoccus, Citrobacter, Lawsonella and Pseudoxanthomonas showed a significantly poor prognosis compared with the other group (HR = 1.65, HR = 2.22, HR = 2.21, HR = 2.36, HR = 2.27, HR = 2.74, HR = 2.50, HR = 3.14, HR = 2.60, HR = 3.48 and HR = 7.41, respectively). On the other hand, these groups with a high relative abundance of Streptococcus, Escherichia, Veillonella and Dialister showed a significantly better prognosis compared with the other group (HR = 0.60, HR = 0.59, HR = 0.50 and HR = 0.35, respectively). In Eggerthella and Corynebacterium, the Fisher exact test revealed that there was significant difference in the Sex between these two groups divided by these threshold values (p = 0.026 and p = 0.046, respectively. Supplementary Table 3). In Staphylococcus and Veillonella, the Fisher exact test revealed that there was significant difference in the Age between these two groups divided by these threshold values (p = 0.021 and p = 0.003, respectively. Supplementary Table 3).
In the pancreatic lesions, the group with high relative abundance of Klebsiella, Veillonella, Acinetobacter, Selenomonas and Paracoccus showed a significantly poor prognosis compared with the other group (HR = 1.82, HR = 2.87, HR = 2.88, HR = 3.49 and HR = 4.04, respectively). On the other hand, these groups with a high relative abundance of Enterococcus, Staphylococcus, Bacteroides, Raoultella and Streptococcus showed a significantly better prognosis compared with the other group (HR = 0.37, HR = 0.28, HR = 0.37, HR = 0.12 and HR = 0.55, respectively). The Fisher exact test revealed that there was significant difference in the ASA score between these two groups divided by these threshold values in Raoultella.

Discussion
Pancreaticobiliary tract cancers, such as PDAC, cholangiocarcinoma, and gallbladder carcinoma, are aggressive malignancies with a high risk of invasion and metastasis. Furthermore, they are resistant to most cytotoxic agents [4,5,22,23]. Due to the lack of sensitive clinical methods to detect these pancreaticobiliary tract cancers, most chemotherapeutic patients are often diagnosed at advanced stages and show an abysmal prognosis [24][25][26]. Thus, it is critical to establish new diagnostic, prognostic, and therapeutic biomarkers. Recent studies of the microbiota in the colorectum have suggested numerous links between these microbial communities and colorectal cancers [13,16]. However, the association of changes in the microbiota in the gallbladder with pancreaticobiliary tract and gallbladder cancer has been rarely reported. In this study, we find a significant association between the relative abundance of certain microbes in gallbladder and the malignancy of lesion. Then, these microbes showed association with good prognosis in cholangiocarcinoma, but with poor prognosis in pancreatic adenocarcinoma, and vice versa.
The microbiota in the normal gallbladder consists of five main phyla: Proteobacteria, Firmicutes, Bacteroidetes, Fusobacteria, and Actinobacteria [18]. We also found an equivalent relative abundance of phyla in the bile juice derived from gallbladders in our study with pancreaticobiliary and gallbladder lesions. NMDS and PCoA analysis did not show distinct clustering, indicating the absence of overall microbiome differences among the types and locations of lesions. However, there was a significant association between trends in these variances and the N-score of the TNM staging system. Moreover, clustering analysis using relative abundance showed 2 cluster have significantly different prognosis in bile duct lesion. These results suggested that microbiota in gallbladder showed association with pancreaticobiliary and gallbladder cancer malignancy. In the cholangiocarcinoma, genes relative abundance analysis and clinical information showed that Canpylobactor, Citrobacter and Leptotrichia increased as the progress of lymphonode metastasis. Similarly, in the pancreatic adenocarcinoma, Phyllobacterium and Sphingomonas increased as the progress of lymphonode metastasis. In contrast, Enterobacter, Hungatella and Mycolicibacterium decreased. Our results support these findings of a wide range of infectious etiologies caused by Campylobactor, Leptotrichia, Phyllobacterium and Sphingomonas at different anatomic sites have been reported in the literature, suggesting its highly pathogenic potential [27][28][29]. Collectively, these results suggest that the progress of pancreaticobiliary tract cancer may affect changes in the gallbladder environment and, in turn, microbiome composition. In particularly, Schaalia, Alloprevotella, Bilophila, Dialister, Eggerthella, Selenomonas and Streptococcus showed a higher relative abundance in invasive intraductal papillary mucinous carcinoma than in intraductal papillary mucinous neoplasms. Thus, the pancreatic lesions affect the gallbladder microbiota, despite the pancreas being anatomically distant from the gallbladder.
The presence of Delftia, Dermacoccus, Haemophilus, Leucobacter, Methylocapsa and Staphylococcus were prognostic factors after adjustment for clinicopathologic variables in bile duct lesion. Delftia, Haemophilus and Leucobacter were common factor between the total Fig. 4 Comparison relative abundance of microbiota between Cluster 1 and Cluster 2. Comparison relative abundance of microbiota between the Cluster 1 and the Cluster 2 on clustering analysis. The color of plot means location of lesion. The red plot, bile duct lesion; the blue plot, gallbladder lesion; the green plot, pancreas lesion; the black plot, other. The plot shape means type of lesion. ■, cancer; ▲, cystic lesion; •, benign; ×, Other cases and bile duct lesion cases. On the other hand, in pancreatic lesion analysis, Abiotrophia, Aureimonas, Flavobacterium, Gemella, Howardella, Klebsiella, Proteus, Pelagibacterium, Sphingopyxis and Williamsia were prognostic factors after adjustment for clinicopathologic variables. Abiotrophia, Flavobacterium, Gemella, Pelagibacterium, Sphingopyxis and Williamsia were common factor between the total cases and pancreatic lesion cases. Naito et.al., indicates that using human enteroids derived from the transverse colon, lipopolysaccharide from crypt-specific core microbiota (i.e., Acinetobacter, Delftia, and Stenotrophomonas) induced an increase in goblet cell-associated proteins such as MUC2 [30]. Unusual expression of mucin in pancreaticobiliary cancer may be responsible for mucosal microbiota. Genus-level analyses showed that four genera (Actinomyces, Atopobium, Fusobacterium, and Haemophilus) were present in significantly high proportions in colorectal cancer [31]. Tumor and the peri-tumoral regions of prostate had a higher relative abundance of Staphylococcus compared to normal areas [32]. There were some patients with cholangiocarcinoma and pancreatic adenocarcinoma that had a very high abundance of these genus, suggesting their direct involvement in malignant progression.
Furthermore, we found that the threshold value of relative abundance of microbiota is a significant marker for the prognosis of pancreaticobiliary tract cancer. These threshold values and clinical conditions, such as sex, age, American Society of Anesthesiologists (ASA) score, stages, and the administration of preoperative chemotherapy were non-confounding factors. In bile duct cancer, Enterococcus, Staphylococcus and Bacteroides were poor prognosis factors, but they were good prognosis factors for pancreatic cancer. Streptococcus was only common good prognosis factor. These results showed no common microbe correlate with poor prognosis between the bile duct lesion and the pancreatic lesion. Thus, the effect of the gallbladder environment was different for each lesion and that the difference in prognosis among the lesions might be caused by the gallbladder

Conclusions
In conclusion, this is the first report to indicate a link between gallbladder microbiota and pancreaticobiliary cancer prognosis. Although we took precautions in collecting biliary fluid from surgically resected gallbladders to prevent contamination by gastric or duodenal juices, it is impossible to completely exclude the likelihood of contamination because fluid samples from the stomach, duodenum, and intestines were not cross-checked in this study. However, our findings demonstrate that the alterations in the gallbladder microbiota population could be used to accurately distinguish the overall survival prognosis in pancreaticobiliary tract cancer patients after surgery.
Further studies of these microbial markers are necessary to facilitate patient counseling, decision making regarding individualized therapy, and follow-up scheduling.

Bile samples and patients
We obtained 244 bile juice samples from surgically resected gallbladders. These bile juice samples were obtained from the gallbladders by a pathologist using a needle and immediately frozen (− 70 °C     . After PCR products were quantified, equimolar ratios from each sample were pooled and sequenced on a MiSeq System (Illumina) according to the manufacturer's instructions. It showed no amplification that the negative control using PBS buffer apply to similarly steps for DNA extraction and several PCR.

Taxonomic assignment
Raw reads obtained from the sequencer were filtered according to the barcode and primer sequences using the MiSeq system. Then, the reads were imported into QIIME2 [34] v2019.4 in Linux, and quality assessment, filtering, and chimera detection were performed using the DADA2 pipeline. Taxonomic classification was assigned to amplicon sequence variants using 99% clustering in SILVA 132 (https:// www. arb-silva. de/ downl oad/ archi ve/) [35]. All samples were rarefied to the lowest reads, i.e., 10,000 reads, to minimize the effects of sequencing depth on alpha and beta diversity measures using "qiime feature-table rarefy" in Qiime2. The amplicon sequence variants were adequately detected for relative abundance (permyriad of total sequences) and alpha-diversity determination (i.e., Shannon Index, Cho1, Faith phylogenetic diversity [PD], and Pielou evenness index).

Statistical analysis
Data were analyzed using the R computing environment v4.0.2 [36]. The normality of the data distribution was evaluated using the Kolmogorov-Smirnov test. Differences between groups were analyzed using the Welch t-test. A nonparametric test of group differences was performed using the Mann-Whitney U test (Wilcoxon rank sum test) and Kruskal-Wallis test. A parametric test of group differences was performed using the student t test. The effect size was calculated on cliff 's delta as nonparametrical and cohen's d as parametrical. Categorical variables were compared by the Fisher exact test. Survival was estimated by the Kaplan-Meier method.
The Cox proportional hazards model and survival curves were compared using the log-rank test. Prognostic factors were adjusted by a Cox regression model. All of prognosis analysis were performed on overall survival. The output data from QIIME2 were mined to determine alpha and beta diversity analyses and composition analysis in R using qiime2R, phyloseq, effesize and tidyverse packages. The clustering analysis on vegan package. A p-value < 0.05 was considered statistically significant.