Characterising the Tasmanian devil (Sarcophilus harrisii) pouch microbiome in lactating and non-lactating females

Wildlife harbour a diverse range of microorganisms that affect their health and development. Marsupials are born immunologically naïve and physiologically underdeveloped, with primary development occurring inside a pouch. Secretion of immunological compounds and antimicrobial peptides in the epithelial lining of the female’s pouch, pouch young skin, and through the milk, are thought to boost the neonate’s immune system and potentially alter the pouch skin microbiome. Here, using 16S rRNA amplicon sequencing, we characterised the Tasmanian devil pouch skin microbiome from 25 lactating and 30 non-lactating wild females to describe and compare across these reproductive stages. We found that the lactating pouch skin microbiome had significantly lower amplicon sequence variant richness and diversity than non-lactating pouches, however there was no overall dissimilarity in community structure between lactating and non-lactating pouches. The top five phyla were found to be consistent between both reproductive stages, with over 85% of the microbiome being comprised of Firmicutes, Proteobacteria, Fusobacteriota, Actinobacteriota, and Bacteroidota. The most abundant taxa remained consistent across all taxonomic ranks between lactating and non-lactating pouch types. This suggests that any potential immunological compounds or antimicrobial peptide secretions did not significantly influence the main community members. Of the more than 16,000 total identified amplicon sequence variants, 25 were recognised as differentially abundant between lactating and non-lactating pouches. It is proposed that the secretion of antimicrobial peptides in the pouch act to modulate these microbial communities. This study identifies candidate bacterial clades on which to test the activity of Tasmanian devil antimicrobial peptides and their role in pouch young protection, which in turn may lead to future therapeutic development for human diseases.

external microorganisms, or exposure to anthropogenic disturbances increase variation between individuals [44][45][46][47][48] .These differences, however, can be accounted for by including location as a covariate in analyses where possible.
Here we build on previous investigations by utilising a larger sample size collected from wild, free-ranging devils.We aimed to characterise and compare the taxa and community diversity between the pouch microbiome of lactating and non-lactating devils.
Swabs were aseptically sealed in a sterile container and kept on ice in the field.They were then stored at -20 °C prior to extraction.A total of 55 samples were used in this study; 25 from lactating females and 30 from nonlactating females.A power analysis was conducted in the pwr package 52 in R (version 4.2.1 53 ) to determine the minimum sample size to detect large differences (equivalent to an effect size of 0.8) in community composition (detected using UniFrac distance) between two equal sized sample groups; 25 samples per group gave a power of 0.8 to detect a change at a significance level of alpha = 0.05.

DNA extraction and sequencing
Microbial DNA was extracted using the DNeasy PowerSoil Pro Kit (Qiagen, Hilden, Germany).Extractions were conducted according to the Qiagen protocol with the following modifications: (i) the rayon end of the swab was cut off into the PowerBead tube, and (ii) the microcentrifuge tube containing the lysis buffer (Qiagen solution CD1) and swab tip were homogenised for 10 min in an Omni Bead Ruptor 4 at 3 m/s to lyse the bacterial cells present on the swab.Negative control extractions in each extraction batch omitted the swab.Extraction products were tested for DNA quantity and purity using a NanoDrop 2000 Spectrophotometer (ThermoFisher Scientific, www.nature.com/scientificreports/Waltham, MA, USA), and DNA concentration was quantified with a Qubit 2.0 Fluorometer using the Qubit dsDNA HS Assay Kit (ThermoFisher Scientific).An in-house PCR was performed on a small subset of samples to ensure that amplifiable microbial DNA was present (protocol outlined in Additional file 2).Extracted DNA was sent to the Ramaciotti Centre for Genomics (University of New South Wales, Kensington, Australia) for library preparation and sequencing.PCR products were indexed using the KAPA HiFi Hot Start Readymix (Roche, Basel, Switzerland), and normalised and pooled for sequencing using the SequalPrep™ Normalization Plate Kit (Invitrogen, Waltham, MA, USA) following the manufacturer's instructions.16 s rRNA V3-V4 amplicon sequencing was performed using the Illumina MiSeq v3 platform (2 × 300 bp) to amplify the V3-V4 hypervariable region in bacteria only.This region was selected for sequencing due to the higher rate of diversity capture as compared to primer pairs targeting other hypervariable regions 54 .Negative extraction controls were also sequenced.The standard V3-V4 primer pair, 341f (5ʹ-CCT ACG GGNGGC WGC AG-3ʹ) and 805r (5ʹ-GAC TAC HVGGG TAT CTA ATC C-3ʹ) 55 , was used to generate sequencing libraries, with sequencing generating 300 bp paired-end reads.

