Association of maternal genetics with the gut microbiome and eucalypt diet selection in captive koalas

Ethics statement

This study adhered to the Animals in Research: Reporting On Wildlife (ARROW) guidelines. Sample collection and behavioral recordings were approved by Hirakawa Zoological Park and Awaji Farm Park England Hill as collaborative projects with Hayakawa Lab of Hokkaido University and Ogura Lab of Kitasato University and were fully performed through noninvasive approaches, except for blood collection. To minimize suffering, blood samples were not collected for the purpose of this study; instead, used residues from routine health examinations were employed. Behavioral recordings were completely noninvasively conducted in the zoo-visitors’ area and did not artificially control koalas’ behavior for the purpose of this study. The koalas were healthily kept for the purpose of public exhibitions in the zoos in the enough space of enclosure (>5 m wide, >5 m depth, and >5 m height). Their environments were enriched in terms of their active movement as follows. They were always provided branches of a variety of eucalypt species scattered in the enclosure to increase their active movement to find the eucalypt diet items. Logs were arranged in a three-dimensional manner for them to enjoy moving in any direction. In Japan, different management groups are fed different types of eucalypts according to their preferences, in accordance with husbandry guidelines issued by Japanese Association of Zoos and Aquariums (Japanese Association of Zoos and Aquariums, 2020). The animal experimentation protocol was approved by the President of Kitasato University through the judgment of the Institutional Animal Care and Use Committee of Kitasato University (Approval No. 21-069).

Animals

To perform sufficient statistical analysis, a total of 15 captive koalas were selected and examined in this study (Table 1, Table S1). Of these, nine koalas obtained care at Hirakawa Zoological Park, whereas the other six obtained care at Awaji Farm Park England Hill. Four of the Awaji koalas were southern koalas in the management grouping (Yuki, Daichi, Nozomi, and Midori), whereas the others (Peta and Umi) and all koalas in Hirakawa Zoological Park were northern koalas in the management grouping. We selected these two zoos because Hirakawa Zoological Park has large number of koalas and Awaji Farm Park England Hill is the only zoo that keeps southern koalas in the management grouping in Japan. To be blind test, behavioral recording was performed by MA, CA and RW, and statistical analysis was performed by KK.

Mitochondrial phylogenetic analysis

Mitochondrial phylogenetic analysis, based on the study by Neaves et al. (2016), was conducted to estimate the mitochondrial maternal origin of each koala. Blood or fecal samples of koalas were collected to extract and purify genomic DNA. The collected blood samples were immediately mixed with anticoagulants (EDTA or heparin). All samples were frozen at −20 °C prior to DNA extraction. Total DNA was extracted from blood samples using Qiagen DNeasy Blood and Tissue Kit (Qiagen GmbH, Hilden, Germany) and from fecal samples using QIAamp Fast DNA Stool Mini Kit (Qiagen GmbH). Next, using TaKaRa Ex Taq Hot Start Version (Takara Bio Inc., Shiga, Japan), the mitochondrial DNA control region (D-loop) was amplified via polymerase chain reaction (PCR) using the following primers: MaL15999M (ACC ATC AAC ACC CAA AGC TGA) and MaH16498M (CCT GAA GTA GCA ACC AGT AG) (Fumagalli et al. 1997). The PCR conditions were as follows: initial denaturation (94 °C for 10 min); followed by 35 cycles of denaturation (94 °C for 10 s), annealing (60 °C for 30 s), and extension (72 °C for 60 s); and final extension (72 °C for 10 min). The PCR products were purified via precipitation with isopropanol. Next, the purified PCR products were directly sequenced using PCR primers for complete coverage in both strand orientations via BigDye Terminator v3.1 Cycle Sequencing Kit and 3130 Genetic Analyzer (Applied Biosystems, Bedford, MA, USA). Chromatograms were imported into FinchTV (Geospiza Inc., Seattle, WA, USA) and analyzed.

The phylogeny of the sequenced D-loop region of mitochondrial DNA was analyzed with the sequences of 48 koalas reported by Neaves et al. (2016). A multiple alignment was constructed by MUSCLE (Edgar, 2004). A phylogenetic tree was reconstructed using the maximum-likelihood (ML) method with 1,000 bootstrap resamplings using MEGA11 (Tamura, Stecher & Kumar, 2021).

Foraging data and gut microbiome

