Variation in the microbiome of the urogenital tract of female koalas (Phascolarctos cinereus) with and without ‘wet bottom’

Koalas (Phascolarctos cinereus) are iconic Australian marsupials currently threatened by several processes. Infectious reproductive tract disease, caused by Chlamydia pecorum, and koala retrovirus infection are considered key drivers of population decline. The clinical sign of ‘wet bottom’, a staining of the rump associated with urinary incontinence, is often caused by chlamydial urogenital tract infections. However, wet bottom has been recorded in koalas free of C. pecorum, suggesting other causative agents in those individuals. Current understanding of the bacterial community of the koala urogenital tract is limited. We used 16S rRNA diversity profiling to investigate the microbiome of the urogenital tract of ten female koalas. This was to produce baseline data on the female koala urogenital tract microbiome, and to undertake preliminary investigations of potential causative agents of wet bottom, other than C. pecorum. Five urogenital samples were processed from koalas presenting with wet bottom and five were clinically normal. We detected thirteen phyla across the ten samples, with Firmicutes occurring at the highest relative abundance (77.6%). The order Lactobacillales, within the Firmicutes, comprised 70.3% of the reads from all samples. After normalising reads using DESeq2 and testing for significant differences (P < 0.05), there were 25 operational taxonomic units (OTUs) more commonly found in one group over the other. The families Aerococcaceae and Tissierellaceae both had four significantly differentially abundant OTUs. These four Tissierellaceae OTUs were all significantly more abundant in koalas with wet bottom. Importance This study provides an essential foundation for future investigations of both the normal microflora of the koala urogenital tract, and better understanding of the causes of koala urogenital tract disease. Koalas in the states of Queensland and New South Wales are currently undergoing decline, and have been classified as vulnerable populations. Urogenital tract disease is a leading cause of hospital admissions in these states, yet previously little was known of the normal flora of this site. Wet bottom is a clinical sign of urogenital tract disease, which is often assumed to be caused by C. pecorum and treated accordingly. Our research highlights that other organisms may be causing wet bottom, and these potential aetiological agents need to be further investigated to fully address the problems this species faces.


