Oral Microbiome Metabarcoding in Two Invasive Small Mammals from New Zealand

: All multicellular organisms host a wide diversity of microorganisms in and on their bodies, which are collectively known as their microbiome. Characterising microbial communities that inhabit di ﬀ erent body niches in wild animals is critical to better understand the dynamics of microbiome diversityand its functional signiﬁcance. The current study is the ﬁrst to apply massively parallel sequencing of 16S rRNA to characterise the microbial diversity and functional content of oral microbiota in two of New Zealand’s most important invasive mammals, the omnivorous common brushtail possum ( Trichosurus vulpecula ) and the carnivorous stoat ( Mustela erminea ). In total, strains of bacteria belonging to 19 di ﬀ erent phyla, 27 classes, 52 orders, 103 families, 163 genera and 51 known species were identiﬁed from the oral cavities of the study species. Strains of the phyla Proteobacteria, Firmicutes, Bacteroidetes, Fusobacteria, and Actinobacteria dominated the core oral microbial diversity in both species, while other taxa were comparatively less abundant. Despite invasive populations typically demonstrating limited genetic variation, intraspeciﬁc variation of the core bacterial taxa in the oral microbiota was considerable. This suggests that a complex interaction between genetic, physiological, and environmental factors determines the diversity of the study species’oral microbiome. phylum TM6 Dependentiae ), thermophilic Thermotogae , Armatimonadetes (OP10) and Acidobacteria from cavities


Introduction
Microbes have dominated Earth's evolutionary landscape for billions of years and have developed a plethora of biochemical adaptations to exploit all conceivable ecological niches on the planet. Some of these microbial metabolic repertoires ultimately became integrated into multicellular organisms [1,2], and have been an essential component in organismal evolution and development ever since [3].
All eukaryotic organisms host a variety of bacteria, archaea, fungi, protists, and viruses in and on their bodies. The combined number of these microorganisms and their genome size, known collectively and genomic adaptations to new environments take place on a comparatively short evolutionary time scale, implying that evolutionary mechanisms other than those originating from the slowly-evolved host genome may be involved [12,42].
In the current study, we used metabarcoding of 16S rRNA by means of massively parallel sequencing to characterise the taxonomic diversity and functional content of oral microbiota in the brushtail possum and the stoat in New Zealand and investigated the extent of intra-and inter-species variation in oral microbial diversity in the subject animals.
This study may serve as a starting point for more elaborate and robust analyses of oral microbiomes in the context of biological invasions. Such knowledge holds unprecedented therapeutic, environmental, and socio-economic promises that reach beyond the realm of wildlife genomics [43].

Animal Ethics Statements
The Lincoln University, Animal Ethics Committee approved oral sample collection and animal handling (SOP 29-14 "Protocol for sedation of possums"; SOP 33-14 "Protocol for anaesthesia of stoats"). All necessary steps were taken to minimise the pain and suffering of the test subjects.

