Gut microbiota of the European Brown Hare (Lepus europaeus)

Diseases of the gastrointestinal tract due to changes in the bacterial flora have been described with increasing incidence in the European brown hare. Despite extensive demographic and phylogeographic research, little is known about the composition of its gut microbiota and how it might vary based on potential environmental or host factors. We analysed the intestinal and faecal microbiota of 3 hare populations by Illumina MiSeq 16S rRNA gene amplicon sequencing. The phyla and OTU abundance composition differed significantly between intestinal and faecal samples (PERMANOVA: P = 0.002 and P = 0.031, respectively), but in both sample types Firmicutes and Bacteroidetes dominated the microbial community composition (45.51% and 19.30% relative abundance). Intestinal samples contained an enrichment of Proteobacteria compared with faecal samples (15.71-fold change, P < 0.001). At OTU level, a significant enrichment with best BLAST hits to the Escherichia-Shigella group, Eubacterium limosum, Sphingomonas kyeonggiensis, Flintibacter butyricus and Blautia faecis were detected in intestinal samples (P < 0.05). In our statistical model, geographic location and possibly associated environmental factors had a greater impact on the microbiota composition than host factors. Population had a significant effect on the composition of abundant intestinal and faecal OTUs, and on the abundance of potential pathogenic bacteria of the family Enterobacteriaceae, regularly associated with intestinal dysbiosis in hares, in faecal samples. Our study is the first to describe the microbiota in brown hares and provides a foundation to generate hypothesis aiming to test the role of gut health in population fluctuations of the species.