Data processing and statistical analysis
QIIME2 (version 2022.2 56 ) was used to process and analyse the resulting raw reads.Demultiplexed paired-end sequence reads were merged, denoised, and quality filtered into ASVs using the DADA2 plugin 57 .Reads with a median quality score < Q30 were trimmed and truncated, producing trimmed reads that were 523 bp long.A feature table and phylogenetic tree were produced using the QIIME2 'feature-table summarize' and 'align-totree-mafft-fasttree' commands, respectively.ASVs were assigned to taxonomy using a V3-V4 amplicon-specific naïve bayes classifier trained on the SILVA 138.1 small subunit (SSU) reference database 58 using the QIIME2 plugin, RESCRIPt (Reference Sequence Annotation and Curation Pipeline 59 ).An amplicon specific classifier was generated due to its improved accuracy in taxonomic classification of sequences 60 .The metadata table along with the resulting feature table, taxonomy table, and unrooted tree were imported and compiled into a phyloseq object using the qiime2R (version 0.99.6 61 ) and phyloseq (version 1.40.0 62 ) packages in R. Features present in the extraction controls were removed using the prevalence method in the decontam package (version 1.16.0 63 ) at a threshold of 0.5 (Additional file 1: Figs.S3-S5).This threshold identifies ASVs as contaminants when their prevalence is higher in control samples than biological samples 63 .Additionally, mitochondrial, chloroplasts, other eukaryotic sequences, and ASVs classified under any kingdom other than bacteria were removed, as they are either not targeted by 16S sequencing and are therefore sequencing errors, or are likely to be artificially amplified while not actually contributing to the normal community composition 64 .Due to the unequal number of reads per sample, the feature table was normalised through rarefaction to a depth of 18,198 for diversity analysis.
Prior to conducting the focal analysis, exploratory analysis was undertaken to investigate any potentially confounding factors present in the dataset.A Welch's t-test was used to assess any significant differences in the age of the female devils between reproductive groups.An exact date of birth cannot be measured for wild devils, so female age in years was used, knowing that sexual maturity is reached at approximately two years of age 49 .A similar data exploration could not be conducted for location due to the categorical nature of the data.Instead, location was included as a covariate in all steps of analysis.
The relative abundance of ASVs were assessed per sample at different taxonomic ranks to visualise differences in dominant taxa and their relative abundance between sample groups.Additionally, taxa were agglomerated separately at each taxonomic level using the 'tax_glom' phyloseq function to extract relative abundance data per sample group.Alpha diversity (i.e., the diversity contained within a single sample) was assessed to compare the diversity of individual samples within each sample group.Observed ASVs 65 , Shannon Diversity Index 66 , and Faith's Phylogenetic Diversity 67 were selected as metrics to measure alpha diversity.The statistical significance of these metrics was compared between reproductive statuses and locations using generalised linear models (GLMs), with the response variable, alpha diversity, tested against the predictor variables of reproductive status and location.
Beta diversity (i.e., the similarity or dissimilarity between two communities) was assessed using weighted and unweighted UniFrac distances and visualised by Non-Metric Multidimensional Scaling (NMDS).UniFrac distance was selected over other beta diversity metrics as it measures the dissimilarity between microbial communities based on the unique and shared evolutionary history of their constituent species 68 .Both UniFrac distance matrices were used to determine if compositional differences were present, and if so, whether they were influenced by the abundance of taxa.Statistical significance of UniFrac dissimilarity between sample groups was tested using a Permutational multivariate analysis of variance (PERMANOVA) test with 999 permutations through the vegan package (version 2.6-2 69 ) in R. The homogeneity of group dispersions was tested using the 'betadisper' function from the vegan package as it may influence PERMANOVA results (Additional file 1: Fig. S1).Additionally, pairwise multilevel comparisons were conducted using the pairwiseAdonis package (version 0.4 70 ) to determine the spatial drivers of community dissimilarity.Due to the unequal sample sizes at some locations, we also conducted a pairwise comparison between reproductive status and beta diversity within the sites with larger sample sizes (i.e.Fentonbury, Kempton, and Narawntapu National Park).This approach was chosen instead of an overall comparison between all locations due to non-homogenous dispersion of unweighted UniFrac between locations.The non-homogenous group dispersions may confound any outcomes of an overall comparison.
Differentially abundant taxa were identified and statistically compared using DESeq2 (version 1.36.0 71 ) in R.This differential abundance analysis was conducted on non-rarefied sequencing data due to the alternate normalisation methods utilised by DESeq2 72 .Testing was conducted to determine differentially abundant taxa between lactating and non-lactating pouch microbial communities at multiple taxonomic ranks; phylum, class, order, family, genus, and ASVs.ASVs were agglomerated into taxonomic groups prior to analysis using phyloseq's 'tax_glom' function, hence ASVs were chosen instead of species to avoid agglomeration of unidentified species.
A total of 62 samples (55 devils + 7 negative extraction controls) were sequenced, resulting in 3,344,445 total reads.After filtering, trimming, decontamination, and removal of non-targeted taxa, 16,106 ASVs were identified throughout an average of 59,940 reads and a median of 56,189 reads per sample.The sequencing of the negative controls yielded a total of 1397 reads (an average of 200 reads per control) following which 12 ASVs were identified as contaminants and removed from the dataset before further analysis (Additional file 1: Table S4).Some of the contaminating genera identified by decontam are common contaminants of lab reagents and disposables, such as Burkholderia, Pseudomonas, and Streptococcus [73][74][75] .