Foraging data were collected from six koalas in Hirakawa Zoological Park for 8 days between 21–30 November 2021. Subjects were housed independently (Boonda and Ito) or in cohabitation with two individuals (Himawari with Kibou and Sora with Itsuki). During foraging and feeding, all koala pairs that use the same enclosure rarely interfere with each other. At Hirakawa, four of the five eucalypt species (E. camaldulensis, CR; E. microcorys, M; E. punctata, P; E. robusta, R; E. tereticornis, T) were fed twice a day (9:00 and 16:00). Since Japan is not a natural habitat of eucalypts, available eucalypts for koala diet are limited depending on the cultivation and logistics (Japanese Association of Zoos and Aquariums, 2020). Eucalypt species fed to koalas in this study were regularly fed to koalas as their typical food as available in these zoos. The combination of eucalypts fed at the same time and the frequency of feeding was counterbalanced. The method of behavioral recoding was determined a priori. Eucalypt species consumed by each koala were observed using the instantaneous focal sampling method in 30-s intervals (Altmann, 1974). The observation time was 1 h each at 9:00 am and 4:00 pm immediately after feeding by caretakers and 1 h each at 11:00 am and 2:00 pm outside of the immediate feeding period. Finally, the observation was carried out for 8 days (a total of 32 h per individual).

Fecal samples were collected during foraging observation periods. Fecal sampling was performed by collecting fresh feces immediately after defecation. In case of cohabiting housed individuals, individual identification was performed by direct observation of defecation. The fecal samples of nonsubject koalas were also collected in Hirakawa and Awaji for comparison. Fresh fecal samples were collected for sampling. Finally, 2–5 fecal samples per koala and a total of 62 fecal samples were collected and stored at −20 °C until DNA extraction.

16S rRNA gene sequencing

According to Hayakawa et al. (2018a); Hayakawa et al. (2018b), we performed 16S-based gut microbiome analysis using collected fecal samples. After quantifying the concentration of purified fecal DNA using Qubit dsDNA HS Assay Kit equipped with Qubit fluorometer (Thermo Fisher Scientific, Waltham, MA, USA), we amplified the V3–V4 region of the 16S rRNA gene using KAPA HiFi Hot Start Ready Mix (Kapa Biosystem, Inc., Wilmington, DE, USA). We used the following primer pair: 1S-D-Bact-0341-b-S-17 (forward), CCT ACG GGN GGC WGC AG; S-D-Bact-0785-a-A-21 (reverse), GAC TAC HVG GGT ATC TAA TCC (Klindworth et al., 2013), with the specific overhang adaptors TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG-[3-6-mer Ns]-[forward primer] and GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA G-[3-6mer Ns]-[reverse primer], where 3-6mer Ns can improve the sequencing quality (Lundberg et al., 2013). After confirming PCR amplification via gel electrophoresis, we purified the PCR amplicons with Agencourt AMPure XP beads (Beckman Coulter, Inc, Brea, CA, USA). Then, we performed index PCR using Illumina Nextera XT Index Kit (Illumina, Inc., San Diego, CA, USA). Additionally, we confirmed the presence and appropriate length of index PCR products via electrophoresis and then purified index PCR products with Agencourt AMPure XP beads. All index PCR products were mixed at the same molarities for constructing a library. After PhiX spike-in (30%), we sequenced the library using the Illumina Miseq platform (Illumina, Inc., San Diego, CA, USA) (2 × 301 bp).

Data analysis

The MiSeq base calls were converted to FASTQ files using configureBclToFastq.pl implemented by the bcl2fastq conversion software v1.8.4 (Illumina, Inc., San Diego, CA, USA) (options: no-eamss, mismatches 0, and use-bases-mask Y300n,Y8,Y8,Y300n). The read pairs were demultiplexed, the primer sequences were trimmed, and those with low-quality index sequences were discarded, where the index sequences included nucleotide(s) with a quality score of <30, using clsplitseq in Claident (Tanabe & Toju, 2013) (option: minqualtag = 30).

