Characterization of shifts of koala (Phascolarctos cinereus) intestinal microbial communities associated with antibiotic treatment

Sample collection

Samples were collected from koalas admitted to the Australia Zoo Wildlife Hospital in Beerwah, Queensland, Australia and the Port Macquarie Koala Hospital, Port Macquarie, New South Wales, Australia. Koalas were lodged in outdoor, single-occupancy enclosures. They were fed a variety of fresh eucalyptus leaves that were collected daily from local forests. Standard disease transmission prevention procedures were used at each hospital. Only individuals that were sampled for at least 21 days were included in this study because we thought that given koalas’ slow metabolisms and digestive mobility, we would still be able to capture immediate changes to the intestinal microbiome in this three-week period. However, this resulted in a limited number of control koalas because those that were not treated with antibiotics were released or died sooner than our 21-day cut off.

Seven koalas from the Australia Zoo Wildlife Hospital were included, six of which received the antibiotic chloramphenicol for the treatment of Chlamydia and one that did not receive any antibiotic treatment. Treatment involved a single, daily treatment of 60 mg/kg for varying durations; mL of the dosage is provided in Table 1. Four koalas were sampled from the Port Macquarie Koala Hospital, three of which were administered a combination of chloramphenicol and enrofloxacin antibiotics for the treatment of Chlamydia and one that did not receive any antibiotic treatment. Treatment involved a single, daily dose of 10 mg/kg of enrofloxacin for 2 weeks, followed by a single, daily dose of 60 mg/kg of chloramphenicol for varying durations (Table 1).

Two animals from Australia Zoo Wildlife Hospital (antibiotic-treated) and one Port Macquarie Koala Hospital animal (also antibiotic-treated) were euthanized during the course of this study. These animals exhibited too poor of a health status that they were determined to be unfit for survival by their respective veterinarians.

A total of 141 fecal samples were collected for analysis in this study from the eleven animals receiving care. These samples were composed of fresh fecal pellets (<15 min since excretion) collected every three days, beginning on the day each koala was admitted (i.e., immediately prior to first dosage of antibiotics) through to the day each koala was either released back to the wild or deceased. Fecal material was collected from the floor or branch with gloved hands and sealed in a small, sterile plastic bag. Samples were transported to a lab at the University of the Sunshine Coast on ice (<1 h) from the zoo and stored at −80 °C (Table 1).

Additionally, 17 built environment samples were collected from individuals’ enclosures. Samples were collected by a 10-second dry swabbing of floors and branches with a gloved hand, and sealed in a small, sterile plastic bag. Leaf samples were collected as well by removing a full leaves from the branches provided for each enclosure and sealing them in a small plastic bag. Samples were transported to a lab at the University of the Sunshine Coast on ice (<1 hr) from the zoo, and stored at −80 °C. These samples are included in Fig. 1, but were otherwise not analyzed in this report.

DNA extraction

Based on examination of previous work, we concluded that fecal pellets were the best possible proxy for the overall intestinal microbial community without the use of invasive collection techniques (Barker et al., 2013). To collect fecal material from the inside of the fecal pellet, the pellet was placed partway into a sterile 2 ml Eppendorf tube, the pellet was broken in half, and the innermost material from the middle of the pellet was removed with clean, sterile forceps. For each sample, approximately 150 mg of the fecal pellet was transferred with sterile forceps into a new sterile 2 ml Eppendorf tube. DNA was then extracted from this inner-pellet material using the Qiagen QIAamp Fast DNA Stool Mini Kit (following the manufacturer’s protocol). DNA was eluted in a final volume of 100 µl in 1.5 ml LoBind Eppendorf tubes and stored at −80 °C.

To confirm that there was bacterial DNA in the samples, we first tested for PCR amplification of 16S rRNA genes using full-length eubacterial 16S rRNA gene PCR primers, 27F (5′-AGAGTTTGATCMTGGCTCAG-3′) and 1492R (5′-TCNGGYTACCTTGTTACGAC TT-3′). Amplifications were carried out in an Eppendorf Mastercycler^® PCR Cycler in 50 μl reactions containing 5 μl of the eluted DNA from each samples’ extraction, 18 μl sterile Mill-Q^® H₂0, 25 μl MyTaq™ HS Mix 2x (Bioline, London, UK), and 1 μl of 0.2 μM concentration of each primer 27F and 1492R (Sigma-Aldrich, St. Louis, MO, USA). The cycling conditions were: (1) 95 °C for 10 min, (2) 35 cycles of 45 s at 95 °C, 1 min at 50 °C, and 2 min 20 s at 72 °C, (3) a final incubation at 72 °C for 7 min and (4) holding at 4 °C upon completion. A negative control was included in every reaction which replaced the 5 μl of DNA with 5 μl of sterile Mill-Q^® H₂0.