Community diversity
While no distinct clustering is visible, the unweighted UniFrac NMDS plot indicates non-lactating samples to be more spread out than lactating samples (Fig. 4A).A PERMANOVA was used to determine if this distinction between the two reproductive groups was statistically significant.Unweighted UniFrac distance was significantly different between reproductive statuses (F = 1.84, p < 0.05) and locations (F = 1.20, p < 0.005; Additional file 1: Table S6A), however, the significant difference in dispersion (reproductive status: F = 5.86, p < 0.05; location: F = 8.31, p < 0.005; Additional file 1: Table S7) means that it cannot be concluded whether the significant difference in unweighted UniFrac distance is due to average community structure or variation in dispersion between groups (Additional file 1: Table S7A; Table S7B).Weighted UniFrac distance, however, was significantly different between locations (F = 1.75, p < 0.005) but not reproductive status (F = 0.87, p > 0.05; Fig. 4B; Additional file 1: Table S6B).The dispersion effects for weighted UniFrac distance were not significantly different (reproductive status: F = 3.67, p = 0.065; location: F = 1.15, p = 0.346; Additional file 1: Table S7).Due to the high NMDS stress score and unclear unweighted UniFrac results, a PCoA was run yet also failed to identify any clear clustering between treatment groups (Additional file 1: Fig. S6).The pairwise comparisons of UniFrac distance between reproductive status within each location separately found there to be a significant difference in the community composition of the pouch microbiome between lactating and non-lactating females at Narawntapu National Park (unweighted UniFrac: F = 1.40, p < 0.01; weighted UniFrac: F = 2.547, p < 0.05; Additional file 1: Table S8, Fig. S7), but the microbial communities in the pouches of lactating and non-lactating females from Fentonbury and Kempton were not dissimilar from each other (Additional file 1: Table S8).
The results from the DESeq2 differential abundance analysis found four orders, three families, one genus, and 25 ASVs to be differentially abundant between lactating and non-lactating pouches (Fig. 5  Table S9).Of the 25 ASVs identified as differentially abundant, 16 were from the phylum Firmicutes, six Proteobacteria, two Fusobacteriota, and one Actinobacteriota (Fig. 5).Clostridium and Tisierella ASVs were most commonly identified as differentially abundant.Clostridium botulinum and Clostridium paraputrificum were more abundant in lactating pouches, whereas Clostridium colicanis was more abundant in non-lactating pouches.Four unidentified Tisierella ASVs were more abundant in non-lactating pouches than lactating.Apart from Clostridium spp., ASVs from the same genus exhibited consistent positive or negative log2 fold changes.All ASVs identified as differentially abundant produced a BH adjusted p-value of < 0.0001 (Additional file 1: Table S9).

