Uncovering the microbiome of invasive sympatric European brown hares and European rabbits in Australia

Sample collection

Faecal samples were collected at necropsy from nine wild European brown hares (Lepus europaeus) and 12 wild European rabbits (Oryctolagus cuniculus) of both sexes (Table S1). Although caecal contents may more accurately reflect the microbiota involved in digestion in these species, these samples were not collected for this study. The hares and rabbits used in this study were mostly adults or young adults, they were outwardly healthy, and were shot as part of routine vertebrate pest control operations in Mulligans Flat, ACT (−35.164, 149.165), between January and September 2016. Sample collection was approved by Barry Richardson, Oliver Orgill, and the ACT Parks and Conservation Service rangers of Mulligan’s Flat Woodland Sanctuary. Mulligans Flat Nature Reserve comprises 781 hectares of box-gum grassy woodland with a sparse to moderately dense ground cover of grasses, herbs and shrubs. The reserve is fenced, preventing immigration and emigration of most animals. Rabbits and hares were shot from a vehicle using a 0.22-caliber rifle targeting the head or chest. Faecal pellets from the descending colon were collected during post-mortem and stored at −20 °C. All sampling was conducted according to the Australian Code for the Care and Use of Animals for Scientific Purposes as approved by CSIRO Wildlife and Large Animal Ethics Committee (approvals #12-15 and #16-02).

Genomic DNA extraction

Genomic DNA was extracted from 50 mg of faecal pellets using the DNeasy Blood and Tissue kit as per manufacturer’s instructions and described in detail at https://www.protocols.io/view/genomic-dna-extraction-from-animal-faecal-tissues-6bbhain (Qiagen, Chadstone Centre, Victoria, UK). Hare and rabbit samples were processed separately and reagent-only controls (ROC) were included with each set of extractions. All samples were quantified and evaluated for integrity using the Qubit dsDNA Broad Range assay kit and the Qubit Fluorometer v2.0 (Thermo Fisher Scientific, Waltham, MA, USA) and Nanodrop ND-1000 (Thermo Fisher Scientific, Waltham, MA, USA). Genomic DNA samples were stored at −20 °C.

Illumina library preparation, sequencing and bioinformatics analysis

We amplified the V3–V4 region of the 16S rRNA gene (~460 bp) from genomic DNA and ROC according to the Illumina protocol for 16S Metagenomic Sequencing Library Preparation using a dual-indexing strategy, with modifications (Illumina, 2013). Briefly, initial PCR reactions (25 μl) were performed for each sample (including ROC) using 2× Platinum SuperFi PCR Master Mix (Thermo Fisher Scientific, Waltham, MA, USA), 0.5 μM of forward (5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG-3′) and reverse primer (5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC-3′) (overhang adapter sequence highlighted in bold) and 4.6 ng of genomic DNA as described at https://www.protocols.io/view/library-preparation-protocol-to-sequence-v3-v4-reg-6i7hchn. A ‘no template control’ (NTC) was included during reaction setup. Cycling conditions were: 98 °C for 30 s, followed by 25 cycles of 98 °C for 10 s, 55 °C for 15 s and 72 °C for 30 s, with a final extension of 10 min at 72 °C. PCR products were confirmed by agarose gel electrophoresis before being purified using AMPure XP beads (Beckman Coulter, Indianapolis, IN, USA) at a 1× ratio as described previously (Illumina, 2013). Each amplicon was then dual-indexed with unique DNA barcodes using the Nextera XT index kit (N7XX and S5XX, Illumina, San Diego, CA, USA) for PCR-based barcoding (Illumina, 2013). For each 50 μl PCR reaction, we used five μl of each index (i7 and i5), and five μl of the first PCR product. Cycling conditions were as described above but limited to eight cycles. Final libraries (including ROC and NTC) were purified with AMPure XP beads (Beckman Coulter, Indianapolis, IN, USA), validated using the Tapestation D1000 high sensitivity assay (Agilent Technologies, Santa Clara, CA) and the Qubit High Sensitivity dsDNA assay kit (Thermo Fisher Scientific, Waltham, MA, USA), pooled at equimolar concentrations, and sequenced on an Illumina MiSeq using 600-cycle v3 chemistry (300 bp paired-end) at CSIRO, Black Mountain, Canberra.

Illumina fastq reads were analysed in QIIME 2-2019.7 software (Bolyen et al., 2019). Raw fastq reads were quality filtered (i.e. filtered, dereplicated, denoised, merged and assessed for chimaeras) to produce amplicon sequence variants (ASV) using the DADA2 pipeline via QIIME2 (Callahan et al., 2016). The DADA2 generated feature table was filtered to remove ASVs at a frequency less than two, and remaining ASVs were aligned against the SILVA_132 (April 2018) reference database using QIIME2 feature-classifier BLAST+ (Bokulich et al., 2018; Camacho et al., 2009; Quast et al., 2013). Fasta sequences from DADA2 were also aligned to the same reference database using BLASTn (2.2.28) through the command line interface (CLI) using an e-value cut-off of 1E−90 (Table S2). We processed and exported BLASTn outputs into QIIME2 to perform taxonomic analysis. All scripts used were deposited at https://github.com/SomaAnand/Hare_rabbit_microbiome. Raw sequence data were deposited in the sequence read archive of NCBI under accession number PRJNA576096.

The microbial diversity and richness in individual samples was estimated using alpha diversity metrics (Shannon index and observed OTU), and the variation in overall composition between rabbits and hares was measured using beta diversity metrics (Bray–Curtis dissimilarity) after rarefaction at a subsampling depth of 149,523 using the q2-diversity pipeline within QIIME2 (Faith, Minchin & Belbin, 1987; Anderson, 2001). This subsampling depth was chosen based on the sample (excluding one ROC (rabbit) and PCR NTC) with the lowest number of reads.

Nanopore library preparation, sequencing and bioinformatics analysis

The entire 16S rRNA gene (~1,500 bp) was targeted for sequencing using the Nanopore MinION sequencing platform. We amplified the 16S rRNA gene from all genomic DNA samples using the universal primers 27F and 1492R (Weisburg et al., 1991). PCR reactions (100 μl) were performed for each sample using 2× Platinum SuperFi PCR master mix (Thermo Fisher Scientific, Waltham, MA, USA), 0.4 μM each primer, and 11.5 ng genomic DNA as described at https://www.protocols.io/view/library-preparation-protocol-to-sequence-full-leng-6j6hcre. Cycling conditions were: 98 °C for 30 s, followed by 28 cycles of 98 °C for 10 s, 55 °C for 15 s and 72 °C for 40 s, with a final extension of 5 min at 72 °C. PCR products were confirmed by agarose gel electrophoresis before being purified with AMPure XP beads (Beckman Coulter, Indianapolis, IN, USA), validated using the Tapestation D1000 high sensitivity assay (Agilent Technologies, Santa Clara, CA, USA), and quantified using the Qubit High Sensitivity dsDNA assay kit (Thermo Fisher Scientific, Waltham, MA, USA). For library preparation we used the ligation sequencing kit 1D (SQK-LSK108) in combination with the native barcoding kit 1D (EXP-NBD103) (Oxford Nanopore Technologies, Oxford, UK) as per the manufacturer’s protocol, except 500 ng of input DNA per sample was used for end preparation (Hu & Schwessinger, 2018). For the addition of barcodes, 80 ng of end-prepped DNA was used as input, and barcoded samples were pooled in equimolar concentrations to obtain at least 400 ng of pooled DNA. We used a final library amount of 200 ng to obtain maximum pore occupancy on a MinION R.9.4.1 flow cell. Two MinION flowcells were used to run all samples and each flow cell had a mix of hare and rabbit samples to control for potential batch effects.

Nanopore raw reads in fast5 format were demultiplexed using deepbinner (Wick, Judd & Holt, 2018). Basecalling, adapter trimming, and conversion into fastq format was performed in Guppy 2.3.7 (Oxford Nanopore Technologies, Oxford, UK). BLASTn was used to locally align 16S sequences with a quality score higher than seven against the SILVA_132 reference database. Top hits were exported in Biological Observation Matrix format and imported into QIIME2 for subsequent analyses. Nanopore data was rarefied at a subsampling depth of 109,435 for diversity analysis based on the sample with the fewest Nanopore reads, in order to mimic the Illumina workflow. This rarefied dataset was then used to analyse the taxonomic results as per the Illumina workflow described above. Species level taxonomic classification was performed through the EPI2ME platform (Oxford Nanopore, Oxford, UK) (Table S2). Raw sequence data were deposited in the sequence read archive (SRA) of NCBI under accession number PRJNA576096.

Statistical analysis

We assessed whether the bacterial diversity of hares and rabbits was statistically different by using permutation-based statistical testing (PERMANOVA) via the QIIME diversity beta-group-significance pipeline within QIIME2 for both Illumina and Nanopore datasets (Anderson, 2001). Statistical significance in alpha diversity was estimated using the nonparametric Kruskal–Wallis test within the R package ‘ggpubr’ (Kruskal & Wallis, 1952). We also evaluated statistical differences between hare and rabbit faecal samples for each observed bacterial phyla using a combination of multiple tests. We performed analysis of variance (ANOVA) to estimate the variance of both populations and Student’s t-test to identify statistical differences between the means of both groups. Estimated p-values were corrected using the Benjamini–Hochberg method to control the False Discovery Rate for multiple hypothesis testing. Pearson correlation coefficients were calculated to assess the correlation between Illumina and Nanopore sequencing results using the scipy.stats package in Python.

Interrogating hare and rabbit faecal microbiota using 16S rRNA sequencing

We performed Illumina short-read sequencing of the 16S rRNA V3–V4 region and Nanopore long read sequencing of the entire 16S rRNA gene on faecal samples of wild hares and rabbits living sympatrically in a periurban nature reserve in Australia. Short-read sequencing of 24 samples through the Illumina MiSeq pipeline produced an aggregate of 14,559,822 reads greater than Q30 (mean 501,908 reads per sample, excluding ROC and NTC). Reads were converted to ASVs, which were assigned using the SILVA_132 reference database via BLAST+ in QIIME2. This produced 6,662 unique features at a frequency greater than two. Long-read sequencing of 21 samples (no ROCs or NTC due to lack of amplification) using the Nanopore MinION platform produced an aggregate of 6,544,770 reads (mean of 311,656 reads per sample) above a quality threshold of seven across two MinION runs. After aligning the reads against the SILVA_132 reference database using BLAST, 47,100 unique features with a frequency greater two were obtained.

The faecal microbiome composition varies significantly between sympatric wild rabbits and hares, and between individual rabbits

We conducted alpha and beta diversity analyses to assess the species richness, evenness, and composition both within individual animals (alpha diversity) and between rabbits and hares (beta diversity). The richness and evenness of bacterial species in individual animals was not significantly different between rabbits and hares, as measured by the Shannon index and Kruskal–Wallis test (Fig. S1). However, the beta diversity between rabbit and hare faecal samples differed significantly (p = 0.01 as measured by Bray–Curtis dissimilarity matrix and PERMANOVA) (Fig. 1). We also observed considerable variance in faecal microbial beta diversity between individual rabbits (Fig. 1). In contrast, the hare faecal samples examined all had a similar bacterial composition with relatively low variance between individuals. We observed these effects in both the Illumina and Nanopore 16S sequencing data sets (Fig. 1). We did not observe clear correlation between bacterial diversity and sex, reproductive status, or season of sample collection, although our statistical power was low due to the small sample size. Although there were more pregnant and/or lactating rabbits than hares, which may lead to biases in the observed bacterial diversity, this difference was not statistically significant (Fisher’s exact test, p = 0.14). Similarly, the variances of bodyweight for each population were not significantly different (F-test, p = 0.14).

Taxonomic classification of both the Illumina and Nanopore sequencing data sets identified Firmicutes and, to a lesser extent Bacteroidetes, to be the two dominant phyla in both hare and rabbit faecal samples (Fig. 2). On average, Firmicutes and Bacteroidetes together comprised more than 85% of the faecal microbiome in hares and rabbits and across both platforms. Hare faecal samples contained a significantly higher ratio of Firmicutes to Bacteroidetes compared to rabbit faecal samples (p = 0.015 as measured by Student’s t-test). In both hare and rabbit faecal samples, Firmicutes were predominantly represented by the families Ruminococcaceae, Christensenellaceae, Lachnospiraceae (genus Roseburia), Eubacteriaceae and Erysipelotrichaceaea, while Bacteroidetes were predominantly represented by the families Rikenellaceae and Barnesiellaceae (Fig. 3; Fig. S2). Within the Firmicutes and Bacteroidetes, the families Marinifilaceae (genera Odoribacter and Butryicimonas), Tannerellaceae (genus Parabacteroides) and Veillonellaceae were more abundant in hares compared to rabbits, while the families Clostridiaceae 1, Enterococcaceae, Planococcaceae, Lactobacillaceae, Peptostreptococcaceae and Bacillaceae were more abundant in rabbits, although not all families were present in all rabbits (Fig. 3; Fig. S2). Phylum Tenericutes (genus Anaeroplasmataceae) was present at a significantly higher relative frequency in rabbit faecal samples (p < 0.01 as measured by Student’s t-test). Phylum Verrucomicrobia (genus Akkermansia) was identified only in rabbit samples and not in hare samples, while phyla Lentisphaerae (genus Victivallis) and Synergistetes (genus Pyramidobacter) were identified only in hare samples. Other phyla such as Proteobacteria (genus Sutterella in both hares and rabbits, family Desulfovibrionaceae in hares only, and genus Enterobacter in rabbits only), Actinobacteria (genus Eggerthellaceae) and Patescibacteria (family Saccharimonadaceae) were present at low abundance (0.5–3%) in both rabbit and hare samples (Figs. 2 and 3; Fig. S2).

We were also able to perform taxonomic classification of Nanopore sequencing data to the species level, due to the longer read lengths associated with this sequencing technology. There was considerable variation in bacterial species between individual faecal samples and many reads remained unassigned at the species level. However, Ruminococcus albus clearly dominated the hare gastrointestinal microbiome, being present in all hare faecal samples at a mean relative frequency of 11.4% (Table 1). This bacterial species was also consistently detected in rabbit faeces, but at a much lower relative frequency (1.5%). Of the ten most abundant individual species, hare faeces were dominated by Ruminococcus sp, Bacteroides sp, Alistipes putredinis, and Roseburia hominis; in contrast, rabbit faeces were dominated by several Clostridium and Paraclostridium species, as well as Paeniclostridium sordellii, Ruminococcus albus, and Akkermansia glycaniphila (Table 1).

The NTC microbiome comprised the phyla Proteobacteria and Bacteroidetes at relative frequencies of 99% and 1%, respectively, from a total of 2,158 raw reads (Fig. S3). When analysing the ROCs, there was obvious genomic DNA contamination of both ROCs, with several phyla shared between faecal samples and ROCs (Fig. S3). Additional bacterial phyla were also detected only in ROCs, likely comprising the ‘reagent microbiome,’ including Plantomycetes, Dependentiae, Chloroflexi and Chlamydiae. The total number of raw reads in the rabbit ROC was 6,279 (compared to the mean 501,908 reads/sample), while contamination was more pronounced in the hare ROC, which produced 265,118 raw reads. Since the phyla-level distribution of reads in the hare ROC strongly reflected that seen in all hare samples, it was assumed that this contamination would not grossly alter the results, particularly because we were not dealing with low-biomass samples.

Geography alters the faecal microbiome of wild hares

We then compared the faecal microbiome of Australian wild hares to that of wild hares in their native range of Europe (sampled and sequenced in a previous study (Stalder et al., 2019)) to investigate the influence of geography on the gastrointestinal microbiome. Strikingly, phylum Spirochaetes was the third most abundant phylum in faecal samples of European origin (Fig. 4), yet it was only present at very low abundance from Australian hares and rabbits. Furthermore, phylum Patescibacteria was below our limit of detection (<0.5% relative abundance) in European hare samples, while in Australian hare samples it was present at a higher relative abundance (1.3%). Other phyla were detected in both European and Australian hare samples at similar relative abundances (Figs. 2 and 4).

Nanopore and Illumina sequencing platforms reveal a similar faecal microbiome independent of platform

We found a strong correlation (as measured by Pearson correlation coefficient) between Illumina and Nanopore sequencing results at the phylum (r = 0.945, p < 0.001), family (r = 0.947, p < 0.001) and genus (r = 0.923, p = <0.001) levels (Figs. 2 and 3; Fig. S2). Most notably, phyla Actinobacteria and Synergistetes were present at a higher relative abundance in the Illumina dataset, while phyla Tenericutes, Lentisphaerae, and to a lesser extent Cyanobacteria, were present at higher abundance in the Nanopore dataset (Fig. 2). At the family level, Thermoanaerobacteraceae, Bacteroidaceae, Clostridiales vadinBB60 group and Atopobiaceae were detected in the Illumina but not the Nanopore data, while Muribaculaceae and Flavobacteriaceae were identified only in the Nanopore data (Fig. 3; Fig. S2). Even at the genus level there was generally good agreement between datasets (Fig. 3). However, genera Butyricimonas, Mailhella, Ruminococcus 2, Erysipelatoclostridium and Sutterella were only detected in Illumina data and Parasutterella, Schwartzia, an unassigned Victivallaceae, and the Ruminococcaceae NK4A214 group were only identified in the Nanopore data.