Gluconobacter dominates the gut microbiome of the Asian palm civet Paradoxurus hermaphroditus that produces kopi luwak

Collection of fecal samples and DNA extraction

Three fresh fecal samples from wild civets were collected in the morning on June 20, 2015 from a coffee plantation (GPS coordinate: latitude: 1.7495934, longitude: 103.387647) in Johor, Malaysia. A field permit was not required since the raw kopi luwak specimens were provided by the owner of the coffee plantation, Mr. Jason Liew. The samples were kept at 4 °C for ≤2 h before being used for DNA extraction. To improve the DNA yield during extraction, ~0.2 g of each fecal sample was first incubated in 200 μl of suspension buffer (10 mM Tris-HCl, pH 8.0, 0.1 mM EDTA) containing 2 mg/ml lysozyme (Sigma Aldrich, St. Louis, MO, USA) at 37 °C for 30 min to remove the bacterial cell wall. The DNA was extracted from feces using the QIAmp^® DNA stool kit (Qiagen, Hilden, Germany). To digest RNA, RNase A (Sigma Aldrich, St. Louis, MO, USA) was added (final concentration of 0.1 mg/ml) to buffer ASL from the stool kit. The quantity and quality of genomic DNA were estimated using a NanoDrop^™ 2000/2000c spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA).

PCR amplification, sequencing and analysis of the 16S rRNA gene