Discussion
Here, we aimed to characterise the microbiome of pouches from lactating and non-lactating female Tasmanian devils and to identify differences in ASV richness and diversity between the two groups.We found lactating and non-lactating pouches to consist of the same dominating taxa, yet lactating pouches contained fewer overall ASVs and were less diverse.We demonstrate the dominating phyla to be Actinobacteriota, Bacteroidota, Firmicutes, Fusobacteria, and Proteobacteria, comprising a combined total of 99% and 96% of the microbiome in lactating pouches and non-lactating pouches, respectively.Previous studies found the same five phyla to comprise ~ 90% of the microbiome in reproductively inactive devils 26 , while the only previous study to investigate the pouch microbiome of lactating devils did not report relative abundances at the phylum-level 28 .Our study had different average relative abundance in non-lactating pouches to previous studies for Firmicutes, 61.3% here compared to 28.4% 28 and 36.2% 26 , and Proteobacteria, 17.2% here compared to ~ 35% previously 26,28 .The relative abundance of Fusobacteria was within 1% of the relative abundance in non-lactating pouches as identified by Cheng et al. yet was ~ 18% lower than the relative abundance identified by Peel et al.The mean relative abundances of Bacteroidota and Actinobacteria in non-lactating pouches were within 3% of values from previous findings 26,28 .The differences observed in abundance of Firmicutes and Proteobacteria may be due to either differing analysis techniques between our study and previous work, or environmental conditions experienced by sampled individuals.Peel et al. 28 and Cheng et al. 26 clustered their sequences into OTUs, whilst we utilised ASVs, a more reproducible and contemporary method of denoising Illumina amplicon data 65 .OTUs and ASVs have been found to be comparable at higher taxonomic ranks, however, and exhibit highly similar alpha and beta diversity  S5.Visualisation was performed using ggplot2 (v 3.4.2).
Vol:.( 1234567890) metrics, producing similar ecological conclusions [76][77][78] .As Peel et al. 28 utilised V3-V4 amplicons sequenced on an Illumina MiSeq, as per this study, differences in analysis are unlikely to be driving the differences observed with the relative abundance of Firmicutes and Proteobacteria.The most likely explanation of the differences observed between our study and previous ones, is environmental variation between sampled individuals.Differences in results may also be influenced by the stage of lactation at which samples were collected.It is well documented that the immunological properties of marsupial milk can change throughout the stages of pouch young development (see 79,80 and references therein), and as such, the expression of antimicrobial peptides in the pouch may also fluctuate throughout lactation, yet the stage at which samples were collected for Peel et al. 28 is unknown.Here we sampled a greater number of individuals across multiple geographic locations, likely to be more indicative of the true patterns of microbial communities within the devil pouch.Cheng et al. 26 found microbial communities present in the pouch of captive devils to be significantly different to wild devils, with significant variation between some wild sites (N = 18).Just as Cheng et al. 26 identified differences in microbial communities of devil pouches across geographic locations, we identified samples from Narawntapu National Park to exhibit higher ASV richness and phylogenetic diversity.Narawntapu National Park was also the only location where the overall Colour represents reproductive status while shape indicates location.Results of statistical analysis can be found in Additional file 1: Table S6.Visualisation was performed using ggplot2 (v 3.4.2).
microbial communities associated with lactating and non-lactating pouches were significantly dissimilar from each other.To truly disentangle spatial differences in pouch microbiome, future studies should sample between geographic locations at the same time point to remove any seasonal breeding effects.Furthermore, it is uncertain whether the lower number of observed ASVs and diversity found in the lactating pouch were driven by lactation or by the age of the female devil.Since sampling took place during the breeding season at in-situ sites with small populations of devils, most of the adult females were reproductively active, making it challenging to collect samples from females of the same age to eliminate this variation.This is a limitation of this study, and future studies may be able to address this by collecting samples for lactating females during the breeding season and for non-lactating females outside of the breeding season.The previous study investigating changes in the pouch microbiome associated with lactation in wild devils found the level of diversity in pouches from both reproductive stages to be 'similar' , yet differences in abundance caused lactating pouches to exhibit a high degree of dissimilarity from non-lactating pouches 28 .This differs from our results which found a significant difference in species diversity between the microbiome of lactating and non-lactating pouches but failed to identify distinct community differences.Similarly, recent studies comparing the pouch microbiome of lactating and non-lactating pouches in koalas and southern hairy nosed wombats found significant differences in species diversity and richness, however, unlike these studies, we failed to identify distinct community differences between lactating and non-lactating pouches 29,31 .The lower species richness and diversity observed here between lactating pouches and non-lactating pouches indicates a reduction in the number and type of species present.The lower phylogenetic diversity means that the bacterial species present in the presence of pouch young are more closely related.This reduction in richness and diversity may be due to the proposed antimicrobial action of pouch secretions 28,[81][82][83] .Peel et al. 28 identified one antimicrobial peptide (cathelicidin, Saha-CATH2) as being expressed more highly in the pouch versus other tissues.This cathelicidin, however, showed minimal antimicrobial activity against the bacterial strains tested, some of which were found in the devil pouch, such as Staphylococcus aureus and Pasteurella multocida.Instead, it is possible that Saha-CATH2 may be active against bacteria that was not tested but is found at higher relative abundances in the non-lactating devil pouch (i.e.Clostridium spp.).In addition, Saha-CATH2 was found to be evolutionarily related to tammar wallaby and opossum cathelicidins 28 .A closely related tammar wallaby cathelicidin, MaeuCath8, was highly expressed in the female's pouch, suggesting that secretions from epithelial cells within the pouch may provide antimicrobial protection to pouch young 84 .The decreased phylogenetic diversity of microbial pouch communities during lactation may be explained by the expression of antimicrobial peptides exhibiting specific antibacterial action against specific bacterial taxa within the pouch.Future studies testing the antimicrobial properties of devil pouch secretions should be guided by the bacterial genera that were found to be significantly lower in lactating devil pouches; Vagococcus, Peptostreptococcus, W5053, Ignatzschineria, Moraxella, Wohlfahrtiimonas, Tissierella, Cetobacterium, and Clostridium sensu strico 1. www.nature.com/scientificreports/Many of the dominating phyla identified in the pouch microbiomes are highly prevalent in the skin microbiome of other vertebrates, with Actinobacteria, Firmicutes, Proteobacteria and Bacteroides dominating 99% of the human skin microbiome 85 .When comparing the devil pouch microbiome to the skin microbiome of other terrestrial mammals, it is evident that Firmicutes comprised a much larger proportion of the devil microbiome (60-70%) as compared to other terrestrial mammalian species such as the wild bank vole (Myodes glareolus; 28% Firmicutes 47 ) or microbat species, which typically saw Actinobacteria as the most abundant phyla 86 .The relative abundance of Firmicutes is high in the Tasmanian devil pouch, even compared to other marsupial species such as the southern hairy nosed wombat (> 10-35% in reproductively active and inactive pouches, respectively) and the koala (10-12% during anoestrus and early lactation).
Multiple genera identified within the pouch, irrespective of reproductive status, are implicated with human and wildlife diseases, posing a potential threat to the immunologically naïve pouch young.Vagococcus contains two species that are pathogenic to domestic mammals, birds, and fish 87,88 , while Peptostreptococcus spp.have been identified in human skin diseases 89 .Ignatzschineria and Wohlfahrtiimonas spp.are pathogenic species that are interestingly both associated with the microbiome of maggot-infested wounds in humans 90,91 .Other differentially abundant taxa have been identified in the mammalian gastrointestinal tract, such as Tissierella and Clostridium colicanis [92][93][94] .Tissierella has been associated with some infections in humans 92 , while the pathogenicity of C. colicanis has not been well-explored.Additionally, increased abundance of W5053 spp. in the mammalian urogenital tract has been associated with adverse reproductive conditions 95,96 .The decreased abundance of this genus in the pouches of lactating devils may illustrate a healthy reproductive environment.While some protection appears to be provided to the pouch young through the hypothesised mechanisms of antimicrobial secretions, some species of bacteria found in the lactating pouch are still potentially pathogenic.Species such as Clostridium botulinum and Corynebacterium urealyticum were more abundant in lactating pouches but are associated with paralysis and urological diseases in other vertebrates, respectively 97,98 .The impact of these species on the health of the immunologically naïve pouch young is unknown, yet with a 75% survival rate of pouch young to weaning 99 , it is likely that these bacteria are either non-pathogenic to this species or that other mechanisms of protection are employed by the female to protect the pouch young, potentially via antimicrobial peptides and/or the transfer of immunological compounds in the milk 79 .