Material and Methods
sample and data collection. A total of 25 brown hares (Lepus europaeus) of both sexes were included in this study. Metadata of all animals included is listed in Supplementary Table S1. The study including the experimental protocol was carried out and approved in accordance with all current laws of Austria and the guidelines and regulations of the institutional Ethics Committee and the institutional Good Scientific Practice Guidelines of the University of Veterinary Medicine, Vienna. Samples were all obtained during the hunting season in December 2015, in order to exclude any seasonal effects. In all hunting, respectively sampling areas, a systematized yearly drive-hunt using rifles according to the respective hunting law was performed. Sampling locations are spatially visualized in Fig. 1. The two Austrian populations are known to genetically diverse with geneflow occurring both locally and regionally 40,41 . The population on Pellworm has more restricted genetic diversity typical of an isolated island population (unpublished data). Sampled individuals were from well separated areas within each hunting ground to avoid the possibility of including individuals from the same family groups within each population sample.
As age could constitute an influential factor on the microbiota of hares, the age was determined by the weight of the dried eye lenses 42 . This method allows for the sensitive discrimination of subadult hares (i.e. young of the year or individuals not older than 9-10 months of age) that can otherwise not easily be distinguished from adult hares based on morphological characteristics. Hares that were younger than one year of age were classified as subadults. Digesta samples from the cranial flexure of the duodenum (intestinal samples) and faecal samples from the ampulla recti were collected in sterile sampling vials during a standardized necropsy by trained pathologists in the field within 20 minutes after death and immediately frozen at −20 °C, transported to the lab and frozen at −80 °C until use.
Heart fat (coronary fat deposit) as a parameter defining the nutritional status was assessed ranging from 1-6, such that 1 defined cachectic, 2 bad, 3 moderate, 4 good, 5 very good and 6 described obese. The parameter gut health was defined based on the following variables assessed by pathological examination of the gastrointestinal tract: 1: no macroscopic and/or no/very mild histopathological lesions (smooth mucosa, no signs of inflammation (0-+), no increase in inflammatory cells); 2: moderate dysbiosis signs of inflammation, ((+-++) inflammatory cells); 3: severe macroscopic and/or histopathological lesions (severe mucosal erythema, severe signs of inflammation, and (++-+++) inflammatory cells). The animals originated from three different geographical locations associated with different land use practices. As hares are highly selective feeders that prefer weeds/ grasses and various crop types while avoiding cereals 14 , we chose the sampling sites according to the different land use types: (a) Grassland-type (Pellworm): Around 80% of the whole area of the island, Pellworm, northern Germany, are used as agricultural land; around 70% of the agricultural land consists of grassland (pastureland and green fodder for domestic species such as cattle and sheep); the remaining 30% are used as cropland; (main crop consists of winter wheat (Triticum aestivum), rapeseed (Brassica napus), as well as to a lower percentage, barley (Hordeum vulgare), oat (Avena sativa) and Flax (Linum usitatissimum) 43,44 . Around 8% are used for maize (Zea mays) plantation for a local biogas facility 45 . (b) Cropland-type (Lower Austria-District Mistelbach): Of the 89% agricultural land of the total area of the district, around 79% are used as cropland, where predominantly various species of grain and maize (Zea mays) are cultivated followed by oil seed, sugar beet (Beta vulgaris), wine and fruit plantation as well as fallow land; to a lower percentage, potatoes (Solanum tuberosum) and protein rich plants are cultivated. Only around 0.4% are defined as grassland 46 . (c) No agricultural use (Military Airport Langenlebarn, 16S rRNA gene sequencing and sequence processing. For amplicon sequencing, we targeted the 16S rDNA hypervariable region V345, using an Illumina MiSeq sequencing platform with a 300-basepair paired-end read protocol. The universal primers 357F-HMP (5′-CCTACGGGAGGCAGCAG-3′) and 926R-HMP (5′-CCGTCAATTCMTTTRAGT-3′) were used to generate amplicons of ~570 base pairs 48,49 . Libraries were constructed by ligating sequencing adapters and indices onto purified PCR products using the Nextera XT Sample Preparation Kit (Illumina) and equimolar amounts of each of the libraries were pooled and sequenced. 16S rRNA gene PCRs, library preparation and sequencing, as well as demultiplexing and read stitching were performed by Microsynth (Microsynth AG, Balgach, Switzerland). process-and sequencing-controls. Prior to sampling, DNA extraction was attempted from blank tubes and a PCR with a long 16S rRNA gene primer panel was performed to verify that sampling tubes were not contaminated. A negative control was also included in the DNA extraction run. This negative control was also sequenced and should be a proxy for reagent contamination and cross-contamination between samples. Sequencing of the negative control revealed 5,000 sequences (average read count levels per sample was 77,000). Three non-template controls and one positive control (human stool sample) were included as internal controls in library preparation. These controls were identified as expected (3 negatives, one positive).The MiSeq Run was controlled by a PhiX control sample (https://emea.illumina.com/products/by-type/sequencing-kits/cluster-gen-sequencing-reagents/ phix-control-v3.html?langsel=/it). www.nature.com/scientificreports www.nature.com/scientificreports/ Microbial community analysis and predictive functional profiling. Analysis of the NGS data was conducted using the QIIME software package v1.9.1 50 . A total of 3,856,491 demultiplexed sequences were produced, of which 2,700,987 (70.04%) passed the quality filter with a phred value of >9 and a minimum read length of 520 base pairs. A total of 7,947 chimeric sequences were filtered out before operational taxonomic units (OTUs) were clustered using 97.00% identity (=0.03 distance) against the USEARCH 6.1 database 51 and de novo. Singletons were removed, and representative sequences were classified via the RDP classifier 2.2 52 . Alpha diversity was calculated by subsampling the OTU table to the read depth of the sample with the lowest amount of sequences (n(seqlowest) = 8,671), resulting in 10 repetitions at 1,000 sequences up to 8,000 sequences/sample. Beta diversity was calculated by subsampling to an even read depth of 6,000 reads/sample, of which weighted UniFrac distance matrices 53 and principal coordinate analysis were calculated 54 . Venn diagrams of OTUs were generated with the Venn diagram plotter (PNNL, Richland, WA) to display the number of OTUs shared by different populations. In addition, the core microbiota OTUs, defined as OTUs shared by over 90% of the samples, were identified. BLASTn analysis (NCBI GenBank database excluding uncultured/environmental sample sequences in the search set) of the top 50 most abundant OTUs and core microbiota OTUs were done using the highest identity matches of each unique OTUs representative sequence.
For qPCRs, DNA samples and negative controls were run in duplicate in a 20 µl reaction. The master mix contained 10 µl 2 × Brilliant III Ultra Fast SYBR Green qPCR Master Mix (Agilent, Vienna, Austria), 2 µl of 2.5 µM primers and 5 µl of water (nuclease-free). 1 µl DNA template including approximately 10 ng DNA (intestinal samples) and 1 ng DNA (faecal samples) was added. The following amplification protocol was used: Initial denaturation at 95 °C for 3 min, 40 cycles of 95 °C for 5 s and by 20 s at 61 °C and 63 °C for general bacteria and Eco1457F -Eco1652R qPCR, respectively. A melting curve (70 °C to 90 °C, with fluorescence measurements at 1 °C intervals), was done after each qPCR and specific amplicon peaks were received. qPCR results were analysed using the Stratagene MxPro software (QPCR Software, version 2.00). statistical analyses. Statistical analysis of the bacterial microbiota composition was implemented using the R statistical computing environment, version 4.3.3 59 . Alpha diversity indices were assessed by calculating the Chao1 index and the observed OTU richness. We used lme models (linear mixed-effect models, R-package, nlme 60 ) to analyse the effects of independent variables, i.e. population (geographical areas (Airport (A), Lower Austria (LA) and Pellworm (P)), sex (male vs. female), age class (sub adult vs. adult), and heart fat as a parameter defining the nutritional status of the animal (ranged 1-6, whereas 1 defined cachectic, 2 bad, 3 moderate, 4 good, 5 very good and 6 described obese) on diversity indices (dependent variable). Beside these fixed effects (geographical area, sex, age, heart fat and gut health), we also used the individual animal IDs (n = 25) as a random effect. All predictor variables in the full model were checked for multicollinearity by calculating the variance inflation factor. We used the Akaike Information criteria corrected for small sample size (AICc) to determine the best model with the most relevant fixed factors for each response variable. None of the models indicated primary collinearity issues (all values were below 2.6 61 ). The residuals of the models were assessed visually via histograms and Q-Q-plots as well as by calculation of the Shapiro-Wilk normality test. The dependent variable Chao1 calculated from the intestinal tract data was log transformed. The modelling procedure was carried out for both the faecal and intestinal datasets, although the effect of the faecal and intestinal samples on diversity indices itself was also calculated in an individual model run. We estimated the proportion of the variance of the fixed and random regressors using Pseudo-R-Squared (pseudo-r 2 -values 62 ). In the next step, contrast coefficients of the explanatory variables were calculated among the response variables diversity indices by using a pairwise comparison of least-squares means, corrected for multiple comparison using the false discovery rate (FDR). The contrast calculation was implemented in R using the package lsmeans 63 . The level of significance was defined at P ≤ 0.050. The normality distribution of the abundance of phyla (n = 23), and most abundant OTUs (n = 50) were checked with the tests for multivariate data with the function manova (R Package stats) and the Cullen and Frey graph with the R package fitdistrplus and logspline 64,65 . Due to the non-normal distribution of the data, model residuals, and different identified optimal theoretical distributions of individual phyla and OTUs, the beta diversity of the microbiota of phyla and OTUs was assessed by applying a permutational multivariate analysis of variance (PERMANOVA) with the ADONIS function and 5000 permutations in R (R Package vegan 66 ). To analyse whether independent variables (i.e. population, sex, age, heart fat and gut health) had also a significant effect on the composition of phyla and OTUs abundances, we calculated Bray-Curtis-dissimilarity matrix with the R function vegdist of the vegan package 66 . The distance matrix was applied as response variable. The modelling approach was carried out for faecal and intestinal samples, although the effect of the two on the composition of www.nature.com/scientificreports www.nature.com/scientificreports/ phyla and OTUs itself was also calculated in an individual ADONIS function. The multivariate homogeneity of group dispersions was performed with the betadisper function, a multivariate analogue of Levene's test, followed by a permutation-based test of multivariate homogeneity of group dispersions with pairwise comparisons of group mean dispersions. The Kruskal-Wallis test was applied to identify potentially statistically significant differences between populations concerning each individual phylum, OTU, and independent variables, followed by a pairwise test for multiple comparison of mean rank sum (Dunn's test), adjusted with Benjamini-Hochberg method (R Package PMCMR, FSA 67,68 ).
Relative abundance plots of microbial communities were calculated with the transform_sample_counts function within the phyloseq package. The corresponding visualisation was conducted with the package ggplot2 in R 69,70 . The tool 'Phylogenetic Investigation of Communities by Reconstruction of Unobserved States' (PICRUSt) was applied to predict the metagenome functional content from amplicon sequencing 71 . Statistical analysis of metagenome predictions was carried out with STAMP 72 by using an ANOVA with Benjamini-Hochberg FDR correction.
Additionally, the Kruskal Wallis test was applied to compare the different populations in terms of qPCR abundance of Eco1457F-Eco1652R and all bacteria in faecal samples and intestinal tract samples. A post-hoc test using Mann-Whitney tests with Bonferroni correction for multiple comparisons was used and the corresponding Z Values were computed to determine the effect size (r).

Results
the microbiota of intestinal and faecal hare samples. In total, 2,700,987 sequences (70.04%) passed our quality filter and were processed together for all downstream analysis. Reads generated for intestinal and faecal samples were distributed homogenously, with 53% reads belonging to intestinal and 47% belonging to faecal samples. At both microbiota sources analysed (intestinal and faecal samples) Firmicutes and Bacteroidetes dominated the microbial community composition (45.51% and 19.30% relative abundance over all samples respectively), followed by Spirochaetes and Proteobacteria (7.76% and 5.73% relative abundance over all samples) (  Table 1). The core microbiota consisted of 66 shared OTUs, which all belonged to the 500 most abundant OTUs across all samples ( Table 2).
Observed diversity richness OTU counts and diversity estimator Chao 1 counts are listed in Supplementary  Table S2. Over all populations, the average observed diversity richness per sample was 711 OTUs in faecal and 510 OTUs in intestinal samples. Chao 1 estimated the diversity on average to contain 1,686 OTUs in faecal and 932 OTUs in intestinal samples. Additionally, a statistically significant effect between intestinal and faecal abundance for the diversity indices were observed (Chao 1, P < 0.001 and observed OTUs P < 0.001). The sampling site could explain approximately 22-32% of the variance of the diversity data, but most of the variance was explained at the individual level (OTU richness: conditional r 2 = 0.57: Chao: conditional r 2 = 0.56).  (Table 3). OTU 296045 and OTU 1108377 (best BLAST hit Bacteroides sartorii and Muribaculum intestinale) had a significantly higher abundance in the Pellworm population compared with the Airport population (Table 3). NR1 (best BLAST hit Murimonas intestine) was higher in the Airport www.nature.com/scientificreports www.nature.com/scientificreports/ population compared with the Pellworm population (Table 3). Heart fat, sex and gut health parameter had no significant effect on intestinal or faecal-associated OTUs (P > 0.050). In this context, no significant effect of population, age, heart fat, sex and gut health on the diversity abundance in intestinal and faecal samples was observed. Our models of alpha diversity indices (Chao1 index and the observed OTU richness) could explain on average 89% of the variance of the data (conditional r 2 = 0.88 and r² = 0.89 for intestinal and faecal samples, respectively). However, much of the variation in the data was due to the difference between individuals, since data variance explained by the fixed factors alone was considerably lower (mean marginal r 2 = 0.08 and r² = 0.14 for intestinal and faecal samples, respectively). Abundant OTUs summed up at family level are shown in Fig. 2B.
The Venn diagram ( Supplementary Fig. S1) displays an unequal amount of OTUs per population and a high number of unique OTUs in each population. In faecal samples, 7.69% of all OTUs, and in intestinal samples 8.97% of all OTUs were shared between the three populations examined. The phylogenetic distance measurement done with weighted UniFrac analysis did not reveal a clear clustering of populations, however the Airport population samples were more similar to each other compared with the other populations for both intestinal and faecal samples. The highest variability within a population was found for the LA population (Fig. 3). The PICRUSt analysis indicated one significant difference on KEGG (Kyoto Encyclopedia of Genes and Genomes) Level 3 in the population comparison: The linoleic acid metabolism was significantly higher in Airport and Pellworm populations compared with LA population (Effect size: 0.36, P < 0.005; Supplementary Table S3).  www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
This is to our knowledge the first study characterising the intestinal and faecal microbiota of the European brown hare based on 16S rRNA gene sequencing.
To date, studies describing the gastrointestinal microbiota are limited to one member of the family Leporoidae, the rabbit 73,74 . As sympatric relatives, often sharing the same habitat as well as feed preferences, the species do share similarities of their gastrointestinal physiology 75 . Leporidae are hindgut fermenters and practitioners of caecotrophy. Nevertheless, hares and rabbits show both morphological differences of their gastrointestinal tract (GIT) and different digestive strategies in terms of nutrient extraction 14,75,76 . Those differences in digestive function and strategies between the two species suggest that microbial communities in the GIT may also differ and potentially show divergent functionality or metabolic pathways. Indeed, we were able to identify that the core microbiota of the hare GIT diverges from that previously described in rabbits. Our study reveals that the gut microbiota of hares includes Clostridium, Bacteroides, and Ruminococcus species, but not Streptococcus and Enterobacter species as previously described to be abundant in rabbit caecum [77][78][79][80][81][82] . Despite these differences we also detected similarities between the gut microbiota of the two species, e.g. Eubacterium and Ruminococcus related OTUs, both found among the 50 most abundant and core microbiota OTUs in this study, have been described in the rabbit GIT as well 83 . Although the highly abundant phylotypes detected in this study were in accordance with phylotypes recently described in rex rabbits 74 , hare intestinal and faecal samples were highly enriched for Spirochaetes compared with rabbits. The majority of the phylum Spirochaetes consisted of the family Sphaerochaetaceae mainly comprising of OTUs affiliated to Sphaerochaeta pleomorpha or Sphaerochaeta globosa. These OTUs showed only 88% sequence similarity to reference 16S rRNA gene sequences but 93% to 96% sequence similarity to 16S rRNA gene sequences of uncultured Spirochete clones found in rumen of ruminants or the oral cavity of canines 84,85 . Interestingly, BLASTn database sequences include a large number of uncultured Sphaerochaeta species residing in a broad range of mammalian GIT or oral cavities. Unlike most other Spirochaetes, which usually have a helical morphology and flagella, Sphaerochaeta are spherical and not motile 86 . They are capable of heterofermentative growth on carbohydrates such as pentose and hexose monosaccharides, disaccharides and soluble starch 86 . In contrast to their closest relatives Spirochaeta, Treponema, and Borrelia, they are not thought to be pathogenic 87 .
Similar to the Sphaerochaeta associated OTUs, other high-abundance OTUs observed in this study had similarity scores to reference 16S rRNA genes below 90%, which is in accordance with a study on bacterial species in the rabbit caecum where the majority of sequences shared less than 97% identity to BLASTn database sequences 83 . Hence our data demonstrate that gastrointestinal tracts of hares harbour multiple bacterial species that that are yet to be described.
Studies in wildlife species are often restricted in sample availability due to limitations of access to free-ranging animals and restrictions in the application of invasive sampling methods. Non-invasive samples such as faeces are in many cases the only available sample rescource. In order to evalutae the representation of the gastrointestinal microbiota composition in faecal samples and therefore its potential use in future studies in hares, we tested the similarity between the faecal and intestinal microbial community composition. www.nature.com/scientificreports www.nature.com/scientificreports/ Firmicutes and Bacteroidetes dominated the microbial community composition in both sample types, followed by Proteobacteria and Spirochaetes. Intestinal samples showed a higher number of phyla with a significantly different abundance pattern compared with faecal samples. Similarly, diversity indices and the abundance of OTUs differed significantly between sampling sites, indicating clear microbiota shifts along the hare's gastrointestinal tract. Faecal samples were not representative for the description of the intestinal microbiota. This is in accordance with recent studies that identified a distinct microbial environment between luminal and mucosal sites and in respective compartments of the gastrointestinal tract of monogastric animals 88,89 . Therefore, future studies should consider using GIT samples representative for testing their specific hypothesis.
Land use practices associated with poorly diversified and restricted food availability is known to effect gut microbiota and individual fitness 11 . Highly specialized feeders, such as the brown hare, might be even more severely affected by a decrease in plant biodiversity 14 . Our three chosen sampling locations differ tremendously in land use and habitat type parameters. The agricultural land use in Pellworm consists predominantly of grassland and to a much lower percentage cropland; whereas the situation in LA is the opposite, with the majority of land used as cropland. The control site Airport consists of fallow land without any agricultural use. As highly selective feeders, hares prefer weeds and grasses, and specifically select for certain plant taxa 14 . Furthermore, they avoid crude fibre and select for high fat and energy content in their diet to provide essential polyunsaturated fatty acids (PUFA) such as linoleic acid and alpha-linolenic acid needed to reproduce and survive. A study by Popescu and colleagues found that European hares are able to selectively absorb PUFA in the gastrointestinal tract and excrete faeces that are highly depleted in PUFA and enriched in saturated fatty acids (FA) 76   metabolism was significantly elevated in the Airport and Pellworm populations compared with the LA population (Supplementary Table S3), probably indicating the variation in food and nutrient variability between the populations in our study.  www.nature.com/scientificreports www.nature.com/scientificreports/ The overall phyla composition was relatively stable in hares (only intestinal samples varied between populations as a trend), but hare location had a significant impact on the OTU composition of faecal and intestinal samples between the different populations. In intestinal samples a number of specific taxa differed in abundance between the different populations. For example, there was a trend for a difference in the abundance of Bacteroidetes between LA and Airport and a significant difference in the abundance of OTU 296045 (best BLAST hit: Bacteroides sartorii), which was enriched in intestinal samples of the Pellworm population relative to the Airport population. Bacteroides species have been detected at high abundances in rabbit caecum in the past and are associated with healthy growth and high-weight gain 74,90 . Overall, the Pellworm population had the lowest OTU diversity. Population genetic data also suggest that Pellworm has lower diversity than the Austrian populations and its restricted geneflow is likely also associated with the low OTU diversity. Further research on associations between OTU diversity on gut health, respectively disease susceptibility and on physiological parameters is warranted. Hares of the Airport population showed lower abundances of the yet uncultured bacterial family S24-7 compared with the Pellworm and LA population. This family is being recognized as a predominant member of the intestinal tract of homeothermic animals, but data on their physiology, metabolic capacity, and interactions with the host are limited 91 . Plant glycan (hemicellulose and pectin) degradation has been described based on genomic characterization of the S24-7 family 91 . Our results suggest that both grassland and cropland agricultural areas are associated with increased abundances of S24-7 affiliated bacteria in the intestinal tracts of hares.
Population had a significant effect on both the composition of intestinal and faecal OTUs. Population was also the explanatory factor for the abundance of potential pathogenic bacteria of the family Enterobacteriaceae.
Many members of the family Enterobacteriaceae, which showed the highest abundance in the LA population and lowest abundance in hares of the Pellworm area, are connected to gastrointestinal disease [92][93][94] . The low abundance of Enterobacteriaceae and the high abundance of OTU 296045 (i.e, likely Bacteroides) in the GIT of Pellworm hares might indicate a healthier microbial composition in the intestinal tract. Host traits such as age, sex, body condition (heart fat) and gut health were not related to the abundance and diversity of the bacterial composition in intestinal and faecal samples of brown hares in this study.
Nevertheless, further investigations that look into associations of agricultural use, plant diversity and occurrence of dysbiosis are desirable. Considering the low number of sampled animals (n = 25) in the presented study, significant impact of independent variables should be interpreted with caution. A larger sample size in future studies would be desirable and would allow for more robust conclusions than those derived from the present study.

Conclusion
The current knowledge of gut microbiota composition and influencing factors in wildlife species under natural conditions is extremely limited. Our aim was to provide a broad foundation to generate specific hypothesis at the intersection of gut health, land use change and population viability for the EBH and other species impacted by rapid habitat modification. This combined information is lacking to date, despite being essential for understanding gut microbial variation within wildlife species to differentiate natural variability from that related to external stressors leading to functional dysbiosis affecting host health. Anthropogenic habitat modifications together with species-specific susceptibility to these modifications might also have an impact on wildlife gut microbiota.
This study generates fundamental baseline data on the diversity and composition of gut microbiota in the brown hare that will complement previous research and create a starting point for future research to understand the physiological and pathophysiological dynamics of gut microbiota as well as the causes of alterations to the intestinal microbial communities in this species.

Data Availability
Nucleotide sequence accession numbers: Sequencing data are available in the European Nucleotide Archive Database under the accession number PRJEB15166.