Quality control and data analysis were performed using QIIME2 v2021.4.0 (Caporaso et al., 2010). To generate amplicon sequence variants (ASVs), DADA2 v2021.4.0 (Callahan et al., 2016) was used to quality filter the sequences with a read cut length of 260 (forward) and 260 (reverse) based on quality control results and denoise chimeric sequences with a read count of 1,000,000 for training the error model. The taxonomy analysis was conducted with the SILVA 138 database (Bokulich et al., 2018; Robeson 2nd et al., 2021). Subsequently, sample depths were rarefied where the value was ≥0.99 for all samples using Good’s coverage (4,105 sequences per sample). Microbial composition of the samples is shown in Fig. S1 at phylum level and at genus level.

We determined alpha diversity (Shannon index, Chao1, Simpson index, Simpson index of evenness, Pielou’s evenness index, Faith’s phylogenetic diversity, and observed features) and beta diversity (unweighted UniFrac, weighted UniFrac (Lozupone & Knight, 2005; Lozupone et al., 2007), Jaccard index, and Bray–Curtis dissimilarity) to analyze differences in the gut microbiome between mitochondrial lineages or between the management groups. The pairwise Kruskal–Wallis test with Benjamini–Hochberg correction was used for comparing alpha diversity, whereas permutational multivariate analysis of variance (PERMANOVA) was used to assess the effect of mitochondrial lineage on gut microbiome similarity. Individual foraging data (proportion of foraging of each eucalypt species) were used to visualize similarities in foraging patterns using nonmetric multidimensional scaling (NMDS) via vegan v2.6.4 in R v4.2.3 (R Core Team, 2023). PERMANOVA was performed using the vegan package to assess the effect of mitochondrial lineage on foraging patterns.

To investigate the relationship between the diversity of individual eucalypt foraging and alpha diversity of the gut microbiome, Spearman’s rank correlation coefficients between the alpha diversity (Shannon index) of the foraging and alpha diversity of the gut microbiome were calculated. In this analysis, the Shannon index was used to indicate the diversity of individual eucalypt foraging. Analysis of composition of microbiomes (ANCOM) (Mandal et al., 2015) was used to determine whether gut bacteria with characteristic relative abundances exist in each mitochondrial lineage or management group. All codes are available in the Supplemental Information.

Mitochondrial phylogeny

We performed phylogenetic analysis of the mitochondrial D-loop to estimate the mitochondrial maternal origins of the analyzed koalas (Fig. 1). A phylogenetic tree was constructed using the sequences of five koala haplotypes from Japanese zoos as well as the sequences of 48 individuals used for phylogenetic analysis by Neaves et al. (2016). Thus, we revealed that koalas from Japanese zoos examined in this study have three different origins (northern-2, N = 2; central, N = 4; and southern, N = 9).

Gut microbiome diversity among mitochondrial lineages

We investigated the mechanism by which the mitochondrial maternal origins of koalas, which were determined using mitochondrial phylogenetic analysis, and management groups influence the gut bacterial microbiome (Dataset S1, S2). Seven indices of alpha diversity were calculated (Tables S2, S3, Fig. 2). The index that considers evenness revealed that alpha diversity was the highest in the southern lineage and lowest in the central lineage (Simpson Index; South = 0.87 ± 0.07, Center = 0.75 ±  0.08, North2 = 0.84 ± 0.04; Simpson Index of Evenness; South = 0.10 ± 0.04, Center = 0.04 ± 0.01, North2 = 0.06 ± 0.01). In contrast, the index considering richness revealed that the southern lineage had the lowest alpha diversity (Chao1; South = 103.55 ± 28.36, Center = 129.84 ±  19.18, North2 = 125.31 ± 9.13; Faith’s phylogenetic diversity; South = 6.78 ± 0.98, Center = 7.81 ±  0.92, North2 = 7.50 ± 0.54; and observed features; South = 91.95 ± 21.58, Center = 115.57 ± 18.18, North2 = 113.25 ± 5.36). This trend was more clearly demonstrated between the management groups (Figs. 2A, 2B) than the mitochondrial lineages (the pairwise Kruskal–Wallis test with Benjamini–Hochberg correction; Shannon Index; South = 4.32 ± 0.22, North = 3.96 ± 0.50, P = 0.002022; Simpson Index; South = 0.90 ± 0.02, North = 0.81 ±  0.09, P = 2.60E-07; Simpson Index of Evenness; South = 0.14 ± 0.03, North = 0.06 ± 0.02, P = 4.12E-10; Chao1; South = 89.18 ± 11.66, North = 123.30 ± 25.96, P = 0.000003; Faith’s phylogenetic diversity; South = 6.28 ± 0.50, North = 7.50 ± 0.98, P = 0.000001; and observed features; South = 79.80 ± 6.94, North = 109.67 ± 20.61, P = 2.28E-07).