Conclusions
This study is one of the largest investigations into the microbial changes associated with lactation in the marsupial pouch.Our findings show that the bacterial communities found in lactating devil pouches are comprised of taxa that are more closely related and contain less rare taxa when compared to pouches of non-lactating devils.We propose that these results may in part be due to the secretion of antimicrobial peptides that are secreted by the pouch skin.Further investigation into the antimicrobial action expressed in the devil pouch against the suite of bacteria characterised here will determine the bacterial taxa targeted by these protective mechanisms, furthering our understanding of the immunological development of marsupial pouch young.

Figure 1 .
Figure 1.Map of sites where pouch swab samples were collected from Tasmanian devils.Samples were collected from the following sites for microbial analysis; Narawntapu National Park (NP; N = 11), Stony Head (SH; N = 5), Fentonbury (FB; N = 18), Kempton (KP; N = 16), and Buckland (BL; N = 4).A breakdown of the number of samples collected for each sampling group can be found in Additional file 1: Table S2.Colour corresponds to month when fieldtrip took place, scale bar represents kilometres, and arrow indicates north.Graphic was created using QGIS (v 3.22) and Australian Bureau of Statistics Digital Boundary Files 100 .

Figure 2 .
Figure 2. Taxonomy bar plots displaying the relative abundance of bacteria present in Tasmanian devil pouches across reproductive stages at three taxonomic levels; (A) Phylum, (B) Family, and (C) Genus.Each column represents one sample.Samples are organised according to reproductive stages; lactating and nonlactating.Location is indicated below the x-axis as follows; BL = Buckland, FB = Fentonbury, KP = Kempton, NP = Narawntapu National Park, and SH = Stony Head.Taxa with a median relative abundance of < 0.25% were merged into 'Other' .Visualisation was performed using ggplot2 (v 3.4.2).