The PCR fragments were visualized by 45 min of 120 V electrophoresis on a 1% agarose gel stained with ethidium bromide. All samples showed a prominent band at the expected length (∼1,465 base pairs) for bacterial 16S rRNA genes; no bands were visible with the negative control of replacing the 5 μl of DNA with 5 μl of H₂0.

Sequencing

PCR amplification of the V4 region of the 16S rRNA gene was performed using primers 515 F and 806 R, as recommended by Caporaso et al. and modified by the addition of a custom barcode system described previously by Lang, Eisen & Zivkovic (2014) and Caporaso et al. (2012). In addition to primers, Invitrogen Platinum SuperMix was also used to perform PCR on 5 ng of DNA for each sample following previously established protocols (Lang, Eisen & Zivkovic, 2014). PCR cycling conditions were (1) 95 °C for 2 min, (2) 30 cycles of 95 °C for 30 s, 55 °C for 30 s, 72 °C for 1 min, and (3) 72 °C for 3 min. PCR cleanup and normalization was performed according the manufacturer protocol for the Invitrogen 96 well SequalPrep Normalization Plate. The resulting DNA elution of all pooled samples was further purified and concentrated according to the manufacturer protocol for the NucleoSpin^® Gel and PCR Clean-up Kit (Macherey-Nagel, Duren, Germany). All DNA was quantified with Qubit^® dsDNA HS Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA). Purified DNA of a final concentration of 34.1 ng/µl was submitted for sequencing at the UC Davis Genome DNA Technologies Core in a multiplexed Illumina MiSeq lane. We generated 8,889,513 250 bp pairs of raw reads and 7,020,375 253 bp merged reads, using a custom script (available in Github: https://github.com/gjospin/scripts/blob/master/Demul_trim_prep.pl) to assign each pair of reads from the custom dual barcode system to individual samples. All subsequent analysis of sequences was done on the merged reads.

Data analysis

Sequences were analyzed using the QIIME (Quantitative Insights Into Microbial Ecology) version 1.9.1 workflow (Caporaso et al., 2010). Quality filtering and chimera removal were performed on the sequencing reads before downstream analysis. Operational taxonomic units (OTUs) were picked at 97% similarity with QIIME’s script (pick_otus_through_otu_table.py) and clustered using the UCLUST algorithm (Edgar, 2010). Taxonomy was assigned using the taxonomy assignment script in QIIME (assign_taxonomy.py) with the open reference database, BLAST and reads that were not assigned a taxonomy via that approach were assigned using the RDP Classifier (Altschul et al., 1990; Wang et al., 2007). OTUs with five or fewer sequence representatives were filtered out. Additionally, the data was rarefied to 5,000 reads per sample. OTUs that did not have a taxonomic assignment at the ‘order’ level were filtered out.

Analysis was performed on samples from 11 individual koalas, nine of which had been treated with antibiotics and two that had not. Samples were collected every three to four days from the day the koala was first administered antibiotics until they were released or deceased; there were different numbers of samples collected from individual koalas because individual treatment time was varied. To normalize the analysis, only seven samples were used from each individual: the first sample, the last sample, and five samples that were selected to provide roughly equal chronological spacing from the treatment period. To ensure we only included OTUs capable of conveying group-level effects, we filtered the relative OTUs abundance table to only include OTUs that are present in two or more samples. All OTUs that did not satisfy this threshold were grouped together (sum of their abundances per sample) so that the underlying distribution of OTU abundances across samples did not change. This resulted in 1,511 OTUs across 77 samples. All negative controls for the DNA extraction kit, PCR, and sequencing were below the minimum 5,000 counts per sample that we set (see above) and thus were excluded from downstream analysis. Koalas R, P, and J were treated with a second antibiotic in addition to the same antibiotic as koalas A, B, D, E, F and G. However, calculations using the software R (R Core Team, 2013) determined no statistical difference existed between these two treatment groups (see ‘Results’) so they were pooled into a single antibiotic group.

Identifying predictive OTUs of fate: random forest analysis

To test whether a given set of OTUs was predictive of koala fate, and thus identify correlations between survival and intestinal microbiome composition, we applied the Random Forest model (available in Python’s scikit-learn package (Pedregosa et al., 2011)), a supervised machine-learning technique, to our dataset. In our model, we assigned OTUs as the features and identified the set of OTUs that can more accurately predict the fate (i.e., lived or died after treatment) of each individual. Using the Random Forest model, we examined the microbiome of all of the individuals in the antibiotic treatment group and determined the set of OTUs that create the highest distinction among the individuals that survived and the individuals that died.

The goal of the Random Forest classifier is to learn dependencies and complex relationships (both linear and nonlinear) among the features (here, relative OTU abundances) and find a set of features that are the most discriminatory among the groups and can improve the predictive accuracy of the model. It is important to note that the set of OTUs found in Random Forest model can be a combination of both abundant and relatively rare OTUs across samples.

Finally, there is an importance score assigned to each OTU that shows how predictive this OTU is in classifying the output. We measured the success of the Random Forest model by k-Fold cross-validation; this includes training the model on a subset of samples and then using the patterns learned to classify the remaining samples that were not used in the training step.

Diversity analysis

For the purpose of tracking antibiotic effects on the intestinal microbiome diversity, the mean alpha diversity of each sample was measured by the Shannon Diversity Index using the vegan package in R to measure variation within samples (R Core Team, 2013; Oksanen et al., 2013) (Fig. 2). Statistical significance of the differences in alpha diversity between groups was calculated using a Welch two-sample t-test in R. To determine variations in microbial communities between samples, principal-coordinate analysis (PCoA) plots were generated for weighted UniFrac, unweighted UniFrac, and Bray-Curtis distances using the QIIME software (Caporaso et al., 2011). Statistical significance of the differences between groups was calculated by either PERMANOVA or ANOSIM with 9,999 permutations in QIIME using the compare_categories.py script (Caporaso et al., 2010).

OTU co-occurrence network analysis

We investigated the potential interactions among various microbial taxa using network analysis. It has been shown that studying OTU co-occurrences patterns using network analysis provides insight into dynamics of complex microbial systems (Barberán et al., 2012; Faust et al., 2015; Faust et al., 2012). Here, we built a network of the OTUs based on the presence/absence patterns of OTUs across koala samples using the first and the last sample of individuals. Using our OTU table, each sample can be defined as a binary vector of OTUs showing presence/absence of OTUs and thus each OTU can be defined as a binary vector of that OTU presence/absence patterns across samples. To obtain co-occurrence patterns between every pair of OTUs and thus building the OTU co-occurrences matrix, we used the dot product of every two OTU vectors. The co-occurrence matrix was further used to build the network.

Each OTU represents a node in the network and there exists an edge between every two nodes if they have co-occurred together. The edge weights are the dot product of each OTU vector pair and are representative of how many times those two OTUs co-occurred together. We built a total of three OTU co-occurrence networks for the following sample subsets: (1) The initial time points for all individuals that were subjected to antibiotic treatment regardless of fate; (2) the final time point for individuals who lived; and (3) the final time point for individuals who died. Note that since we have smaller number of individuals in the deceased group, we expect to observe more nodes in the network built from this group because it is more likely to see an OTU present in three samples than it is to observe it in six individuals (case of released group) or nine individuals (network built based on the initial time point of all individuals). Once the networks were built, we obtained clusters of OTUs that co-occurred together in each of these networks separately. We looked at the differences among clusters of OTUs co-occurring in the initial sample and compared that to OTU clusters that co-occurred in the last sampling. We then identified the intersection of the co-occurring OTUs among these three networks.

Diversity analysis

Diversity analysis was performed on samples all time points for every koala as outlined in Table S1, with outliers removed where applicable, unless otherwise stated.

Random forest analysis

Using the Random Forest model with a 3-fold cross-validation, we determined that individuals in the antibiotic treatment group were classified accurately into two groups, those that survived and those that died, with an accuracy of 92 percent and a total of 516 predictive OTUs. The OTU that was most predictive of fate was classified (using the QIIME taxonomy assignment script (assign_taxonomy.py) (Caporaso et al., 2010)) as Lonepinella (sp. koalarum). The distribution of the feature importance was scored with the gini impurity criterion (measures the randomness of false label assignment). Using this distribution, we identified the most predictive top 20 percent OTUs and their distribution across samples (Fig. 3).

OTU co-occurrence network analysis

The resulting bacterial networks for co-occurring OTUs all had a density of one, meaning that there existed an edge between every two nodes (see ‘Methods’). Because all three networks (antibiotic treated koalas (N = 9), koalas that lived (N = 6), and koalas that died (N = 3)) had a density of one, they were defined as cliques and therefore considered to be one cluster. In more detail, the network that was built based on the initial time point from nine individuals resulted in 34 nodes with 561 edges. The average weighted degree was 297. Similarly, for the network built for the final time point for three individuals that were deceased, there were 100 nodes with 4,950 edges with an average weighted degree of 297. For the network of six released individuals, there were 55 nodes with 1,485 edges and average weighted degree of 324. As a final step, we looked at the intersection of co-occurring OTUs among the networks we built. Twenty-four of the OTUs that were co-occurring at the initial time point still existed in all of the samples in both deceased and released individuals at the final time point. There were four OTUs we determined to be of most interest, New.ReferenceOTU131, New.ReferenceOTU23, New.ReferenceOTU271 and New.ReferenceOTU265, which existed in all samples at initial time point and existed in all released samples. However, they were absent from several samples of koalas that died. We note, the taxonomic identification suggests that New.ReferenceOTU265 is most closely related to sequences that were annotated as being in the Lonepinella genus. For more information on the OTUs that resulted from network analysis, see Table S1.

Diversity analysis

Random forest and co-occurrence network analysis

To identify the most predictive OTUs correlated with whether or not a koala survived antibiotic treatment, we performed a Random Forest analysis of OTU relative abundance and presence/absence. The Random Forest analysis revealed that an OTU identified as Lonepinella koalarum is the most predictive OTU of whether or not a koala lived or died following the administration of antibiotics. Koalas that died after antibiotic treatment had much lower relative abundance (sometimes even zero), of this L. koalarum OTU compared to koalas that survived antibiotic treatment (Fig. 3). We believe that this correlation of the abundance of L. koalarum in koala intenstines with surviving antibiotic treatment is potentially a key finding, as L. koalarum is known to be a tannin-degrading microbe in the koala intestinal (Osawa et al., 1995). Our finding that the set of OTUs predictive of koala survival classified with an accuracy of 92% suggests there is a strong correlation between intestinal microbiome composition and koala prognosis, which is consistent with our initial hypothesis.

In addition to the Random Forest analysis, an OTU co-occurrence network analysis also revealed that L. koalarum was correlated with whether or not koalas lived after antibiotic treatment. L. koalarum was one of only four OTUs that was present at the beginning and at the end of antibiotic treatment in koalas that lived, but absent in at least one koala that died after antibiotic treatment.

Another finding from our Random Forest analysis was that the OTU “Cyanobacteria_YS2” was in the top 30 most predictive of whether or not a koala lived following antibiotic treatment. Recently, it was shown that the group to which this OTU was assigned to is actually Melainabacteria, a new phylum of non-photosynthetic Cyanobacteria that has been found in numerous mammalian intestines (Soo et al., 2014; Di Rienzi et al., 2013).

Overall, the results of Random Forest and co-occurrence network analyses support our hypothesis that the administration of antibiotics, regardless of the combination of the treatment, was associated with a change in the presence and relative abundance of L. koalarum, a known tannin-degrader. This is not only important for veterinarians to consider when administering antibiotics to koalas, but also for the development of Chlamydia infection treatments that do not impact this critical intestinal microbe.

Potential confounding variables that may have influenced our results are sample storage (Lauber et al., 2010), DNA extraction method (Salter et al., 2014), PCR (Brooks et al., 2015; Brown et al., 2015), and sequencing (Schirmer et al., 2015). Due to not having full medical histories of the koalas in this study, it is also possible that unknown confounding variables (e.g., immune status or prior environment exposures), rather than antibiotics, contributed to our results. Other environmental variables that may be confounding factors that we were unable to test for include the effects of being handled by hospital staff, koala compliance with antibiotic treatment, differences between wildlife hospitals’ environments and procedures, and day-to-day care.

Conclusions

Our findings are consistent with numerous other papers that suggest antibiotic treatments can cause a disturbance to intestinal microbial communities e.g., Jakobsson et al., 2010; Dethlefsen et al., 2008; Jernberg et al., 2007; Sullivan, Edlund & Nord, 2001. Such disturbances are likely to be particularly important in species like koalas, where it is thought that the intestinal microbial community may be required for survival (in this case, detoxifying their food). This is in contrast to many other animals where the intestinal microbiome is important but not necessarily essential for survival. In addition to showing differences in antibiotic treated and untreated koalas, we also found differences in richness, evenness, and structure of intestinal microbial communities in antibiotic treated koalas for those that survived versus those that did not. In particular, we found that the relative abundance of some key OTUs are correlated with survival. This is consistent with our hypothesis that koala survival after antibiotic treatment is related to whether or not key OTUs persist after the antibiotic treatment. We believe this suggests that there may be a need to develop alternative treatments for koala Chlamydia infections without the use of antibiotics and/or supplementing the koala diet with probiotics during antibiotic treatment. Furthermore, our conclusions may be transferable to other species that consume high-tannin diets. Future studies about this topic could be more comprehensive in scale (e.g., sample size) and depth (e.g., metagenome rather than 16S rRNA gene sequencing). It would be valuable to track koalas over a longer period of time rather than just their time in captivity. Furthermore, an important area to investigate is the role of pap consumption in the colonization and microbial community structure changes of koala intestines. Given our findings, it would be particularly important to investigate the pap of antibiotic treated mother koalas and the impacts this has on joey health.