Mitochondrial lineage have a significant impact on gut microbiome similarity (beta diversity), both qualitatively (Unweighted UniFrac) and quantitatively (weighted UniFrac) Management groups have a significant impact on gut microbiome similarity, both qualitatively and quantitatively (pairwise PERMANOVA tests; Unweighted UniFrac, pseudo- F = 16.698158, q = 0.001; number of permutations = 999, Fig. 3A; weighted UniFrac, pseudo- F = 27.86, q = 0.0015; number of permutations = 999, Fig. 3B). In the same way, mitochondrial lineages also have a significant impact (pairwise PERMANOVA tests; Unweighted UniFrac, Central vs Northern 2: pseudo-F = 8.03, q = 0.001, Central vs Southern: pseudo- F = 8.13, q = 0.001, Northern 2 vs Southern: pseudo- F = 7.29, q = 0.001; number of permutations = 999, Fig. 3C; weighted UniFrac, Central vs Northern 2: pseudo- F = 5.78, q = 0.0015, Central vs Southern: pseudo- F = 11.22, q = 0.0015, Northern 2 vs Southern: pseudo- F = 5.09, q = 0.002; number of permutations = 999, Fig. 3D).

Based on differences in gut microbiome similarity between mitochondrial lineages, we investigated whether there were bacteria that differed in relative abundance with mitochondrial lineages or management groups using ANCOM. The ANCOM results revealed significant differences in the relative abundance of 12 bacterial genera between the two management groups (Table S4, Fig. 3E). Of them, eight genera showed higher relative abundance in southern koalas, whereas the four other genera showed higher relative abundance in northern koalas. Further, two bacterial genera showed significant differences among mitochondrial lineages (Table S5, Fig. 3F). These two genera clearly showed differences among mitochondrial lineages; they are unknown genera belonging to the families Tannerellaceae (significant in southern, W = 119) and Rikenellaceae (significant in northern 2, W = 111).

Relationship between eucalypt foraging and the gut microbiome

We recorded foraging data on zoo koalas to investigate the relationship between eucalypt foraging and the gut microbiome (Fig. 4A; raw data in Dataset S2). The proportion of each eucalypt species in foraging varied greatly from individual to individual. These foraging data were used to visualize the similarity of NMDS foraging patterns (Fig. 4B). PERMANOVA showed that mitochondrial lineage significantly influences foraging patterns (F = 5.88, R² = 0.68, P = 0.014, number of permutations = 719).

We investigated the relationship between the diversity of eucalypt foraging and the gut microbiome diversity via Spearman’s rank correlation coefficients (Table S6, Fig. 4C). The result shows a positive correlation between the diversity of the foraging and diversity of the gut microbiome (Spearman’s rank correlation; ρ = 0.89, P = 0.009, Shannon entropy). There was a correlation with the indicator that considers evenness but not with the indicator that considers richness (Table S6), i.e., significant in Shannon, Simpson, Simpson’s index of evenness, and Pielou’s evenness index, but not in Chao 1, Faith’s phylogenetic diversity and observed features (Spearman’s rank correlation; Shannon entropy: ρ = 0.886, P = 0.009, Pielou evenness: ρ = 0.770, P = 0.036, Simpson: ρ = 0.830, P = 0.021, Simpson e: ρ = 0.830, P = 0.021, Chao1: ρ = 0.030, P = 0.479, Faith pd: ρ = 0.540, P = 0.133, Observed features: ρ = 0.030, P = 0.479).

We investigated whether gut bacteria at the genus level vary in relative abundance depending on the foraging proportion of the eucalypt species. We revealed that the relative abundance of some gut bacteria increased or decreased depending on the foraging proportion of each eucalypt species (Table S7, Fig. 4D). The number of gut bacterial genera that showed a significant positive correlation was one with E. microcorys and two with E. punctata (Spearman’s rank correlation with Benjamini–Hochberg correction; q < 0.05). Conversely, the number of gut bacterial genera that showed a significant negative correlation was four with E. camaldulensis, two with E. microcorys, and one with E. robusta (Spearman’s rank correlation with Benjamini–Hochberg correction; q < 0.05).