Figure 3 .
Figure 3. Alpha diversity metrics of pouch microbial compositions per reproductive status.Whisker box plots depict (A) Observed Amplicon Sequence Variant Richness, (B) Shannon Diversity, and (C) Faith's Phylogenetic Diversity from lactating (purple; n = 25) and non-lactating (orange; n = 30) devils.Results of statistical analysis can be found in Additional file 1: TableS5.Visualisation was performed using ggplot2 (v 3.4.2).

Figure 4 .
Figure 4. Differences in microbial community composition observed between reproductive status and location.Nonmetric Multidimensional Scaling (NMDS) plots of (A) unweighted, and (B) weighted UniFrac distances.Colour represents reproductive status while shape indicates location.Results of statistical analysis can be found in Additional file 1: TableS6.Visualisation was performed using ggplot2 (v 3.4.2).

Figure 5 .
Figure 5. Differentially abundant amplicon sequence variants (ASVs) between reproductive groups as identified by DESeq2 analysis.Each point represents a single differentially abundant ASV, grouped according to genus.The comparison was made between lactating versus non-lactating pouches, hence a positive Log 2 FoldChange indicates a higher prevalence in lactating pouches and a lower prevalence in non-lactating pouches.All ASVs visualised above were found to be significantly different between reproductive groups using a Kruskal-Wallis test.Visualisation was performed using ggplot2 (v 3.4.2).
All research involving the capture, handling, and collection of wild samples was carried out in accordance with guidelines under the University of Sydney Animal Care and Ethics Committee and all collection protocols were approved by the University of Sydney Animal Care and Ethics Committee approval number 2019/1562.As the subjects of the sample collection were wild animals no informed consent was required.Samples collected for this study were collected by the Tasmanian Department of Natural Resources and Environment's Save the Tasmanian Devil Program under their scientific permit.