Sample Collection and Genomic Library Preparation
Five adult brushtail possums (three females and two males) and five stoats (three females and two males) were live-captured from different populations on Bank Peninsula and West Coast, New Zealand. The animals were immediately transferred to the wild animal husbandry, Johnstone Memorial Laboratory at Lincoln University. After arrival at the facility, the health status of each animal was assessed by the Lincoln University wildlife veterinary team. To reduce the effects of diet on oral microbiome diversity, individuals of each species were fed the same diet prior to microbiome sample collection. The microbiome samples were collected between 9 a.m. and 12 p.m. from fasted subject animals within 36 h of arrival at the facility.
Stoats were anaesthetised using halothane and isoflurane inhalation protocols [44] for a maximum of 2 min [45], whereas the brushtail possums were lightly sedated by intravenous administration of an appropriate dosage of tiletamine-zolazepam (Zoleti ® ) [46]. After complete immobilisation was achieved, each animal's outer surfaces of the gingiva in the upper and lower dental arches, the inner cheek lining, the soft and hard palates, the surface of tongue, and the floor of the mouth were carefully swabbed using a sterile cotton swab. Special attention was paid to ensuring that sample collection was consistent between individuals. The cotton swabs were air-dried in a sterile cabinet for a maximum of 5 min, transferred to a sterile cylindrical plastic pouch, and stored at −80 • C until DNA was extracted, within a week.The health condition of each animal was closely monitored during and after the sampling procedure until full recovery was observed.
Metagenomic DNA of high molecular weight was extracted from the buccal swabs [47], and the hypervariable regions (V1-V2) of the 16S rRNA were amplified using the primer pair 8F (5 -TCG TCG  GCAGCG TCA GAT GTG TAT AAG AGA CAG AGA GTT TGA TCC TGG CTC AG-3 ) [48] and  338R (5 -GTC TCG TGG GCTCGG AGA TGT GTA TAA GAG ACA GGC TGC CTC CCG TAG GAG  T-3 ) [49]. Primer sequences contained overhang adapters appended to the 5 end for compatibility with Illumina sequencing platforms. Each PCR reaction contained 12.5 ng of template DNA, 0.2 µM of each primer and 2 × KAPA HiFi HotStart Ready Mix (KAPA Biosystems, Wilmington, NC, USA). Each batch of PCR reactions included three positive and three negative control reactions. The thermal profile for PCR amplifications started with an initial denaturing step at 95 • C for 3 min, followed by 25 cycles of denaturing at 95 • C for 30 s, annealing at 55 • C for 30 s, and a 30 s extension at 72 • C, followed by a final 5 min extension at 72 • C. PCR amplicons were purified using the AMPure XP reagent (Beckman Coulter, Indianapolis, IN, USA).This was followed by eight cycles of barcoding PCRs, each with an initial denaturing step at 95 • C for 3 min, annealing at 55 • C for 30 s, 30 s extension Demultiplexed sequences were processed using the Qiime2 v11.2017 package [50]. Illumina sequencing adapters, contaminants, and over-represented sequences were filtered using the Qiime2 Cutadapt v1.18 plugin [51]. The Qiime2 DADA2 plugin [52] was used to remove low-quality sequences (Phred score < 25), detect potential chimeric reads, and de-noise demultiplexed sequences into unique sequence features. To assign a taxonomy rank to the sequence features, the Qiime2 "blast-consensus feature classifier plugin" was used, with the number of best-accepted hits set to four, the minimum percentage identity to 70%, and the p-value to a maximum of 0.001. Quality-filtered sequences were searched against the Greengenes v13_8, 99% OTUs full-length database [53]. Sequences that were unassigned or lacked taxonomy IDs in the NCBI taxonomy database (https://www.ncbi.nlm.nih.gov/taxonomy) were discarded. The taxonomy ID of each consensus match was extracted using an in-house Python3 script, and a circular phylogenetic tree was reconstructed using the PhyloT online server [54] and visualised in ITol v.33 [55].
To estimate the phylogenetic distance necessary for computing some diversity indices, an approximate maximum-likelihood phylogenetic tree was constructed using the package FastTree2 [56]. Before running the core diversity script in Qiime2, the sequence counts across different individuals were randomly subsampled to an estimated equal rarefaction depth of 79,690 sequences. Rarefaction at this sequencing depth resulted in the maximum number of sequences (67.7% of the total sequences) being retained in nine individuals. To minimise the bias in diversity metrics estimation, one stoat with a lower number of sequences than the above threshold was discarded from the diversity analysis.
The statistical significance for differences in α-diversity (evenness and Faith phylogenetic indices) and β-diversity (unweighted UniFrac) between two species were tested by performing group-wise non-parametric Kruskal-Wallis [65] and permutational multivariate analysis of variance (PERMANOVA) tests [66].
To find modules of co-occurring bacterial taxa in each species, the SCNIC package (https://github.com/shafferm/SCNIC) was used to create a sparse positive co-occurrence correlation network of different bacterial taxa by setting the minimum pairwise correlation value to 0.30 (R = 30), and the network files were visualised in Gephi v0.9.2 [67].
Bacterial taxa whose abundance differed between two species were identified using the Gneiss balances method [68]. First, a correlation-clustering balance tree was constructed for each species. In the correlation-clustering algorithm, an unsupervised Ward's hierarchical clustering algorithm [69] groups different bacterial taxa interactively based on their co-occurrence in each species' oral microbiota. Then, a multivariate linear regression analysis [70] was performed independently on each node (balance) of the resulting tree and tested for significant changes in bacterial abundance between pairs of species. In the regression analysis, the description (animal ID) and taxonomic infraclass (marsupial versus placental) were selected as explanatory variables, and taxon abundance was set as the response variable. To evaluate the explanatory power of a single covariate, a leave-one-variable-out method was applied. The variable whose removal caused the most significant changes in the fitted model was reported as having the strongest effects on oral microbiome diversity. Overfitting of the regression model was Diversity 2020, 12, 278 5 of 18 monitored by dividing the data into ten partitions and evaluating the prediction accuracy of the model built based on nine partitions on the remaining partition.
The PICRUSt2 [71] pipeline was executed to predict metabolic pathways for the different bacterial taxa present in each species' oral cavities. Alternatively, the biochemical functions of differentially abundant taxa between pairs of species were predicted using the BioCys online server [72]. Functional enrichments of the oral microbial communities for different biochemical functions were tested using Microbiome Analyst tools [73]. Raw sequences generated for this study are available on the NCBI SRA database under BioProject (accession number PRJNA61279).

Results
None of the animals showed any signs of major acute or systematic pathology during physical examination. Only minor gum inflammation, typical of wild populations, was observed. Negative control reactions in all PCR setups did not produce any products. Moreover, genomic DNA libraries prepared from negative control reactions failed quality checking prior to Illumina sequencing and were discarded.
The sequencing run yielded 1219,643 paired-end, quality-filtered sequences. Each buccal swab produced 121,964 sequences on average. Brushtail possum oral swabs yielded more raw sequences compared with stoat oral swabs. The DADA2 plugin collapsed the raw sequences into 1083 unique sequence features. The blast-consensus feature classifier assigned a taxonomy rank to 590 sequence features in brushtail possums and to 549 sequence features in the stoats, with a mean posterior confidence level of 92%.
Principal coordinates analysis of the oral microbiota based on Bray-Curtis, Jaccard's, and unweighted Unifrac metrics demonstrated that the oral microbiota of the brushtail possums clustered densely together, while stoat microbiota showed higher intra-individual variation and clustered more loosely. This pattern is less evident in weighted UniFrac metrics where the presence of low-abundance taxa most likely confounded the clustering pattern ( Figure 1).    The linear regression summary in Gneiss indicated that approximately 45% of the variance in oral bacterial diversity were explained by our model. Taxonomic infraclass (brushtail possum versus stoat) had significantly more influence on bacterial diversity than intraspecies variation (Rdif of 0.188 versus 0.0796). The 10-fold cross-validation confirmed that within-model error was higher than prediction accuracy in all replicates, and that model overfitting was rejected.
Bacterial taxon abundance in one specific balance in the correlation tree differed significantly between possums and stoats (false discovery rate corrected coefficient p-values of 0.001 and 0.030 for intercept and a taxonomic rank, respectively). This balance, and its descending leaf nodes, mainly consisted of strains of Streptococcus spp., Leptotrichia spp., Heamophilus spp., and unidentified strains belonging to Pasteurellaceae and Aerrococcaceae families.  In the modular network analysis, graph density was similar between brushtail possums and stoats (0.12). The average degree of node connectivity estimated for stoats was higher than that for brushtail possums (16. Metabolic pathway analysis of oral microbiome diversity revealed that the microbial taxa in the oral cavities of both species are functionally enriched for biochemical pathways involved in the biosynthesis of various metabolites (e.g., vitamins, aminoamides and organic carbon compounds), starch and sucrose metabolism, steroid hormone degradation, and antimicrobial activities (e.g., streptomycin biosynthesis) ( Table 2). Functional content analysis of the bacterial taxa that differ in abundance between the two species illustrated that the majority of these bacteria are involved in the biosynthesis of different metabolites, biochemical degradation and assimilation, energy production, respiration, and detoxification.

Discussion
Recent advances in comparative genomics, biochemistry, and developmental biology have made it possible to demarcate species boundaries with an unprecedented level of accuracy. However, the interdependence between hosts and a wide diversity of microorganisms that reside in and on an individual has provided a new set of challenges in defining what constitutes an individual.
In the current study, we have described and discussed the taxonomic diversity and functional content of the oral microbiome of brushtail possums and stoats in New Zealand. The number of sequence features identified in the oral cavities of both species resembles that reported for humans [74]. However, the sequence features that were assigned the taxonomy rank of genus or species were A B

Discussion
Recent advances in comparative genomics, biochemistry, and developmental biology have made it possible to demarcate species boundaries with an unprecedented level of accuracy. However, the interdependence between hosts and a wide diversity of microorganisms that reside in and on an individual has provided a new set of challenges in defining what constitutes an individual.
In the current study, we have described and discussed the taxonomic diversity and functional content of the oral microbiome of brushtail possums and stoats in New Zealand. The number of sequence features identified in the oral cavities of both species resembles that reported for humans [74]. However, the sequence features that were assigned the taxonomy rank of genus or species were limited, and potentially represented a high level of concealed microbial diversity that is only broadly known to science.
The oral microbiota in stoats and brushtail possums shared common features with those reported from other mammals. In both species, the oral microbiota were dominated by a small number of taxa, while other taxa were comparatively rare. The abundant bacterial taxa have been ubiquitously reported in different herbivorous, carnivorous, and omnivorous mammal species (Figure 3) [23][24][25][26][27][28][29][30].
Diversity 2020, 12, x FOR PEER REVIEW 10 of 17 that act as fermentation chambers [80]. Firmicutes bacteria in the oral cavity of brushtail possums can initiate an early digestion step that ameliorates overall energy uptake efficiency. The absence of enrichment for biochemical pathways involved in vitamin B6 (pyridoxine) biosynthesis in the brushtail possum's oral microbiome is also consistent with a diet rich in plants. Similarly, lack of metabolic pathways to biosynthesise essential aminoacids, such as valine, leucine, and isoleucine,in the stoats' oral microbiota may reflect a high intake of animal protein in the diet. Our results are amongst the first to report the presence of culture-independent candidate phylum TM6 (Dependentiae), thermophilic Thermotogae, Armatimonadetes (OP10) and Acidobacteria from the oral cavities of wild mammals. Reduced genome size and a limited biochemical repertoire in TM6 (Dependentiae) bacteria is indicative of an obligate endosymbiotic or parasitic lifestyle in the oral cavity [81]. Strains of Thermotogae, Armatimonadetes (OP10), and Actinobacteria have been reported from the microbiota in the gut of pregnant rats (Rattus norvegicus) [82], in the human oral cavity [83], and in the intestines of rice frogs (Fejervaryalimnocharis) [84]. The origin and functional significance of these taxa in the oral microbiota of the study species remain to be discovered.
The considerable intraspecific variation in the oral microbial diversity of each study species (Figures 4 and 5) suggests the influence of variables other than diet that were more difficult to control. These may include genetic and physiological factors, as well as differences in the environment prior to capture. These need to be taken into consideration to better understand the complexity of microbiome diversity [85][86][87][88].  Diet has been suggested as one of the major determinants of oral microbial diversity [75][76][77][78]. However, consuming a similar diet is not necessarily a good predictor of oral microbiome diversity. For example, the microbiome diversity of pig (omnivore), cat (carnivore) and wallaby (herbivore) is quite similar, whereas the closely related chimpanzees and bonobos show considerable differences. In the present study, oral microbiome samples were obtained from several different niches within the oral cavity. Previous studies have shown that each niche may harbour a unique microbial community, and direct comparisons with such studies are thus limited by the fact that we pooled samples from multiple niches [79].
Our data suggest a potential link between diet and the presence of some bacterial taxa, and the biochemical functions they may perform in the oral cavity. For instance, a plant polysaccharide metabolising taxon, Treponema sp., was completely absent from the oral microbiota of the carnivorous stoat. High diversity of Firmicutes in brushtail possum oral microbiota most likely reflects the presence of staple plant-based polysaccharides in their omnivorous diet (Figures 4 and 5). The distal digestive tract in this species consists of a well-developed caecum and an enlarged proximal colon that act as fermentation chambers [80]. Firmicutes bacteria in the oral cavity of brushtail possums can initiate an early digestion step that ameliorates overall energy uptake efficiency. The absence of enrichment for biochemical pathways involved in vitamin B6 (pyridoxine) biosynthesis in the brushtail possum's oral microbiome is also consistent with a diet rich in plants. Similarly, lack of metabolic pathways to biosynthesise essential aminoacids, such as valine, leucine, and isoleucine, in the stoats' oral microbiota may reflect a high intake of animal protein in the diet.  Microorganisms in the alimentary tract represent a subset of the microbial diversity found in the environment. The exact mechanisms by which the host and its resident microbiome differentially select specific microbial diversity and functional content from the environment is unknown. However, earlier studies have shown that microbial diversity can be strongly predicted by the genetic make-up of the host species [15,[89][90][91][92][93][94][95][96][97]. For examples, polymorphisms in the immunity-related gene interleukin-2 [90] and genes involved in core metabolic pathways, such as nucleotide-binding oligomerisation domain-containing protein 2 [91] and fucosyltransferase 2 [94], have been linked to a specific microbiome structure in the host species.
Despite considerable intra-and interspecific variations in oral microbiota, both species nonetheless display significant similarities in the metabolic functions of their oral bacterial communities. The observed pattern substantiates the idea that different hosts may assemble taxonomic diversity and functional content of their microbiome in different ways, but with a high level of taxonomic and functional redundancy that prevents a host's dependency on a limited number of bacterial taxa [98][99][100]. Sharpton [101] suggested that the functional content of the microbiome can be more heritable than its taxonomic diversity.
Species in the current study have been subjected to intensive population control measures using a combination of toxins and trapsfor more than 60 years. While numerous studies have underlined the critical role that gut microbiota play in drug and xenobiotic substance metabolism [102], the role of microbial communities in the proximal alimentary tract is yet to be understood.
As conservation efforts to eradicate invasive species intensify around the globe, focusing on the microgenome rather than the slowly-evolving host genome in genetically deprived invasive species presents an under-appreciated means of studying evolutionary responses to altered environments over a short time-scale. Future comparative studies should focus on studying microbiome diversity in a larger sample of individuals from both native and invaded ranges, and investigate the significance of bacterial communities in modulating the host's physiological and behavioural responses. Such studies would benefit from a considerably more complex study design than the one described here, as they would require that the roles of biotic and abiotic factors be detangled, while controlling for the confounding effects of the host's diet, physiological condition, sex, age, and whether individuals were sampled in the wild or had spent some time in captivity. Our results are amongst the first to report the presence of culture-independent candidate phylum TM6 (Dependentiae), thermophilic Thermotogae, Armatimonadetes (OP10) and Acidobacteria from the oral cavities of wild mammals. Reduced genome size and a limited biochemical repertoire in TM6 (Dependentiae) bacteria is indicative of an obligate endosymbiotic or parasitic lifestyle in the oral cavity [81]. Strains of Thermotogae, Armatimonadetes (OP10), and Actinobacteria have been reported from the microbiota in the gut of pregnant rats (Rattus norvegicus) [82], in the human oral cavity [83], and in the intestines of rice frogs (Fejervaryalimnocharis) [84]. The origin and functional significance of these taxa in the oral microbiota of the study species remain to be discovered.

Conclusion
The considerable intraspecific variation in the oral microbial diversity of each study species (Figures 4 and 5) suggests the influence of variables other than diet that were more difficult to control. These may include genetic and physiological factors, as well as differences in the environment prior to capture. These need to be taken into consideration to better understand the complexity of microbiome diversity [85][86][87][88].
Microorganisms in the alimentary tract represent a subset of the microbial diversity found in the environment. The exact mechanisms by which the host and its resident microbiome differentially select specific microbial diversity and functional content from the environment is unknown. However, earlier studies have shown that microbial diversity can be strongly predicted by the genetic make-up of the host species [15,[89][90][91][92][93][94][95][96][97]. For examples, polymorphisms in the immunity-related gene interleukin-2 [90] and genes involved in core metabolic pathways, such as nucleotide-binding oligomerisation domain-containing protein 2 [91] and fucosyltransferase 2 [94], have been linked to a specific microbiome structure in the host species.
Despite considerable intra-and interspecific variations in oral microbiota, both species nonetheless display significant similarities in the metabolic functions of their oral bacterial communities. The observed pattern substantiates the idea that different hosts may assemble taxonomic diversity and functional content of their microbiome in different ways, but with a high level of taxonomic and functional redundancy that prevents a host's dependency on a limited number of bacterial taxa [98][99][100]. Sharpton [101] suggested that the functional content of the microbiome can be more heritable than its taxonomic diversity.
Species in the current study have been subjected to intensive population control measures using a combination of toxins and trapsfor more than 60 years. While numerous studies have underlined the critical role that gut microbiota play in drug and xenobiotic substance metabolism [102], the role of microbial communities in the proximal alimentary tract is yet to be understood.
As conservation efforts to eradicate invasive species intensify around the globe, focusing on the microgenome rather than the slowly-evolving host genome in genetically deprived invasive species presents an under-appreciated means of studying evolutionary responses to altered environments over a short time-scale. Future comparative studies should focus on studying microbiome diversity in a larger sample of individuals from both native and invaded ranges, and investigate the significance of bacterial communities in modulating the host's physiological and behavioural responses. Such studies would benefit from a considerably more complex study design than the one described here, as they would require that the roles of biotic and abiotic factors be detangled, while controlling for the confounding effects of the host's diet, physiological condition, sex, age, and whether individuals were sampled in the wild or had spent some time in captivity.

Conclusion
Our study highlights the critical role that a combination of genetic, physiological, and environmental factors plays in shaping oral microbiome diversity. The integrated anatomical, metabolic and immunological fitness of the host and its associated microbiome depend on a combination of concurrent selective and random changes in the genomes of both the host and its microbiome in response to changing environmental conditions.