Study sites and habitat
Captive
Captive SHNWs accessed for this study were held by an experienced wildlife carer in the Adelaide Hills, approximately 10km south-west of Adelaide, South Australia. These were rescue animals brought in from across the state, with most being hand-reared from a young age. The animals were housed in separate, custom-built enclosures containing burrows, and fed Barastoc Complete Performer (primarily cooked and flaked barley and lupins, lucerne and cereal chaff). The diet was bolstered in the wintertime with sunflower, oats, carrots, and sweet potato. Animals had access to water ad-lib.
Degraded habitat - Kooloola Station and Brookfield Conservation Park
Wild SHNWs from degraded habitat were located at Kooloola Station and Brookfield Conservation Park (Figure 1). Both sites were historic sheep grazing properties and are located approximately 2.5 hrs drive north east of Adelaide, in South Australia’s Murraylands region. The annual rainfall in this area is ~270mm (Bureau of Meteorology -- BOM, 2018). The area consists of remnant mallee eucalypt woodland with a chenopod and grassland understory composed of a variety of introduced weeds, native grasses and chenopods. The diet of wombats in this area is made up largely of introduced weed species including: thread iris (Moraea setifolia), wards weed (Carrichtera annua) and burr medic (Medicago minima), along with remnant chenopods and a few native grass species [31]. Except immediately after rain, no free water is available at these sites.
Intact habitat - Wonga Station
Wild SHNWs from intact habitat were located at Wonga Station, a sheep grazing property of ~130,000 acres, located about 3.5 hours drive north east of Adelaide, in South Australia’s Murraylands region (Figure 1). This area receives ~290 mm of rainfall annually (BOM, 2018). The vegetation on this property is predominantly composed of intact, native grassland habitat interspersed with remnant mallee eucalypt woodland. Active weed management is undertaken at this site. Native grasses and forbs, the dominant component of wombat diet in this region, flourish at this site and include many Stipia, Hyalosperma, Silene, Rytidiosperma and Sida spp, and a wide variety of Chenopodiaceae (saltbush - Atriplex, Enchylaena and Rhagodia spp, and bluebush - Maireana spp) among others [31]. Except immediately after rain, no free water is available at this site.
Sample collection
Once faeces have been deposited by an animal, the microbes inside can continue to grow and distort the true representation of the microbial community as it was inside the host. To counteract this potentially confounding variable, faecal samples are best collected fresh and preserved. Freezing is a commonly used preservative method, but is difficult to use in a fieldwork context. Instead, we opted for preservation of faecal samples by immersion in 95% ethanol, which has been demonstrated to reliably preserve faecal microbial community composition at room temperature [32]. We collected samples in the evening/early morning when SHNWs are most active, offering the best opportunity to collect fresh faecal samples.
The entrance and soil mound area around SHNW warrens were searched for fresh faecal pellets. Old pellets and fresh pellets were qualitatively distinguished by gently squeezing with a freshly gloved hand. Once found, fresh pellets were placed on a piece of aluminium foil and cut in half using the cutting edge of a spatula. The pellet cores were extracted using the opposite end of the spatula and placed in 2 mL plastic screw-top tubes prefilled with 1.5 mL of 95% ethanol, before being shaken vigorously to ensure mixing of sample with ethanol. Gloves were changed between samples, and spatulas were decontaminated using 5% bleach (sodium hypochlorite) followed by an ethanol wash.
Samples from Kooloola Station were collected in the evening of 13 March 2019; Brookfield samples were collected in the morning of 14 March 2019 (mean overnight temperature of 7.2°C; Eudunda weather station, Bureau of Meteorology). Samples from Wonga Station were collected overnight and into the morning of 19 June 2019 (mean overnight temperature of 1°C; Eudunda weather station, Bureau of Meteorology). Scat samples from captive animals were collected between 13 May 2019 and 5 June 2019, with some individuals being sampled at up to three different time points. A soil sample from each wild study site was also collected for analysis.
DNA extraction
All DNA extractions were performed in freshly decontaminated Perspex hoods in a pre-PCR laboratory to prevent contamination with amplicons [33]. DNA was extracted using the QIAGEN DNeasy PowerSoil kit (formerly MO BIO PowerSoil) according to the manufacturer’s protocol. To reduce contamination, all buffers and tubes needed for the various steps were aliquoted prior to opening any sample tubes. Faecal samples were centrifuged at 10,000 g for 5 minutes before pouring the ethanol off. Because SHNW faeces are very dry, only ~150 mg from each sample was used. To minimise batch effects samples in extraction groups were randomized, and to account for laboratory related contamination extraction blank controls from each extraction group was included and carried through to DNA sequencing [33].
Amplicon library preparation, quantification, and DNA sequencing
All samples were PCR-amplified and uniquely barcoded for High-Throughput Sequencing (HTS) using primers targeting the V4 region of the bacterial 16S rRNA gene [34]. We used forward primer: 515F (AATGATACGGCGACCACCGAGATCTACACTATGGTAATTGTG-TGCCAGCMGCCGCGGTAA) and barcoded reverse primer 806R (CAAGCAGAAGACGGCAT-ACGAGATnnnnnnnnnnnnAGTCAGTCAGCCGGACTACHVGGGTW TCTAAT) – the 12 n’s represent unique barcode sequences. The PCR reactions were prepared in a pre-PCR laboratory in a 5% bleached-cleaned and UV irradiated hood. Single reactions [35] of 2.5 µL X10 HiFi buffer, 0.1 µL Platinum™ Taq DNA Polymerase (ThermoFisher), 19.2 µL dH2O, 0.2 µL 100 mM dNTP mix, 0.5 µL each of 10 µM forward and reverse primer and 1 µL DNA. DNA was amplified using an initial denaturation at 94 °C for 3 min, followed by 35 cycles of denaturation at 94 °C for 45 sec, annealing at 50 °C for 1 min, elongation at 68°C for 90 sec, with final adenylation for 10 min at 68 °C, in line with the Earth Microbiome Protocol [36].
Gel electrophoresis was carried out for each PCR reaction on a 3.5% agarose gel to ensure the samples contained library constructs of the desired length (~390 bp). For each sample, 1 µL amplified DNA was mixed into 199 µL Qubit® working solution (diluted Qubit® dsDNA HS Reagent 1:200 in Qubit® dsDNA HS Buffer) and quantified using a Qubit® 2.0 Fluorometer. Samples were pooled equimolar and cleaned using AxyPrep™ (Axygen) following the manufacturer’s instructions. Because negative controls contained little DNA, they were pooled separately and spiked into the final pool at a flat volume [33]. The final pool was quantified and quality checked using an Agilent TapeStation. DNA sequencing was performed on an Illumina MiSeq (v2, 2 x 150 bp) at SAHMRI (South Australian Health and Medical Research Institute).
Data processing and statistical analyses
DNA sequencing data were processed and analysed using QIIME2 v2020.2 [37]. An interactive Jupyter notebook containing all the QIIME2 code used is available (SI File 1). Forward reads (R1) were imported into the QIIME2 and denoised into amplicon sequence variants (ASVs) using the deblur [38] plugin with a trim length of 150 bp. Representative sequences were assigned taxonomy using the QIIME2 feature-classifier plugin (naive bayesian approach) on the pre-trained SILVA [39] 132 V4 region classifier [40]. A phylogenetic tree was created using the SATé-enabled phylogenetic placement (SEPP) technique [41] within the fragment-insertion QIIME2 plugin. Alpha diversity (Observed OTUs and Faith’s Phylogenetic diversity [42]) and Beta diversity (weighted [43] and unweighted [44] UniFrac metrics) were calculated using the QIIME2 diversity plugin with a rarefaction depth of 36,346 sequences per sample. Tests for differentially abundant microbes between populations were performed using ANCOM [45]. SINA (SILVA Incremental Aligner) [46] was used to search specific ASVs against the SILVA 138 database. We used SourceTracker2 (https://github.com/biota/sourcetracker2) [47] to assess whether microbes from the soil were a source for faecal samples. Figures were constructed using ggplot2 [48] in R Studio [49]. The map was generated using Geoplaner (www.geoplaner.com) with data from OpenStreetMaps (wiki.openstreetmap.org). Venn diagrams were created using InteractiveVenn (www.interactivenn.net).
Removal of outlier samples
We observed 7 samples from Brookfield Conservation Park that appeared to be outliers in our dataset. PCoA of unweighted UniFrac distances clearly separated these samples from other wild and captive samples (SI Figure 1). Additionally, the diversity in these outlier samples was higher than other Brookfield samples (mean of 1,323 observed ASVs vs. 831). The taxonomic composition of these outlier samples was also substantially different to other samples (SI Figure 2). An interactive view at different taxonomic levels can be observed by dragging SI File 2 into the https://view.qiime2.org/ webpage. The most stark difference between these outlier samples and the other samples in the dataset is the loss of Spirochaetes, which was present in all other wild samples with a mean relative abundance of 11.6% and captive samples at 5.1%. A possible cause for these differences could be taphonomy—perhaps these faecal samples were older than they seemed, which could have resulted in a shift of observed community composition. We chose to exclude these outlier samples from subsequent analyses.
Frequency-based filtering to mitigate cross-sample contamination
Samples undergoing DNA extraction and library preparation in the same batch can be affected by cross-sample contamination, whereby cells (from DNA extraction) or DNA (from DNA extraction and/or library preparation) can cross from a sample to another [33,50]. Such cross-sample contamination can therefore yield false-positive detection of microbes within samples, with samples of higher microbial biomass being more likely to spillover [50]. Because we randomized sample order within DNA extractions and library preparations (to reduce batch effects), there is a possibility of spillover between faecal samples from different populations. This would confound attempts at identifying ASVs that are unique to a particular population (i.e. Figure 3B, Figure 5). To try and mitigate the effect of cross-sample contamination for population-specific ASV analysis, we applied a conservative frequency-based filtering approach, whereby population-specific ASV tables were filtered to remove ASVs with a relative abundance of <0.00015%). This threshold was based on the frequency of the most abundant ASV in the dataset (3c44df3672100a011a334b67eea24366), which was found in all wild samples with a total frequency of 278,742 reads (5.8% of total reads in wild samples). In contrast, this ASV was only found in the captive-only table with a total frequency of 131 reads (0.0001% of the captive table abundance). Assuming that the most abundant ASV would result in the greatest number spillover events, setting a threshold based on this ASV (0.00015%) would account for the spillover of less abundant ASVs. See the interactive Jupyter notebook for more details and the exact code ran (SI File 1). The issue of cross-samples contamination in microbiome studies clearly deserves further research.