A fragment of the 16S rRNA gene was amplified using the prokaryotic universal primer set 27F (5′-AGAGTTTGATCCTGGCTCAG-3′) and 338R (5′-TGCTGCCTCCCGTAGGAGT-3′) and Ex Taq DNA polymerase, hot-start version (Takara Bio, Shiga, Japan). PCR was carried out in 50-μl reaction volumes containing forward and reverse primers (1 μM each), 20 μl of template DNA, 5 μl of 10× Ex Taq Buffer, and 4 μl of dNTP mix from the Ex Taq kit. The PCR program consisted of 95 °C for 5 min followed by 31 cycles of 95 °C for 30 s, 53 °C for 30 s and 72 °C for 30 s, and a final extension at 72 °C for 3 min. The PCR products were individually concentrated and purified using a 2% E-Gel SizeSelect agarose gel (Thermo Fisher Scientific, Waltham, MA, USA), quantified using the Quant-iT dsDNA HS Assay kit (Thermo Fisher Scientific, Waltham, MA, USA), and qualified using the High Sensitivity DNA kit (Agilent Technologies, Santa Clara, CA, USA). Then, we checked a negative-control sample, in which PCR produced no visible bands on an agarose gel. Sequencing was carried out with the Ion PGM Sequencing 400 kit (Thermo Fisher Scientific, Waltham, MA, USA). Raw data for 16S rRNA gene sequences were analyzed with the VITCOMIC2 web application (http://vitcomic.org/) to derive the genus composition with the option “Conduct 16S rRNA gene Copy number normalization?—No”. VITCOMIC2 can estimate microbial community composition based on the sequence data for the 16S rRNA gene obtained from both metagenomic shotgun and amplicon sequencing (Mori et al., 2018). The sequence data are available in the DDBJ DRA database (https://www.ddbj.nig.ac.jp/dra/index-e.html) under accession DRA006640.

16S rRNA gene sequence–based phylogenetic analysis

The 16S rRNA gene sequences from bacteria in civet cat feces were checked using Trimmomatic version 0.33 with parameters “SE LEADING:17 TRAILING:17 AVGQUAL:25 MINLEN:200” to filter out low-quality sequences (Bolger, Lohse & Usadel, 2014). The three most abundant Gluconobacter sequences (termed civet cat004:Uniq1, civet cat005:Uniq1, civet cat006:Uniq1) were retrieved from each sample using the fastx unique command implemented in USEARCH version 10.0.240 (Edgar, 2010). To compare our Gluconobacter 16S rRNA gene sequences with those of other Gluconobacter species, we collected 16S rRNA gene sequences of Gluconobacter and Acetobacter from the SILVA Living Tree Project version 128 (Yilmaz et al., 2014). The sequences were aligned with MAFFT version 7.313 (Nakamura et al., 2018) and used as input to construct a phylogenetic tree using the neighbor-joining method with a Kimura 2-parameter model of nucleotide substitution using MEGA7 (Kumar, Stecher & Tamura, 2016). Acetobacter was used as an outgroup as the most similar taxonomic group to Gluconobacter species.

Collection of animal fecal 16S rRNA gene sequence data from other studies

Data for the bacterial genera identified in feces of various animals were collected from MicrobeDB.jp (MicrobeDB.jp Project Team, 2017). MicrobeDB.jp is a public database that stores the taxonomic composition data computed by VITCOMIC2 for most 16S rRNA gene sequences in the International Nucleotide Sequence Database Collaboration-Sequence Read Archive (INSDC-SRA), with integrated sampling site information such as environment and host categories. After downloading all the taxonomic composition data, all data with alignment counts less than 1,000 were filtered out. We listed the number of 16S rRNA gene sequence samples used in a comparative analysis after quality control of the raw sequence data (Table 1). Next, we converted alignment counts to mean relative abundance for each sample. After selecting samples that had an annotation of feces (MEO_0000054) by Metagenome and Microbes Environmental Ontology (MEO), we grouped samples into host categories and calculated their average relative abundances to serve as a reference microbiome composition for the gut of all animals in the analysis (MicrobeDB.jp Project Team, 2017).

Comparative genomics of the major microorganisms in feces from civet cat and other animals

To provide a plausible explanation for the known metabolite profiles of kopi luwak, we compared the major bacterial species in civet cat feces (G. frateurii and G. japonicus) with the top nine dominant bacterial genera in feces from other animals (Bacteroides, Prevotella, Barnesiella, Lactobacillus, Oscillibacter, Citrobacter, Streptococcus, Faecalibacterium and Enterococcus). We downloaded complete amino acid sequences for each bacterial species from the National Center for Biotechnology Information (NCBI) RefSeq collection (Table S1) and annotated them with the reference prokaryote Kyoto Encyclopedia of Genes and Genomes (KEGG) orthology (KO) amino acid sequences (KEGG FTP Release 82.0) using DIAMOND version 0.9.3 with the parameter “blastp—max-target-seqs 5—e value 1E−2—id 70—min-score 40” (Buchfink, Xie & Huson, 2015). We used the top-hit results with e-values less than 1E−8 as the KO annotation for the amino acid sequences.

To compare KO terms, we calculated the median number of paralogs for each KO in each bacterial genus, including Gluconobacter from the NCBI RefSeq collection (Table S1). For the Gluconobacter genus, we used G. frateurii and G. japonicus because these two species are most closely related to the Gluconobacter species of the civet cat fecal microbiome (Table S2). To determine the KO(s) that was over-represented in Gluconobacter but not in major fecal microbial genera of other animals, we calculated the $K{O}_{\left(i,ratio\right)}$ as follows.

Here, the $K{O}_{\left(i,Gluconobacter\right)}$ is calculated as the median for the number of $K{O}_i$ homologs in each Gluconobacter strain. Conversely, $K{O}_{\left(i,majorfecalmicrobialgenera\right)}$ is the mean of the median for the the number of $K{O}_i$ homologs in the top nine major bacterial genera in other animal feces (see Table S2). To estimate the over-represented functional units (KEGG module/KEGG pathway) in Gluconobacter (G. frateurii and G. japonicus), we performed Fisher’s exact test for each functional unit based on the number of over-represented KOs belonging to each functional unit in Gluconobacter. The $K{O}_{rati{o}_i}$ was converted to the log₂ fold-change. We calculated (A) the number of KOs in a specific functional unit i, (B) the number of KOs out of the functional unit i, (C) the number of KOs with a log₂ fold-change greater than 1 in specific functional unit i and (D) the number of KOs with a log₂ fold-change greater than 1 in functional unit i. These four types of numbers for KO were used to perform Fisher’s exact test with R version 3.3.2 (R Core Team, 2019; Fisher, 1962).

For a P value ≤ 0.05 with Fisher’s exact test, the generated functional unit i can be considered a special functional unit of G. frateurii and G. japonicus.

To identify functional units in G. frateurii and G. japonicus with log₂ fold-changes > 1 relative to the fecal microbes from other animals, we mapped each functional unit (P value ≤ 0.05) and the positive ratio of G. frateurii and G. japonicus for each KO using the FuncTree2 web service (https://bioviz.tokyo/functree2/) (Uchiyama et al., 2015). This web service is a hierarchical visualization tool for KEGG pathways/modules/KO.

Comparison of fecal microbial communities between civet cat and other animals

We analyzed the 16S rRNA gene sequences of microbes in each of three feces samples obtained from three different civet cats, which yielded 28,120,904 raw reads. Among all samples, the five most average relative abundant genera were Gluconobacter (66%), Citrobacter (14%), Acetobacter (7%), Enterobacter (2%) and Clostridium XI (1%) (Table S3). To determine the bacterial genus that is, predominant in the civet cat fecal microbiome, we compared the civet cat fecal 16S rRNA gene sequencing results with those obtained from the fecal microbiome of other animals. The data for the genus composition of the fecal microbiome of 14 animals were obtained from MicrobeDB.jp (MicrobeDB.jp Project Team, 2017) (http://microbedb.jp). Microbial composition data derived from the same computational pipeline were used for the civet cat fecal microbiome 16S rRNA gene sequence analysis. Among all fecal samples from all animals, the largest sample size was that of human feces (n = 3,003), with the second largest being mouse feces (n = 1,346). Most of the human feces samples have been acquired previously for large-scale projects such as the Human Microbiome Project (Turnbaugh et al., 2007) and the Metagenomics of the Human Intestinal Tract project (Ehrlich, 2011) (Table 1).

We determined the microbial genus composition of feces for each animal species by averaging the relative abundance among all fecal samples for each animal. A comparison of these derived compositions with that determined for civet cat feces revealed a marked difference in bacterial genus composition (Fig. 1; Table S3). For example, Gluconobacter was found to be the predominant genus, representing 60% of the bacterial abundance in all civet cat samples. By contrast, Gluconobacter was rarely detected in the fecal microbiomes of the 14 animals used in the other studies. Therefore, civet cat apparently has a unique fecal microbiome compared with other animals.

Phylogenetic analysis of Gluconobacter species in civet cat feces

To determine the Gluconobacter species that are most common in the civet cat fecal microbiome, we created a phylogenetic tree for Gluconobacter species and Acetobacter species as the outgroup in civet cat fecal samples (Fig. 2). The tree showed that Gluconobacter comprises two groups, namely group 1 and group 2. The species in civet cat fecal samples belonged to group 2 (Fig. 2). The tree also showed that the Gluconobacter species in the civet cat fecal microbiome belonged to the group including G. cerinus, G. frateurii, G. japonicus, G. nephelii, G. thailandicus and G. wancherniae (Fig. 2; Fig. S1). To determine the KO terms for Gluconobacter species, we calculated the KO terms for the amino acid sequences of 39 Gluconobacter strains from the Reference Sequence (RefSeq) database of NCBI. We found that RefSeq genomes for G. nephelii and G. wancherniae were not available in NCBI; hence, these two species were not included in the KO analysis. Next, using hierarchical clustering of the KO terms for each Gluconobacter strain, we showed that strains with similar 16S rRNA gene sequences clustered together, with distinct separations among the clades (Fig. 3). These results demonstrated that the Gluconobacter species of the civet cat fecal microbiome were most likely to be G. frateurii and G. japonicus, and that the KO terms of these two species differed from those of other Gluconobacter species, especially species in Gluconobacter group 1.

Comparative genomics of genera of the feces microbiome

We next compared the KO terms for G. frateurii and G. japonicus, which are the most similar species to the Gluconobacter species of the civet cat fecal microbiome, with the KO terms of the dominant genera in feces of other animals. Other dominant bacterial genera in the feces of other animals were Bacteroides, Prevotella, Barnesiella, Lactobacillus, Oscillibacter, Citrobacter, Streptococcus, Faecalibacterium and Enterococcus (Fig. 1). This comparative genomics analysis provided a list of candidate microbial genes that may be involved in the fermentation of kopi luwak in the civet cat digestive tract.

To identify the KO terms of each genus, we downloaded the amino acid sequences of RefSeq genomes for these genera (Table S1). These sequences were mapped to prokaryote KO amino acid sequences that are publicly available in the KEGG. The KEGG database contains information about most of the bacterial genome sequences (Kanehisa, 2000). Finally, we calculated the median number of KOs for each genus. Likewise, we calculated the median number of KOs for G. frateurii and G. japonicus because these two species were the predominant Gluconobacter species of the civet cat fecal microbiome (Fig. 4; Table S2). Interestingly, certain KO terms were more abundant in G. frateurii and G. japonicus than in other Gluconobacter species or were present in only these two species. We next determined the KO terms that were enriched in G. frateurii and G. japonicus compared with other major genera of the various animal fecal microbiomes in our analysis. We calculated log₂ fold-change values, which were derived from the median number of KO terms of G. frateurii and G. japonicus per mean number of KO terms of other dominant bacterial genera in the feces of other animals; this allowed us to determine enrichment, and we counted the number of KO terms with log₂ fold-change greater than 1.0. These results showed that at least 664 KOs of the 3,623 total had a log₂ fold-change value ≥ 1.0 (Table S4). To minimize false-positive results, we performed Fisher’s exact test for each KEGG module and pathway (see “Methods”). The following pathways and modules were found to be more abundant in G. frateurii and G. japonicus than in other animal fecal microbes: porphyrin and chlorophyl metabolism (map00860); cobalamin biosynthesis, cobinamide => cobalamin (M00122); flagellar assembly (map02040); methionine salvage pathway (M00034); histidine biosynthesis, PRPP => histidine (M00026); assimilatory sulfate reduction, sulfate => H₂S (M00176); cysteine and methionine metabolism (map00270) (Fig. 5). Moreover, we analyzed over-represented KO terms of the TCA cycle (map00020) in G. frateurii and G. japonicus because a previous study demonstrated that the levels of malate and citrate produced from the TCA cycle are higher in kopi luwak beans than other coffee beans (Jumhawan et al., 2016). Although Fisher’s exact test could not detect any differences in the KOs for the TCA cycle among different coffee beans (map00020), this test showed that certain KOs of the TCA cycle (map00020) were over-represented in G. frateurii and G. japonicus (Fig. 5; Table S2).