23
Koalas (Phascolarctos cinereus) are iconic Australian marsupials currently threatened by 24 several processes. Infectious reproductive tract disease, caused by Chlamydia pecorum, and 25 koala retrovirus infection are considered key drivers of population decline. The clinical sign 26 of 'wet bottom', a staining of the rump associated with urinary incontinence, is often caused 27 by chlamydial urogenital tract infections. However, wet bottom has been recorded in koalas 28 free of C. pecorum, suggesting other causative agents in those individuals. Current 29 understanding of the bacterial community of the koala urogenital tract is limited. We used 30 16S rRNA diversity profiling to investigate the microbiome of the urogenital tract of ten 31 female koalas. This was to produce baseline data on the female koala urogenital tract 32 microbiome, and to undertake preliminary investigations of potential causative agents of wet 33 bottom, other than C. pecorum. Five urogenital samples were processed from koalas 34 presenting with wet bottom and five were clinically normal. We detected thirteen phyla 35 across the ten samples, with Firmicutes occurring at the highest relative abundance (77.6%). 36 The order Lactobacillales, within the Firmicutes, comprised 70.3% of the reads from all 37 samples. After normalising reads using DESeq2 and testing for significant differences (P < The koala (Phascolarctos cinereus) is an iconic marsupial species endemic to Australia. 56 Northern koala populations, in the states of Queensland and New South Wales, are currently 57 declining due to impacts from disease and increased urbanisation. A significant pathogen of 58 koalas, Chlamydia pecorum, has been a main focus of koala infectious disease investigations 59 since its discovery. C. pecorum infection of the conjunctiva or urogenital tract can lead to 60 blindness and infertility in koalas, respectively, greatly impacting population fecundity and 61 survivability (1, 2). C. pecorum is commonly associated with the clinical sign known as 'wet 62 bottom' or 'dirty tail' (3). This staining or scalding of the rump is associated with cystitis due 63 to C. pecorum infection in some populations (4), but recently samples from a large number of 64 koalas from Victorian populations with mild wet bottom were negative via qPCR for C. 65 pecorum (5). In particular, koalas in a population considered at the time to be free of C. 66 pecorum (6) had a similar prevalence and severity of wet bottom to populations where C. 67 pecorum occurred in more than 35% of koalas tested. Further analysis demonstrated that 68 whilst wet bottom could be significantly linked to the detection of C. pecorum infection in 69 male Victorian koalas, this relationship did not exist in females (7). It may be that an 70 unidentified organism is causing these mild clinical signs of disease in koalas. To date there 71 has not been extensive research to determine the normal flora of the koala urogenital tract, 72 making it difficult to use traditional microbiological techniques to determine species of 73 interest. Modern sequencing technology, specifically 16S rRNA biodiversity profiling, was 74 used to improve our understanding of the microbiome of the urogenital tract of koalas, and to 75 undertake preliminary comparisons of the microbiome of female koalas with and without 76 mild wet bottom. 77

Clinical status of koalas 79
Urogenital samples previously collected from ten koalas as a component of population health 80 monitoring were selected from an archive of samples available at our institute (7,8 In total, 13 phyla were detected in the ten samples (Table 1), with Firmicutes occurring at the 104 highest relative abundance (77.61%). Just over a third of the OTUs were classified as 105 Firmicutes (95/254), followed by Proteobacteria (59/254) and the Bacteroidetes (35/254). 106 When samples were split into the two groups, koalas without wet bottom had 89.3% of reads 107 classified as Firmicutes, followed by OTUs which could not be assigned using the 90% 108 similarity threshold (5.2%) and Actinobacteria (3.5%). Koalas with wet bottom had 68.2% 109 reads assigned to OTUs classified as Firmicutes. The next two most prevalent phyla were 110 Proteobacteria (12.5%) and Bacteroidetes (12.2%), however these phyla were over-111 represented in two samples, biasing the total relative values. Deferribacteres were detected in 112 only one sample (Koala 70, wet bottom present) and Acidobacteria were only detected in two 113 (one clinically normal koala and one displaying wet bottom). Armatimonadetes was detected 114 in three koalas without wet bottom, but in none of the five diseased koalas. These three phyla 115 were detected at the lowest relative abundance across the ten samples. Data for relative read 116 abundance for OTUs that could be taxonomically assigned to a genus level and occurred at a 117 percentage of 0.01% or more in either group can be found in Table 2. This shows that the 118 order Lactobacillales, and within that the genus Aerococcus, had the highest proportion of 119 relative reads. 120

Richness and diversity 121
Species richness within each sample is described in Table 1 Table 3 for individual alpha diversity 131 values). Results detailing abundance for all OTUs detected in koala urogenital samples are 132 recorded in supplemental material Table S1. 133 Fewer than half of the OTUs detected across the two sample groups were shared between 134 them (112/254) (Figure 2). At a read depth of 160,000 there was a significant difference 135 between the microbial communities in koalas with wet bottom compared to those without, 136 based on the results of a 10,000 permutation PERMANOVA using Bray-Curtis dissimilarity 137 (F = 4.92, P = 0.019) and unweighted (qualitative) UniFrac distances (F = 1.62, P = 0.031). 138 There was no significant difference detected when using weighted (quantitative) UniFrac 139 distances (F = 1.51, P = 0.061). 2D and 3D principal coordinate analysis (PCoA) graphs 140 comparing koalas with and without wet bottom are shown in Firmicutes as the most dominant phylum across all samples. This was followed by 146 Proteobacteria and Bacteroidetes ( Figure 4). Overall there were 25 OTUs with significant 147 (Benjamini and Hochberg (9) adjusted P < 0.05) over-representation or under-representation 148 in wet bottom affected koalas, in comparison to clinically normal koalas, based on these 149 normalised read counts (Table 4). Of those OTUs, when assessing absolute read count, six 150 occurred only in koalas with wet bottom, whilst eight occurred only in koalas without wet 151 bottom ( Table 4). All normalised read values can be found in supplemental material Table  152 S2, and all statistical comparisons in supplemental material Table S3. 153

154
Previous assessment of the koala microbiome has focused on the digestive system of koalas 155 comparing either two free ranging animals from northern populations (10) or two captive 156 koalas in Europe (11), from which the ocular microbiome was also assessed. This study is the 157 first investigation of the microbiome of the urogenital tract of the female koala using modern 158 high-throughput techniques, and only the second to assess the urogenital tract of a marsupial, 159 with the tammar wallaby (Macropus eugenii) investigated previously using terminal 160 restriction fragment length polymorphism analysis (12). The majority of reads in our sample 161 set were classified in the order Lactobacillales (72.1%). This dominance of Firmicutes 162 mirrors what has been seen in the human vaginal microbiome (13). In humans, the acidic pH 163 of the genital tract is maintained by these lactic acid producing bacteria, which in turn is 164 thought to play a role in preventing pathogenic infection (14). It appears from our sample set 165 that koalas have a different family within the Lactobacillales, possibly performing a similar 166 role. The most common family within our classified OTUs, in terms of either relative or 167 normalised read abundance, was Aerococcaceae, whilst in humans the Lactobacilli dominate 168 the reproductive tract. Within the Aerococcaceae, the genera Aerococcus and Facklamia 169 were both represented in the top four most abundant OTUs. For all four significantly 170 differentially abundant Aerococcus spp. OTUs, the same OTU could be detected in at least 171 4/5 (80%) of the converse sample group in absolute reads. For example, OTU 4, an 172 Aerococcus sp. occurred in all ten koala samples, but was present in significantly higher 173 quantities in clinically normal koalas after normalisation (P = 0.004). Whether specific 174 Aerococcus spp. that are over or under-represented are an important factor in terms of disease 175 presence requires further investigation. The production of hydrogen peroxide by commensal 176 Lactobacillus species is thought to play a role in reducing the successful establishment of 177 sexually transmitted diseases in humans (15, 16), and it has been shown that Aerococcus spp. 178 are involved in hydrogen peroxide production (17, 18). In humans Aerococcus spp. have also 179 been associated with disease, such as Aerococcus urinae, which can cause urinary tract 180 infections (19) and septicaemia (20). Investigations into the urinary microbiome of women 181 with and without 'urgency urinary incontinence' found that Aerococcus spp. were detected 182 more frequently in cases where disease was present (21). In our study, the four Aerococcus 183 spp. OTUs that had significantly different normalised abundance were evenly split, with two 184 having higher abundance in koalas with wet bottom and two in koalas without wet bottom. 185 The role of organisms within this family as opportunistic pathogens in koalas cannot be ruled 186 out. 187 The Aerococcus were the most common genus amongst those OTUs with significant 188 differential abundance after normalisation using DESeq2. The representative sequences of 189 these four OTUs did not match known species within the Aerococcus genus, using the 190 Greengenes database, with an identity greater than 90%, suggesting that these may represent 191 novel species. This is not unexpected, as the culture of organisms from the koala urogenital 192 tract has been limited to only a small number of studies, with the majority focused on 193 diagnosing what was later deemed to be chlamydial infection (22)(23)(24). Efforts in culturing 194 novel bacteria from koalas have focused primarily on its gut microbiome (25) The average number of OTUs detected in our samples is difficult to compare to other 210 publications investigating koala microbiomes. This is both due to the impact that sample site 211 differences would have on OTUs present, as well as the method of OTU classification used. 212 For instance, previous research on the koala intestinal microbiome used QIIME for analysis 213 of 454 pyrosequencing reads (10) and detected 1855 OTUs, after removal of chimeras and 214 singletons, from caecum, colon, and faecal samples. Similarly, an Illumina based study of 215 microbiomes from ocular, oral, rectal and faecal samples from two captive koalas found OTU 216 numbers ranging between 597 to 3,592, with a median of 1,456 (11). The average raw read 217 numbers per sample assessed in these projects ranged from 12,831 (454 pyrosequencing) to 218 323,030 (Illumina). Our own average raw reads per sample were within that range (229,561), 219 suggesting that the OTU differences between our studies are either associated with the 220 sample site (urogenital versus digestive tract) or clustering methodology used. We employed 221 UPARSE due to its demonstrated ability to correctly identify OTUs in a mock community 222 and minimise spurious OTUs (30). Whilst there did not appear to be any strong clustering on 223 our 2D or 3D PCoA plots, comparisons of the beta-diversity between groups highlighted that 224 the makeup of the communities was significantly different when assessing both Bray-Curtis 225 dissimilarity and unweighted UniFrac distances. These metrics assess presence/absence of 226 OTUs between groups, with UniFrac also considering phylogenetic distance between OTUs 227 present. Weighted UniFrac distances, which considers the abundance of individual OTUs, 228 were not significantly different between groups. Therefore, koalas with and without wet 229 bottom appear to have a significant difference in which OTUs are present in the samples, but 230 not necessarily the abundance of OTUs between samples. Two samples had widely different 231 OTU profiles (koala 49 and 70). This finding may support the hypothesis that wet bottom in 232 female koalas without C. pecorum may be caused by more than one aetiological agent (5, 31). 233 Further investigations to examine this hypothesis are indicated but require access to a large 234 number of appropriately collected and stored samples. Such sample sets are currently not 235 available for this species. 236 It could be argued that the skewed relative abundance of Proteobacteria and Bacteroidetes in 237 the samples from koala 49 and 70, respectively, could be a result of swab contamination with 238 faecal material, which would impact diversity inferences. The human microbiome project 239 identified that reads from stool samples were predominately from the Bacteroidetes phylum 240 (32), and the most recent assessment of the koala rectal microbiome found these two phyla to 241 be the most abundant in samples taken from both koalas assessed (11). In koalas, the 242 urogenital tract is accessed through the cloaca, which also contains the rectal opening. This 243 makes faecal contamination difficult to avoid during sample collection. Future studies of the 244 urogenital tract microbiome would benefit from either taking control samples from the 245 rectum of the koala being sampled, or inverting the cloaca so that the urogenital opening is 246 more easily accessible, as described previously for the tammar wallaby (12). In that study, 247 approximately a quarter of phylotypes (26/96) were detected in both the urogenital and rectal 248 samples, suggesting that bacteria being detected at multiple sites in marsupials is not unusual. Science Animal Ethics Committee, application ID:1011687.1, and all sample collection was 284 conducted following the Australian code for the care and use of animals for scientific 285 purposes, 8th edition (36). Wet bottom score was assessed using a scoring system as 286 previously described (37). These wet bottom scores grade the clinical findings relating to wet 287 bottom from 0 (absent) to 10 (most severe). For the purpose of this study, koalas were 288 grouped into wet bottom 'present' and wet bottom 'absent' categories. After screening all 289 samples for Chlamydiaceae using a previously described qPCR (8, 38), we selected ten 290 samples from female koalas where no Chlamydiaceae was detected. We used five samples 291 collected from koalas showing no clinical signs of urogenital disease and five samples 292 collected from koalas that showed clinical signs of wet bottom (Table 1). 293

Amplification and sequencing 294
DNA extraction and amplification from the swab samples was performed commercially by 295 The Australian Genome Research Facility (Australia). Variable regions three and four of 296 bacterial 16S rRNA were amplified using primers 341F (5' CCTAYGGGRBGCASCAG 3') 297 and 806R (5' GGACTACNNGGGTATCTAAT 3'). Sequencing was performed on the 298 Illumina MiSeq platform to produce paired end reads of 300 bp (2 × 300 bp). 299

Quality filtering and OTU assignment 300
Quality filtering and operational taxonomic unit (OTU) assignment was undertaken using a 301 mixture of scripts and algorithms available in the programs USEARCH 8.1 (39) and QIIME 302 1.9.1 (Quantitative Insights Into Microbial Ecology) (40). Unless otherwise stated, default 303 settings were used for all scripts. Read processing to reduce errors was undertaken as 304 described by Edgar and Flyvbjerg (41). The forward and reverse 300 bp paired-end reads for 305 each swab sample were merged using the USEARCH script fastq_mergepairs. In this 306 process, the Phred score of overlapping bases is recalculated to improve error calling. Bases 307 with the same nucleotide called in both the forward and reverse reads have an increased 308 recalculated score, and those with disagreements are reduced. This increases confidence in 309 the calculated error probability of the merged reads. Primers were then trimmed from the 5' 310 and 3' ends of the merged reads using seqtk (https://github.com/lh3/seqtk). Trimmed reads 311 were filtered for quality using the USEARCH script fastq_filter. This script filters reads 312 using the maximum expected errors per merged read. The number of expected errors is 313 obtained by the sum of the Phred derived error probability. If the expected number of errors 314 is less than one, then the most probable number of errors is zero (41). We utilised a maximum 315 expected error threshold of 1, resulting in reads with an error probability of 1 or greater being 316 removed. In addition to using the number of expected errors for filtering, trimmed reads 317 shorter than 400 bp were discarded. Unique reads within the entire sample set were assigned 318 OTUs using the USEARCH algorithms derep_fulllength and cluster_otus (30), with a 319 minimum identity of 97% for clustering, or a cluster radius of 3.0. Chimeras are filtered from 320 the sample set within the cluster_otus command using the UPARSE-REF maximum 321 parsimony algorithm (30). Singletons were excluded from OTU clustering due to the high 322 likelihood that they contain errors (30, 41). The merged/trimmed reads from each swab 323 sample, including the previously excluded singletons were matched with the clustered OTUs 324 using USEARCH script usearch_global, with a threshold of 97% identity to group a read 325 into a specific OTU. The taxonomy of each OTU was determined by using the QIIME script 326 assign_taxonomy.py in conjunction with the Greengenes taxonomy database (version 13_5, 327 97% clustered OTUs) (42). This script utilises the UCLUST algorithm (39) to identify a 328 consensus taxonomy of the reads within an OTU against the curated database, based on a 329 similarity of 90% and a minimum consensus fraction of 0.51. Chloroplast and mitochondrial 330 OTUs were removed from the dataset using the QIIME script 331 filter_taxa_from_otu_table.py. 332

Read normalisation and analysis 333
Read data was assessed using three different methods. Relative abundance was utilised to 334 compare basic phylum presence in each sample. Rarefaction of reads was undertaken, using 335 multiple_rarefactions.py QIIME script, to assess alpha and beta diversity at a set read level. 336 Negative-binomial normalisation of reads, using DESeq2 (43) as recommended by 337 McMurdie and Holmes (44), was performed using the QIIME script normalize_table.py. For 338 rarefactions, reads within each sample are subsampled (without replacement) every 5000 339 reads, from 5000 to 250,000 reads. This represented the maximum number of reads present in 340 the sample with the most reads (rounded down to the nearest value divisible by 5,000). At 341 each step, 100 permutations were undertaken. Alpha-diversity metrics (including species 342 richness, Chao1 (45), phylogenetic diversity and Shannon's diversity (46)) were generated 343 for each step. Comparisons of these values were undertaken using values obtained after 344 subsampling to a depth of 160,000. This equalled the sample with the fewest reads (rounded 345 down to the nearest value divisible by 5,000). Non-parametric comparisons of mean alpha 346 diversity metrics between the two sample groups (wet bottom present or absent) were 347 undertaken with the compare_alpha_diversity.py QIIME script. This script utilised a non-348 parametric two sample t-test with 10,000 Monte Carlo permutations to determine whether the 349 alpha diversity was significantly different between the two groups (wet bottom 350 present/absent) at a depth of 160,000 reads. Beta-diversity was assessed at the same depth as 351 above (160,000 reads) using the beta_diversity_through_plots.py QIIME script, in which 352 both unweighted and weighted UniFrac distances (47) were assessed. Bray-Curtis 353 dissimilarity (48) between samples was also assessed. The analysis of beta-diversity required 354 a phylogenetic tree. For this, an alignment of representative sequences of each OTU was 355 created with PyNAST (49) and UCLUST using the align_seqs.py QIIME script. A tree was 356 produced from this alignment using FastTree (50), and used as input for beta-diversity 357 analysis. beta_diversity_through_plots.py produced distance matrices for each of the tests 358 (UniFrac and Bray-Curtis), from which principal coordinates and eigen values could be 359 calculated. PCoA plots using the 2 or 3 most influential principal coordinates were drawn 360 from the resulting distance matrices either using either the make_2d_plots.py QIIME script, 361 or within the beta_diversity_through_plots.py script using EMPeror 9.51 software (51), 362 respectively. Distance and dissimilarity metrics were used to compare the microbial 363 communities between the two groups by utilising the permutational ANOVA 364 (PERMANOVA) method within the compare_categories.py QIIME script, with 10,000 365 permutations. Statistical comparisons of the differential abundance of OTUs between koalas 366 with and without wet bottom utilised DESeq2 within the QIIME script 367 differential_abundance.py. These comparisons aimed to determine OTUs which were over-368 represented in either group. Statistically significant results, using DESeq2s negative binomial 369 Wald test, were based on P-values < 0.05, and were adjusted for false discovery within the 370 script, using the method described by Benjamini and Hochberg (9). 371 The NCBI nucleotide database (52)  Highly efficient Aerococcus viridans L-alpha-glycerophosphate oxidase production in 436 the presence of H2O2-decomposing agent: purification and kinetic characterization.  Table 1. Koala wet bottom score, read metrics and relative abundance data from ten samples submitted for 16S rRNA amplicon sequencing. All koalas were female and sampled from French Island, Victoria, Australia in 2011.    # P value are from negative binomial Wald test, adjusted using the false discovery rate calculation described by Benjamini and Hochberg (9) * OTU was detected with significantly higher normalised read counts in koalas with (WB present) or without (WB absent) wet bottom ^ Organism with the lowest e-value detected using a MegaBLAST search of the NCBI nucleotide database, the nucleotide identity compared to the representative sequence, and the accession number of the hit    . DESeq2 normalised read counts of phyla detected in koala urogenital swab samples. Phyla with fewer than 2% relative reads within each sample have been excluded for clarity. Reads were characterised into taxanomic groups using QIIME (40), utilising Greengenes (42) as a reference database. Koalas 1 -5 were clinically normal (wet bottom absent), whilst koalas 31 -70 had wet bottom