Sympatric Lepus spp. in the central Italian Alps host significantly different gut microbiotas

Sympatric Lepus spp. in the central Italian Alps host significantly different gut microbiotas | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Article Sympatric Lepus spp. in the central Italian Alps host significantly different gut microbiotas Lara Marinangeli, Barbara Crestanello, Nadine Praeg, Theresa Rzehak, and 5 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-7375012/v1 This work is licensed under a CC BY 4.0 License Status: Published Journal Publication published 04 Apr, 2026 Read the published version in Scientific Reports → Version 1 posted 5 You are reading this latest preprint version Abstract The mountain hare (Lepus timidus) is an arctic-alpine species with relictual populations in the Italian Alps, typically occurring at elevations above 2000 m a.s.l. This species is threatened by habitat loss and fragmentation, and declining snow cover due to climate warming. Moreover, as treelines shift upward, the European brown hare (L. europaeus) is expanding its distribution into areas previously dominated by the mountain hare, potentially leading to resource competition, and loss of local adaptation through hybridization and inter-specific gene flow. In particular, the consequences of sympatry on diversity and composition of prokaryote and fungal communities of the gut microbiota, which are critical to individual health, are currently unknown. Here, we compared the gut microbiota of these two hare species in an area of overlap in the central Alps by analysing fresh faecal pellets collected from Val Mazia/Matschertal, Italy along an elevational gradient (1000 to 2500 m a.s.l.). For the first time, we describe the prokaryote diversity and composition of L. timidus, and the fungal gut communities (mycobiota) of both Lepus species. Species identity was confirmed for 95 samples via mtDNA barcoding, while gut microbiota richness and composition were investigated using amplicon sequencing, targeting the V3-V4 region of the prokaryote 16S rRNA gene and fungal ITS2 regions. Distinct prokaryote and fungal communities were observed for each species, even in sympatry, indicating differences in their functional diversity. Interestingly, for both Lepus species, elevation influenced fungal but not prokaryote diversity. Therefore, sympatry appears to have had minimal impact on gut microbiota composition of either species thus far. Given the expected upward range shift of L. europaeus under climate warming and its continued restocking for hunting, our findings provide an important baseline for assessing the health and adaptability of L. timidus as well as the effectiveness of conservation efforts aimed at protecting L. timidus. However, expanding this research to other areas of sympatry will be essential to understand if gut microbial composition is indicative of L. timidus conservation status across its range. Lepus europaeus Lepus timidus metataxonomy mtDNA 16S rRNA gene ITS2 Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Introduction One of the most dramatic and rapid impacts of climate warming in the Alps is the shifting of vegetational zones to higher elevations [ 1 , 2 , 3 ], accompanied by upward distributional changes of animal species as they track their preferred habitats [ 4 , 5 , 6 ]. An emblematic example, estimated from 30 years of hunting bag data [ 7 ], is the mountain hare ( Lepus timidus ), a species generally associated with boreal zones above the treeline [ 8 , 9 ]. However, results in [ 7 ] also indicated that the European brown hare ( Lepus europaeus ), a steppe-adapted species [ 10 ], is also shifting upward even more rapidly, which may lead to a reduction in mountain hare habitat. In addition, sympatry of the two hare species could result in their hybridization, as already reported by several authors for northern European as well as Alpine populations (e.g. Fennoscandia: [ 11 ]; Swiss Alps: [ 12 ]; western Italian Alps: [ 13 ] and references therein). These risk factors could threaten the long-term persistence of the mountain hare and its unique genetic pool in Alpine habitats. Interspecific hybridization and range shifts (associated with dietary changes) could also lead to the loss of local adaptation, for example by altering the gut flora, which in hybrids can be intermediate between the two parental lineages (e.g. equids: [ 14 ]; suids: [ 15 ]). It is now well-recognized that an intact gut microbiota (including prokaryotes, fungi, and viruses) is essential for immune and metabolic function, as well as for development of the neuroendocrine system [ 16 , 17 , 18 ], impacting host health and resilience. It has been suggested that the presence of specific microorganisms in the gut may be important for host fitness and should be considered for the conservation and management of wild animal species [ 19 , 20 ]. While host phylogeny has a strong influence on gut flora for many mammals, with various authors reporting species-specific compositions [ 21 , 22 ], extrinsic factors such as habitat preference and diet are also expected to have an impact on gut flora alpha and beta diversities [ 23 , 24 ]. Since the European brown hare is rapidly colonizing higher elevations in the Alps due to climate warming, increasing sympatry with the mountain hare is inevitable. Thus, defining gut microbiota of the two hare species is useful for understanding the potential consequences of their interaction. Here, for the first time, we describe and compare the gut microbiota composition of the two Lepus species in an area of sympatry in Val Mazia/Matschertal, Province of Bolzano, Italy using MiSeq amplicon sequencing of the 16S rRNA and ITS2 genes. This research aims to understand how host species and elevation influences the microbiota of these lagomorphs. Although the prokaryote gut microbiota of the European brown hare has previously been reported by Stalder et al. [ 25 ], no similar data exist for that of the mountain hare. In addition, the mycobiota of both species has not been studied to date, despite its potential relevance for host health [ 26 ]. Importantly, no studies have yet examined the gut microbiota of sympatric populations. Materials and methods Study area and sample collection Single fresh pellets of Lepus spp. were collected at least 20 m apart, from 12 sampling areas, three from each of four elevations (1000, 1500, 2000 and 2500 m a.s.l.) in July 2019 and 2020 in the Vinschgau Valley LTSER, South Tyrol, Italy (site code LTER_EU_IT_097, 46.6928°, 1.6157°; https://deims.org/11696de6-0ab9-4c94-a06b-7ce40f56c964 ; Fig. 1 ) as part of the EUREGIO project MICROVALU. Freshly deposited faecal pellets were identified by their shiny (mucous smeared) external surface. Fieldworkers wore N95 masks and used sterile gloves and tweezers to collect the pellets, which were placed in certified sterile DNA/DNAse-free 15 ml tubes (Sarstedt, Germany) and stored at -20°C for up to 24h before being transferred to the Animal, Environmental and Antique DNA Platform at the Fondazione E. Mach (FEM), where they were archived at -80°C until further processing. Sample details are provided in the Supplementary Tab. S1. Meteorological data were collected on site using Campbell Scientific CR1000/CR1000X loggers and provided by Zandonai et al. [ 27 ], including daily average temperature in the three weeks prior to sampling and sum of precipitation of the entire vegetation period (May 1st to Sept 30th ). Temperature was recorded with three loggers per elevation, while precipitation was recorded by one meteorological station at each elevation. Species richness was estimated using a square plot with a grid of 100 cm 2 (10x10 cm) units to calculate the number of plant species in each elevational zone, while plant richness was obtained from Hilpold et al. [ 28 ]. Botanical data of the 2500 m a.s.l. sites was then added to the previously published elevation surveys. DNA extraction Under a sterile Class II biological safety cabinet (BSL2, Telstar, Spain) and using sterile DNA/DNase-free equipment and consumables, 50 mg of faecal material were taken from the centre of each pellet to avoid the outer layer that had potentially been in contact with soil. Total DNA was extracted from each 50 mg subsample of a single faecal pellet with the NucleoSpin Soil mini kit (Macherey-Nagel, Germany) using the lysis buffer SL1 in combination with 50 µl of Buffer SX, as suggested by Praeg, Pauli, and Illmer [ 29 ]. Each sample was co-extracted with a 0.15 dose of ZymoBIOMICS™ Spike-in Control I (High Microbial Load; EuroClone, Milan, Italy), before the initial lysis step following Galla et al. [ 30 ], unless otherwise stated (Tab. S1). This mock community was added as an in situ positive control to validate the accuracy and reliability of sequencing and bioinformatics pipeline. Negative DNA extraction controls were included up to the sequencing step to monitor environmental and reagent contaminations. The purity and quantity of all DNA extracts were assessed by checking the UV/VIS spectra of each extract with a Spark multimode microplate reader (Tecan, Switzerland). Following quantification, all DNA extracts were diluted to a final concentration of 3 ng/µl using nuclease-free water and subsequently used for species identification and amplicon sequencing. Host species identification To determine the origin of each pellet ( L. europaeus or L. timidus ), a ~ 500 bp-fragment of the mitochondrial DNA D-loop region was amplified using GoTaq DNA polymerase by Promega with primers (H) GTTGCTGGTTTCACGGAGGTAG and (L) TCCTACCATCAGCACCCAAAGC’ previously designed by Wilkinson & Chapman [ 31 ]. Cycling conditions consisted of an initial denaturation step at 94°C for 3 min followed by 40 cycles of 94°C for 30 s, 61°C for 30 s and 72°C for 45 s, with a final extension step at 72°C for 5 min. Negative PCR controls (PCR-grade water instead of DNA template) were included to check for contamination. Amplified DNA was purified using ExoProStarTM 1-Step (Cytiva, USA) and sequenced at the FEM Sequencing and Genotyping Platform with both primers H and L using the ABI Prism BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems). Sequencing reaction products were further purified with the Dynabeads Sequencing Clean-Up Kit (ThermoFisher Scientific) and subsequently analyzed on an ABI 3130XL DNA sequencer (Applied Biosystems, Foster City, CA, USA), 96 capillaries (50 cm) and POP07 Performance Optimized Polymer (Applied Biosystems). The sequences were processed with Sequencher software (Gene Codes Corporation-USA) to delete any automatic base assignment errors. BLASTn was used to identify sequences with ≥ 99% identity as defined using the Standard Nucleotide method. Moreover, D-loop sequences were trimmed to a 334 bp region shared by all sequences as well as 147 D-loop sequences deposited in NCBI by Melo-Ferreira et al. [ 32 ], Pecchioli et al. (unpublished PhD thesis) [ 33 ] and Zachos et al. [ 34 ] and aligned with MEGA v. 10.0.5 (option Muscle; [ 35 ]). A Neighbor-Joining phylogenetic tree was constructed in MEGA, using the Phylogeny option, based on pairwise genetic distance. Characterization of V3-V4 16S rRNA gene and ITS2 faecal communities The amplification of the V3-V4 regions of the prokaryote 16S rRNA gene was performed using the KAPA HiFi HS ReadyMix (Roche) in a 25 µl reaction volume containing 1X KAPA HiFi HS ReadyMix Buffer, 0.3 µM each of primers 341F_ILL (CCTACGGGNGGCWGCAG) and 805R-2_ILL (GACTACNVGGGTWTCTAATCC) modified with Illumina overhang adapters [ 36 , 37 ] ( https://support.illumina.com/documents/documentation/chemistry_documentation/16s ) and 9 ng of DNA. The thermal profile for 16S V3-V4 amplification reactions was: 3 min at 95°C, followed by 31 cycles of 30 sec at 95°C, 30 sec at 55°C, 30 sec at 72°C, and a single final extension step of 5 minutes at 72°C. The amplification of the fungal rRNA Internal Transcribed Spacer ITS2 was performed in a 25µl reaction volume containing 1X FastStart High Fidelity Reaction Buffer (Roche Applied Science), 0.4 µM each of primers gITS7 (Illumina_forward_overhang-GTGARTCATCGARTCTTTG) and ITS4 (Illumina_reverse_overhang-TCCTCCGCTTATTGATATGC) [ 38 , 39 ], 200 µM each dNTPs, 9 ng of DNA, and 1.5 U of FastStart High Fidelity Enzyme Blend (Roche Applied Science, Germany). The thermal profile for ITS2 amplification reactions was 3 min at 95°C, followed by 35 cycles of 45 sec at 95°C, 45 sec at 50°C, 45 sec at 72°C, and a single final extension step of 5 minutes at 72°C. All amplifications were performed on a Veriti 96-Well Fast Thermal Cycler (Applied Biosystems, USA). Non-template controls were included in each amplification reaction. Amplicons were visualized by high-resolution capillary electrophoresis using the QIAxcel Advanced System (Qiagen, Hilden, Germany). Quantification of individual amplicon libraries, normalization at equimolar concentrations, pooling of indexed libraries and high throughput sequencing by Illumina technology were performed at the FEM Sequencing and Genotyping Platform. Libraries were sequenced on an Illumina MiSeq platform using Standard Flow Cells (PE300), targeting a minimum sequencing depth of 30,000 reads per amplicon. Amplicon sequencing data analysis Bioinformatic analyses were performed in R v. 4.3.1. DADA2 v. 3.14 [ 40 ] was used for sequence quality inspection, primer trimming, chimera removal and ASVs taxonomy assignment following the standard operating procedure. Taxonomic classification of 16S rRNA gene and ITS2 ASVs was performed using Silva v. 138.1 [ 41 ] and UNITE v. 9.0 [ 42 ] reference databases, respectively. Multiple sequence alignment of 16S rRNA gene ASVs was done with msa v. 1.32.0 [ 43 ] using ClustalW and default parameters. With the Phangorn package v. 2.11.1 [ 44 ], a phylogenetic tree was constructed with 16S rRNA gene ASVs, which was then used for further analysis. Sequencing results of negative controls from both DNA extraction (N = 12) and PCR amplification (N = 13) were used to identify potential contaminants using decontam v. 1.20.0 [ 45 ] with the prevalence method. All mock community-related ASVs were removed from the dataset before performing subsequent analyses. Statistical analyses of 16S rRNA gene and ITS2 ASVs were performed using phyloseq v. 1.50.0 [ 46 ] and microeco v. 4.3.1 [ 47 ]. Images were generated using ggplot2 v. 3.3.5 [ 48 ], using RColorBrewer v. 1.1-2 [ 49 ] and viridis v. 0.6.5 [ 50 ]. The impact of species and elevation on taxonomic composition (alpha and beta diversity) of faecal microbiota were investigated by focusing on elevations where the two species live in sympatry (Tab. S1). L. timidus and L. europaeus datasets were then analysed separately to assess intraspecific variation in faecal microbiota across all elevations where they were sampled. Taxonomic profiling of 16S rRNA and ITS2 regions Taxonomic classification of 16S rRNA gene and ITS2 ASVs was estimated on non-rarefied data and expressed as relative abundances. Visualizations were generated using microeco, focusing on the 10 most abundant phyla and 16 most abundant classes and families. The relative abundance of prokaryote and fungal taxa was compared between the two species with LEfSE [ 51 ]. Alpha and beta diversity estimation To avoid bias in richness estimates due to different sequencing depths, alpha diversity indices Chao1, inverse Simpson, Shannon and Faith’s Phylogenetic Diversity (PD; this latter index for prokaryotes only) were estimated using 16S rRNA gene and ITS2 datasets rarefied to 10,564 and 10,804 ASVs per sample, respectively. Alpha diversity estimates were compared between Lepus species and elevations using Wilcoxon rank-sum tests, with microeco. Venn diagrams representing the number of ASVs shared between L. timidus and L. europaeus were generated using the rarefied data, microeco. Beta diversity estimates were computed on non-rarefied ASV abundances normalized as relative abundances using Bray-Curtis, unweighted Unifrac (for 16S rRNA gene only) and Jaccard (for ITS2 only) indices. Ordination based on distance/dissimilarity matrices was done using non-metric multidimensional scaling (NMDS). Multivariate homogeneity of group dispersions, was tested using betadisper with Vegan v. 2.6.4 [ 52 ]. To compare beta-diversity estimates between groups, permutational multivariate ANOVA (PERMANOVA), and pairwise-PERMANOVA were performed with 999 permutations using the function adonis2 implemented in microeco. The same package was used for the analysis of similarities (ANOSIM). Pairwise-PERMANOVA p-values were adjusted using the false discovery rate (FDR) method to control for multiple testing. Impact of environmental variables on beta diversity Significant differences in daily temperature, pre-sampling total precipitation and plant species richness between elevations were inspected using a Wilcox rank-sum test (data not shown). Autocorrelation among environmental variables was assessed using Spearman's Rank-Order Correlation using microeco. For both 16S rRNA gene and ITS2 datasets, Mantel tests employing Spearman rank order and Pearson correlation coefficients were used to examine associations between environmental parameters and beta diversity estimates. Spearman rank-order correlations were calculated to evaluate relationships between environmental parameters and the abundance of each 16S rRNA and ITS2 ASVs, also using microeco. Characterization of prokaryotic functional diversity Prokaryote functional diversity was inspected using the PICRUSt2 pipeline v. 2.5.3 [ 53 ] implemented on Galaxy ( https://usegalaxy.eu/ ) with default settings. Estimated MetaCyc pathway abundances were normalized to relative abundances. Alpha diversity (Shannon index) and beta diversity (Bray-Curtis dissimilarity) were estimated from normalized MetaCyc pathway abundances using microeco. Differences in MetaCyc pathway abundances between Lepus species and elevations was tested using Wilcoxon rank-sum text with FDR p-value adjustment. Differential pathways were defined by a Log 2 FC cutoff of ≥ 1 and an adjusted p-value of ≤ 0.01. Plots were generated with ggplot2. A heatmap clustering samples and pathways based on MetaCyc pathways abundances was generated using ClustVis [ 54 ] with default parameters and further formatted using GIMP v. 2.10.18 [ 55 ] ) . Results Lepus species identification Of 108 total fresh faecal pellets, 95 provided reliable mtDNA D-loop results, identifying 72 L. europaeus and 23 L. timidus samples. As expected, L. timidus was restricted to higher elevations, with nine and 14 samples found at 2000 m and 2500 m a.s.l., respectively. In contrast, L. europaeus ranged across the study area, with 27, 26, 16 and 3 samples collected at 1000, 1500, 2000, and 2500 m a.s.l., respectively (Fig. 1 ). Taxonomic composition of L. timidus and L. europaeus faecal microbiota The profiling of L. europaeus and L. timidus faecal microbiota identified 6145 ASVs for the prokaryote 16S rRNA gene and 2219 ASVs for the fungal ITS2 region. Of these, 22 and 16 ASVs were identified by decontam as potential contaminants and removed from the 16S rRNA gene and ITS2 datasets, respectively. Additionally, 24 ASVs were classified as belonging to the two taxa composing the mock community: Allobacillus spp. (7 ASVs) and Imtechella spp. (17 ASVs); all of these ASVs were also removed from the dataset before statistical analyses. The overall median relative abundance of these mock community-related ASVs was 3.2% (ranging from 0.2% to 15%) and was comparable between Lepus species and elevations (Kruskal-Wallis, p-value: 0.3487). Following rarefaction, 16S rRNA gene and ITS2 datasets accounted for 4877 and 1514 ASVs, respectively (Fig. S1 a-b). Rarefaction resulted in the loss of seven 16S rRNA gene (5 L. europaeus and 2 L. timidus ) and five ITS2 libraries (three L. europaeus and two L. timidus ). Both species were characterized by private ASVs accounting for approximately 82% of 16S rRNA and 12% ITS2 sequence reads. In contrast, only about 6% of prokaryote and 9% of fungal ASVs were shared by faecal samples of both species. Private L. europaeus prokaryote ASVs prevailed in the gut community which consisted of about 70% of sequenced reads. In contrast, dominant fungal ASVs were shared, accounting for approximately 84% of generated sequence reads (Fig. S1 a-b). Restricting the analysis to elevations where the two species were sympatric produced results consistent with those above, with approximately 45% of prokaryote ASVs found only in L. europaeus samples, 35% of prokaryote ASVs uniquely found in L. timidus samples (Fig. S1 c) and only a small fraction of prokaryote ASVs found in both species with relatively high abundance (e.g. 10 ASVs detected in both species and elevations accounted for 10% of sequence reads). Regarding fungi, most ASVs were shared by both species at all elevations (Fig. S1 d). The taxonomic classification of prokaryote phyla and fungal classes found in L. europaeus and L. timidus faecal pellets is shown in Fig. 2 . At all elevations, the faecal microbiota of L. europaeus was dominated by the prokaryote phyla Firmicutes (median: 63.3%, range: 0.28% − 80.8%) and Bacteroidota (median: 16.4%, range: undetected − 33.0%) followed, but with much lower percentages, by Actinobacteriota (median: 2.0%, range: undetected − 70.3%), Spirochaetota (median: 1.2%, range: undetected − 17.1%) and Patescibacteria (median: 1.0%, range: undetected − 30.5%). In addition, about 15% of the 72 samples collected across all four elevations were characterised by a notably high abundance of Proteobacteria (median: 0.7%, range: 0.1% − 88.0%) (Fig. 2 a). In contrast, L. timidus faecal microbiota was dominated by Proteobacteria (median: 62.9%, range: 4.6% − 78.9%), Bacteroiodota (median: 22.3%, range: undetected − 29.2%) and, to a lesser extent, Actinobacteriota (median: 5.7%, range: 0.2% − 25.8%) and Acidobacteriota (median: 1.2%, range: undetected − 4.2%) (Fig. 2 a). Firmicutes was detected in L. timidus faecal samples as well, but in lower percentages (median: 0.3%, range: 0.1% − 26.2%). The faecal mycobiota of both species was characterized by a high relative abundance of Dothideomycetes ( L. europaeus median: 51.8%, range: 6.3–87.9%; L. timidus median: 43.6%, range: 2.8–75.1%), Leotiomicetes ( L. europaeus median: 15.5%, range: 1.0% − 74.9%; L. timidus median: 22.1%, range: 0.4% − 81.3%) and Sordariomycetes ( L. europaeus median: 6.0%, range: undetected − 49.1%; L. timidus median: 0.9%, range: 0.1–20.3%). L. europaeus faecal samples also had non negligible occurrences of Pezizomycetes (median: 1.2%, range: undetected − 66.6%) and Microbotryomycetes (median: 0.1%, range: undetected − 40.9%) (Fig. 2 b), whereas L. timidus faecal samples were notable for the abundance of Tremellomycetes (median: 14.9%, range: 0.1% − 53.4%), this latter class being particularly abundant in samples collected at 2500 m a.s.l. (Fig. 2 b). The taxonomic composition of faecal microbiota was significantly different between the two study species where they were found to be sympatric (2000 and 2500 m a.s.l.). As shown in Fig. 3 and Tab. S2, L. europaeus faecal microbiota was significantly enriched with the following prokaryote phyla: Firmicutes (LEfSE; LDA = 5.47, adjusted p-value < 0.0001), Spirochaetota (LEfSE; LDA = 4.11, adjusted p-value < 0.0001), Verrucomicrobiota (LEfSE; LDA = 3.57, adjusted p-value < 0.0001). Instead, L. timidus faecal microbiota had significantly higher abundances of Proteobacteria (LEfSE; LDA = 5.42, adjusted p-value < 0.0001), Actinobacteriota (LEfSE; LDA = 4.05, adjusted p-value < 0.0001), Acidobacteriota (LEfSE; LDA = 3.93, adjusted p-value < 0.0001) and Myxococcota (LEfSE; LDA = 3.38, adjusted p-value < 0.001), among others (Fig. 3 a, Fig. S2 , Tab. S2). We also found differences in relative abundance of fungal classes between the two host species at 2000 and 2500 m a.s.l., mainly related to taxa belonging to Ustilaginomycetes (phylum Basidiomycota; LEfSE; LDA = 4.00, adjusted p-value = 2.57E-3) and Pezizomycetes (phylum Ascomycota; LEfSE; LDA = 3.86, adjusted p-value = 4.13E-3), that were all more abundant in L. europaeus (Fig. 3 , Tab. S2). Instead, the two Basidiomycota fungal classes Tremellomycetes (LEfSE; LDA = 4.93, adjusted p-value = 9.34E-4) and Cystobasidiomycetes (LEfSE; LDA = 3.67, adjusted p-value = 5.57E-3) had a higher relative abundance in L. timidus (Fig. 3 b, Fig. S2 , Tab. S2). Interspecific variation in diversity of L. timidus and L. europaeus faecal microbiota Wilcoxon rank-sum tests indicated no significant differences in faecal microbiota richness (Chao1), diversity (Shannon, Inverse Simpson), or Faith’s phylogenetic diversity (PD) between the two Lepus spp. species (Fig. S3), regardless of whether all four elevations were included or only the two where both species were found in sympatry (i.e. 2000m, 2500 m a.s.l.). However, NMDS based on beta diversity estimates highlighted a clear separation between the microbial communities of L. europaeus and L. timidus (Fig. 4 a, c). Consistently, the PERMANOVA analysis performed on samples collected from both species at 2000 and 2500 m a.s.l. which indicated significant differences in prokaryote community composition (PERMANOVA: unweighted Unifrac: R 2 = 0.20, p-value = 0.001; Bray-Curtis: R 2 = 0.27, p-value = 0.001; Table 1 ). These findings were further corroborated by the lack of differences in group dispersions across species and elevations (betadisper; p-value > 0.05 in all cases except L. timidus at 2000 m a.s.l. vs L. timidus at 2500 m a.s.l.). The clear differentiation in faecal prokaryotic communities between L. europaeus and L. timidus at 2000 and 2500 m a.s.l. was further confirmed by ANOSIM (Bray-Curtis: R = 0.886, p-value = 0.001; unweighted Unifrac R = 0.675, p-value = 0.001). Table 1 Results of PERMANOVA analysis comparing the faecal microbiota of L. europaeus and L. timidus at sympatric elevations. Microbial diversity was estimated using Bray-Curtis, Unweighted Unifrac (for 16S rRNA gene only) and Jaccard index (ITS2 locus only). For each microbial community and diversity index, statistical tests were performed by considering species, elevation and their interactions. Microbial community (marker) Diversity Index Variable R 2 F p-value Sig. Prokaryotes (16S rRNA gene) Bray-Curtis Species 0.273 15.207 0.001 *** Elevation 0.022 1.252 0.203 ns Species:Elevation 0.024 1.336 0.171 ns Residual 0.681 ns Unweighted Unifrac Species 0.200 9.903 0.001 *** Elevation 0.017 0.855 0.528 ns Species:Elevation 0.016 0.774 0.655 ns Residual 0.767 ns Fungi (ITS2 locus) Bray-Curtis Species 0.240 13.659 0.001 *** Elevation 0.052 2.967 0.007 ** Species:Elevation 0.041 2.332 0.022 * Residual 0.667 ns Jaccard Species 0.133 6.266 0.001 *** Elevation 0.033 1.534 0.046 * Species:Elevation 0.029 1.364 0.085 ns Residual 0.806 ns p-value ≤ 0.001: ***; p-value ≤ 0.01: **; p-value ≤ 0.05: *; ns: not significant. For the fungal community, NMDS clustering of L. timidus and L. europaeus faecal samples collected at sympatric elevations revealed distinct clusters associated with species (Fig. 4 b, d). This was supported by PERMANOVA analyses, which highlighted significant differences in the fungal community composition between L. timidus and L. europaeus (PERMANOVA: Jaccard: R 2 = 0.133, p-value = 0.001; Bray-Curtis R 2 = 0.240, p-value = 0.001), as well as across the two considered elevations (PERMANOVA: Jaccard: R 2 = 0.033, p-value = 0.046; Bray-Curtis R 2 = 0.052, p-value = 0.007). A weaker but still significant species-by-elevation interaction was observed (Jaccard: R 2 = 0.029, p-value = 0.085; Bray-Curtis R 2 = 0.041, p-value = 0.022; Table 1 ). Additional support for species-specific differences in mycobiota at sympatric elevations was provided by ANOSIM (Bray-Curtis: R = 0.438, p-value = 0.001; Jaccard R = 0.517, p-value = 0.001). Intraspecific variation in diversity along an elevational gradient NMDS clustering and PERMANOVA analysis of L. europaeus samples across the entire elevational gradient (Fig. 5 a, b; Tab. S3) provided limited evidence for an association between elevation and prokaryote community composition (PERMANOVA: Bray-Curtis: R 2 : 0.064, p-value = 0.011; Unifrac: R 2 : 0.053, p-value = 0.082). This result was also confirmed by the ANOSIM results which showed no significant differences across elevations (Bray-Curtis: R = 0.046, p-value = 0.075; unweighted Unifrac R = 0.032, p-value = 0.129). Conversely, the same ordination and statistical analysis of L. europaeus fungal communities revealed a clear association between elevation and composition of faecal mycobiota (PERMANOVA: Bray-Curtis R 2 = 0.153, p-value = 0.001. Jaccard R 2 = 0.098, p-value = 0.001; ANOSIM; Bray-Curtis: R = 0.244, p-value = 0.075; Jaccard: R = 0.3075, p-value = 0.129) (Tab. S3, Fig. S4b). Similarly, no association with elevation was detected for L. timidus prokaryote beta diversity estimates (PERMANOVA and ANOSIM: p-value > 0.05; Tab. S3; Fig. S4c). In contrast, mycobiota composition showed a significant association with elevation (PERMANOVA: Bray-Curtis R 2 = 0.118, p-value = 0.018. Jaccard R 2 = 0.079, p-value = 0.014; ANOSIM; Bray-Curtis: R = 0.339, p-value = 0.001; Jaccard: R = 0.360, p-value = 0.005; Tab. S3, Fig. S4d). Consistently, pairwise PERMANOVAs performed separately for L. europaeus and L. timidus revealed significant differences in fungal beta diversity between faecal samples collected at all elevation pairs, including those from neighbouring elevations (pairwise PERMANOVA: adjusted p-value < 0.05 for both species in all pairwise comparisons; Tab. S4). Conversely, no significant differences were observed in prokaryote communities across elevations for either species, including those between the most distant sites (pairwise PERMANOVA: adjusted p-value > 0.05 for both species in all pairwise comparisons; Tab. S4). Influence of environmental factors on Lepus spp. microbial communities No significant correlations were detected between prokaryote diversity and temperature, pre-sampling precipitation, or plant richness for either L. europaeus or L. timidus (Mantel test: Spearman's ρ adjusted p-values > 0.01; Tab. S5). Instead, fungal diversity in faecal samples of L. europaeus , but not L. timidus , was significantly associated with all three environmental variables: temperature (Mantel test; Bray-Curtis: Spearman's ρ = 0.288, adjusted p-values = 0.002; Jaccard: Spearman's ρ = 0.262, adjusted p-values = 0.001), pre-sampling total precipitation (Mantel test; Bray-Curtis: Spearman's ρ = 0.118, adjusted p-values = 0.006; Jaccard: Spearman's ρ = 0.255, adjusted p-values = 0.001) and plant richness (Mantel test; Bray-Curtis: Spearman's ρ = 0.123, adjusted p-values = 0.002; Jaccard: Spearman's ρ = 0.279, adjusted p-values = 0.001) (Tab. S5). Furthermore, several fungal ASVs were found to be correlated with pre-sampling precipitation (four ASVs, positively correlated), temperature (15 ASVs positively correlated; seven ASVs negatively correlated), and plant richness (five ASVs positively correlated; 17 ASVs negatively correlated) (Fig. 5 e). Functional diversity of prokaryote microbiota in Lepus spp. The Wilcoxon Rank-Sum on predicted metaCyc pathway abundances suggested there were significant differences in Shannon (functional) diversity between faecal samples of sympatric L. timidus and L. europaeus , with L. timidus showing higher Shannon estimates than L. europaeus (Fig. 6 a). However, no intraspecific difference in functional diversity estimates was found between elevations (Fig. 6 b). The clustering of faecal prokaryote communities using NMDS with beta diversity estimates (i.e. Bray-Curtis dissimilarity) on predicted metaCyc pathway abundances highlighted a clear separation between the two Lepus spp. species (PERMANOVA: Bray-Curtis R 2 = 0.498, p-value = 0.001; ANOSIM: Bray-Curtis R = 0.327, p-value = 0.001; Fig. 6 c). The clustering of Lepus spp. samples based on pathway abundances (Fig. 6 d) identified the species as the main clustering variable. Consistently, differential abundance testing comparing the predicted metaCyc pathway abundances across Lepus spp. species and elevations identified 75 pathways, all showing significant differences between L. europaeus at 2000 m a.s.l. and L. timidus at 2500 m a.s.l. (Fig. 6 d, Tab. S7). Of note, 61 out of these 75 pathways displayed significant differences in abundance between L. europaeus and L. timidus at 2000 m a.s.l. as well. Additionally, the contrast between L. europaeus and L. timidus at 2500 m a.s.l. identified 24 metaCyc pathways with significant differences between the two species, the majority of which (70.8%) were also included in the contrast between L. europaeus at 2000 m a.s.l. and L. timidus at 2500 m a.s.l. Differential pathways were involved in different classes of biological processes (e.g. superclasses), with Cofactor, Carrier, and Vitamin Biosynthesis (18 pathways, 24% of differential pathways), Fatty Acid and Lipid Biosynthesis (7 pathways, 9% of differential pathways), Carbohydrate Degradation (7 pathways, 9% of differential pathways), Amino Acid Degradation (5 pathways, 7% of differential pathways), and Nucleoside and Nucleotide Degradation (5 pathways, 7% of differential pathways) being the most represented (Fig. 6 d, Tab. S7). Of note, while most of these superclasses displayed a relatively balanced proportion of pathway abundances between the two species, all differential pathways within the superclasses Fatty Acid and Lipid Biosynthesis and Amino Acid Degradation were significantly enriched in L. timidus compared to L. europaeus . Discussion Gut microbiota diversity and composition are known to significantly impact animal health and survival [ 56 , 57 ]; thus, investigating gut microbial communities in species potentially threatened by rapid biotic and abiotic environmental changes may enhance our understanding of their conservation status and extinction risk. In this study, we characterized the gut microbiota of the mountain ( L. timidus ) and European brown ( L. europaeus ) hares using field-collected faecal samples to investigate the consequences of climate change-driven sympatry on their gut prokaryote and fungal communities. Our data confirms the effectiveness of non-invasive samples for the monitoring of gut microbial diversity and composition in elusive species such as Lepus spp. Unexpectedly, despite overlapping in the upper 1000 m of their vertical distribution, our results demonstrated that L. europaeus and L. timidus have distinct faecal microbiota and mycobiota profiles. For the first time, we also identified a significant association between elevation and fungal diversity in both species. The composition of L. timidus and L. europaeus faecal microbiotas are distinct The prokaryote faecal microbiota of the two Lepus spp. were characterized by comparable alpha diversity estimates (Fig. S3), but marked differences in community structure (Fig. 4 ), with up to 27% of observed variation associated with host species (Table 1 ). Accordingly, prokaryote communities displayed distinct taxonomic profiles (Fig. 2 , 3 ), with only a limited number of shared abundant ASVs (i.e. ~6% ASVs accounting for ~ 21% sequence reads; Fig. S1 ). L. europaeus faecal samples were primarily dominated by the phylum Firmicutes, and rich in Bacteroidota and Spirochaetota. Instead, the most abundant phylum detected in L. timidus was Proteobacteria, followed by Bacteroidota and Acidobacteriota. A dominance of Firmicutes and Bacteroidota has already been reported for L. europaeus faecal samples [ 22 ] as well as L. granatensis [ 58 ] and L. americanus [ 59 ]. Instead, a high relative abundance of Proteobacteria, as seen in L. timidus , was previously reported for L. sinensis [ 60 ]. The clear differentiation observed between the microbial communities of the two species, particularly at higher elevations where both species were present, suggests that any contaminants from environmental sources may have had only marginal effects on the diversity or composition of faecal bacterial communities. While Firmicutes, Bacteroidota and Proteobacteria are both very commonly found in the mammalian gut, the observed differences in prokaryote composition may reflect co-evolution of host species and associated microbiota after speciation and/or adaptation to different ecological niches [ 61 ]. For example, L. europaeus is primarily a steppe species (with the earliest fossil evidence in the Alpine region dating to approximately 8,000 years ago [ 62 ]), while L. timidus is a boreal-adapted species [ 63 ]. While there are no other studies addressing microbial diversity in sympatric Leporidae, the greater impact of phylogeny compared to other environmental factors, including diet, has been already documented in other mammals [ 64 , 65 , 66 ]. Therefore, despite their coexistence at higher elevations in our Alpine study area, in keeping with their evolutionary history, these species probably exploit distinct resources, with L. europaeus grazing on the high fat parts of weeds/grasses and various crop types richer in fats and proteins, as noted by Reichlin et al. [ 67 ] and Schai-Braun et al. [ 68 ]), while L. timidus may exploit significantly more ligneous plants and Ericaceae [ 69 ]. Future studies combining longitudinal vegetation surveys, diet metabarcoding, and metataxonomic investigations are needed to clarify the links existing between seasonal vegetation variation, diet and gut microbial diversity in the two species. Although the effects of seasonal dietary variation over the course of the year on gut microbiota are still largely unknown for these species, the impact of dietary niche on the gut microbiota has been well-documented for many mammal species including humans [ 70 ], laboratory mice [ 71 ] and wild mammals (see review in: Ley et al. [ 58 ]; Alessandri et al. [ 72 ]). Additionally, the contrasting taxonomic composition of L. europaeus and L. timidus gut microbiota implies different strategies for energy resorption adopted by the two microbial communities [ 64 , 73 , 74 ]. In fact, our PICRUSt2 analysis supports the hypothesis that the two host-associated microbiotas have alternative biosynthetic and metabolic potential, as indicated by the differences in the number and relative abundance of predicted functions (Fig. 6 ). Although a high abundance of Ericaceae has been reported in the L. timidus diet [ 69 ], we found no differences in metabolic pathways associated with the degradation of phenolic compounds which are abundant in these plants. Interestingly, all differential pathways involved in fatty acid / lipid biosynthesis and amino acid degradation were significantly enriched in L. timidus , which may again reflect differences in diet between the two species [ 68 , 69 ]. We also speculate that since host-mediated biosynthesis of fatty acids and lipids support metabolic processes that generate heat, thereby aiding in thermoregulation and overall survival in colder climates [ 75 ], microbial mediated biosynthesis of these molecules may help L. timidus to cope with the cold temperatures typically found in its natural home range. Overall, these results suggest that the composition of the gut microbiota may be associated with the adaptation of L. timidus to its boreal habitat. However, again, additional studies on the actual diet of the individuals sampled in this study are needed to confirm this hypothesis. While the prevalence and relative abundance of Proteobacteria in L. timidus were consistent across faecal pellets, most L. europaeus samples were characterized by relatively low abundance of this phylum. However, 15% of investigated faecal samples displayed unexpectedly high abundances (e.g. > 80% of sequence reads) of Proteobacteria, mostly belonging to Enterobacteraceae, Erwiniaceae and Oxalobacteraceae (Fig. S2 ). Enterobacteriaceae is the most studied family of these three and includes both commensal bacteria contributing to the maintenance of the gut anaerobic environment, the production of secondary metabolites (e.g. vitamins) and the protection against gut pathogens, as well as opportunistic disease-causing pathogens [ 76 , 77 ]. Although Enterobacteriaceae are commonly found in low abundance in the gut microbiota of healthy mammals, their proliferation in humans is also considered to be a biomarker of gut dysbiosis, which has been associated with several inflammatory bowel diseases [ 78 ]. Because only non-invasive faecal samples were collected, we do not know if the observed abundance of Proteobacteria in a limited, but not negligible, number of L. europaeus samples could be attributable to a state of inflammation of the host. Changes in microbiota composition or diet may also lead to favourable environmental conditions for Proteobacteria (e.g. variation in carbon sources or decreased hypoxia, among others [ 78 ]). This unusually high abundance of Proteobacteria should be further investigated, since a disrupted gut microbiota, or microbiota consisting of a number of potentially pathogenic taxa could impact the health of Lepus individuals. Furthermore, future studies should also plan detailed investigations targeting specific Enterobacteriaceae taxa as biomarkers, which could be useful for health monitoring of these species. Similarly, both diversity and taxonomic composition of fungal taxa between L. europaeus and L. timidus diverged significantly, (Fig. 4 ), with up to 24% of observed variation, associated with host species (Table 1 ). The fungal classes Tremellomycetes and Cystobasidiomycetes were significantly enriched in mountain hare gut microbiota, while higher abundances of Ustilaginomycetes and Pezizomycetes were detected in the European brown hare. The class Tremellomycetes is a nutritionally heterogeneous group comprising saprotrophs, animal parasites, and fungicolous species previously isolated from various habitats including alpine soil [ 79 ], subalpine grasses [ 80 ], snow collected at high elevations [ 81 ] and the Antarctic Polar Plateau [ 82 ]. Due to the ability of some Tremellomycetes to survive in mountainous or cold habitats, we hypothesize that the high abundance of Tremellomycetes found in L. timidus resulted from these animals having more opportunities to feed on these fungi or have contact with them at higher elevations. However, this fungal class has also been linked to high lipid production [ 83 ], which may contribute to energy homeostasis and thermoregulation [ 84 ] in L. timidus . On the other hand, Ustilaginomycetes, found here to be more abundant in L. europaeus , are well-known to infect vascular plants, especially grass families like Poaceae [ 85 ]. Indeed, L. europaeus are primarily grazers and have a diet that heavily relies on grasses [ 69 ], possibly explaining the higher abundance of this fungal class, as well as Pezizomycetes (phylum Ascomycota) in their gut compared to L. timidus . Again, molecular diet analysis of the faecal pellets would allow us to confirm that the patterns in fungal diversity and composition are mainly due to differences in feeding preferences. In this study, host species identification was conducted using mitochondrial DNA (mtDNA) barcoding, which is currently regarded as the most reliable approach for species assignment from faecal samples, owing to the high copy number of mitochondrial genomes per cell and their relatively high mutation rate which allow for good taxonomic discrimination, even between closely related species [ 86 , 87 ]. However, because mtDNA barcoding does not allow discrimination between hybrids and their maternal parental species, we cannot exclude the possibility that a small number of the analysed samples may correspond to hybrid individuals. However, the occurrence in northern European and Alpine regions has been estimated to occur at a frequency of only 2–5% [ 34 , 88 , 89 ]. Therefore, while hybridization is frequently associated with substantial shifts in the composition and diversity of gut microbial communities [ 90 , 91 , 92 ], we do not know yet how the observed differentiation in composition and predicted functional potential between L. europaeus and L. timidus could be translated into the gut microbiota of their hybrids. L. europaeus and L. timidus hybrids might exhibit transgressive microbiome configurations, not directly resembling the composition and the predicted functional potential between parental species as seen for hybrids between subspecies Mus musculus musculus × M. m. domesticus [ 90 ], intermediate microbiome configurations (e.g. Centropyge spp.; [ 91 ]), or asymmetric, parent-biased microbiomes as seen for Cervus elaphus [ 92 ]. Therefore, at present we are unable to predict the most likely microbiota configuration of Lepus hybrids, or its influence on host nutritional behavior, disease susceptibility, or adaptive potential to changing environments. Genotyping of faecal samples could resolve the issue by identifying hybrid individuals; however, many more samples would be needed since the success rate of genotyping from hare faecal samples is only about 50%, and the occurrence of hybrids is so low and unpredictable. Within-species variation in microbial diversity across an elevational gradient Elevation did not appear to affect the alpha diversity of prokaryote and fungal communities hosted by the two Lepus species, and only marginally impacted (Tab. S3) beta diversity estimates of L. europaeus prokaryote gut communities. The resilience of the hare gut prokaryote microbiota to variation in elevation and associated environmental variables was unexpected based on studies performed on other mammal species like wild house mice [ 93 ], pikas [ 94 ], ungulates [ 95 ]; macaques, humans and domestic dogs [ 96 ]. However, the beta diversity of the prokaryote community structure was in striking contrast to that observed in the gut fungal communities of both species; in fact, marked differences in fungal microbiota composition were detected using both Bray-Curtis (e.g. accounting for ASV abundances) and Jaccard (e.g. based on ASV presence/absence) indices, indicating that variation in fungal composition across elevations reflected changes in both the shared abundant and rare components of the environmental fungal community. Furthermore, we found elevation to be a significant driver of fungal diversity even in pairwise comparisons between neighbouring elevations (e.g. 1000m vs 1500m a.s.l. for L. europaeus and 2000 m vs 2500 m for L. timidus ). Consistently, Mantel tests highlighted low but significant correlations between the composition of overall fungal communities and pre-sampling precipitation, temperature and plant richness. Moreover, several fungal ASVs in L. europaeus , including seven matching the coprophilous genus Sporormiella (phylum Ascomycota), that uses herbivore dung as a primary substrate [ 97 ], was significantly associated with these environmental variables. Interestingly, variation in the relative abundance of Sporomiella in the gut of the lagomorph Ochotona curzoniae (i.e. plateau pikas) captured at different elevations was also reported by Tang et al. [ 98 ]. However, only a few fungal taxa are considered gut residents [ 99 , 100 ] and environmental fungi ingested with diet are almost certainly more exposed to abiotic variation. Consistently, differences in soil fungal communities driven by elevation have been reported frequently [ 79 , 101 , 102 ]. Additionally, the survival of many fungal taxa relies on interactions with plants and soil [ 103 , 104 , 105 ]. Therefore, the high number of taxa common to both species (80% of generated sequence reads, Fig. S1 ), together with the marked association with elevation and related environmental variables (Fig. 5 , Tab. S5), indicate a strong contribution of diet (and possibly geophagy) to the gut fungal assembly in these two Lepus species, as suggested by Grieneisen et al. [ 106 ] for a primate hybrid zone (see also [ 107 ]). The clear differentiation between prokaryote and fungal communities in their resilience to environmental parameters could have interesting implications in relation to the observed elevation shifts in the distribution of L. europaeus and the associated potential reduction in habitat availability for L. timidus . The lack of differentiation in prokaryote communities across elevations suggests that gut functionality was being maintained along the entire gradient. Our data also suggested that L. europaeus and L. timidus had different feeding strategies (e.g. selecting different plant species) maintained along the same gradient, allowing them to co-exist in the same areas. However, we expect a selective feeding strategy to represent a potential vulnerability for L. timidus , particularly if climate change continues to alter the distribution of plant taxa in the L. timidus diet. Furthermore, the unique components of the L. timidus gut microbiota biodiversity also risk disappearing, with unknown impacts on the Alpine ecosystem. Future studies addressing the health impacts of diet, geophagy and other environmental factors across elevations are needed to disentangle the associations of these factors with the taxonomic and functional diversity of L. europaeus and L. timidus in the Alps, and consequently, of the implications for their survival during climate warming. Conclusions Our findings demonstrate that sympatric populations of L. europaeus and L. timidus maintained distinct, species-specific prokaryote and fungal gut microbiota in terms of taxonomy and structure, but not diversity, at least within our study area. Our data act as a baseline for future monitoring of the process of adaptation of these species to climate warming and associated environmental changes. Future studies should prioritize the characterization of gut prokaryote and fungal composition in hybrids of the two lagomorph species across multiple populations across the Alps. Moreover, future research should explore how climate-driven changes in soil microbiota and food plant composition may shape the future diversity and structure of host-associated gut microbiota in these Lagomorph species. Declarations Ethics approval and consent to participate Not applicable. Competing interests The authors declare no competing interests. Consent for publication Not applicable. Additional information The authors declare no competing interests Funding The EUREGIO project: MICROVALU - Evaluating microbiodiversity in alpine pastures (Project ID: IPN94), awarded to PI, JS and HCH, was funded by the "Euregio Tirolo-Alto Adige-Trentino" Interregional Project Network. This study was also partially carried out with funding to HCH at the Fondazione E. Mach (Project BIOALPEC) under the National Biodiversity Future Centre (NBFC) Project (code CN_00000033, Concession Decree No. 1034 of 17 June 2022 adopted by the Italian Ministry of University and Research, CUPD43C22001280006), funded under the National Recovery and Resilience Plan (NRRP), Mission 4 Component 2 Investment 1.4 - Call for tender No. 3138 of 16 December 2021, rectified by Decree n.3175 of 18 December 2021 of Italian Ministry of University and Research funded by the European Union – NextGenerationEU. This manuscript reflects only the authors’ views and opinions, neither the European Union nor the European Commission can be considered responsible for them. Author Contribution B.C., H.C.H., and G.G. conducted the sampling. L.M., B.C. and G.G. completed the laboratory analyses. Data analyses were provided by L.M. and G.G., with support from N.P., T.R., J.S., P.I. and F.N.M.; L.M., H.C.H. and G.G. drafted the manuscript. All authors edited and approved the final manuscript. Acknowledgement The authors wish to thank the Fondazione E. Mach for access to facilities, and the staff of the Sequencing and Genotyping Platform for their outstanding support. Data Availability Sanger sequences have been deposited at NCBI GenBank with accession numbers PX122685 - PX122779. Sanger sequences of the mitochondrial D-loop generated in our study were made available to the editor and reviewers with the uploaded file: Submission2992130.txt.gz. The raw amplicon-sequencing data has been deposited at NCBI Sequence Read Archive (SRA) under the BioProject ID PRJNA1304890. Reviewers can access BioProject and associated SRA metadata at [https://dataview.ncbi.nlm.nih.gov/object/PRJNA1304890?reviewer=rsuoanlfj6uo274kij1dik96rn] . References Lenoir, J., Gégout, J. C., Marquet, P. A., de Ruffray, P. & Brisse, H. A. significant upward shift in plant species optimum elevation during the 20th century. Science 320 , 1768–1771 (2008). Harsch, M. A., Hulme, P. E., McGlone, M. S. & Duncan, R. P. Are treelines advancing? A global meta-analysis of treeline response to climate warming. Ecol. Lett. 12 , 1040–1049 (2009). Noce, S., Cipriano, C. & Santini, M. Altitudinal shifting of major forest tree species in Italian mountains under climate change. Front. For. Glob. Change 6 , 1250651; https://doi.org/10.3389/ffgc.2023.1250651 (2023). Hickling, R., Roy, D. B., Hill, J. K., Fox, R. & Thomas C. D. The distributions of a wide range of taxonomic groups are expanding polewards. Glob. Change Biol. 12 , 450-455 (2006). Popy, S., Bordignon L. & Prodon, R. A. weak upward elevational shift in the distributions of breeding birds in the Italian Alps. J. Biogeogr. 37 , 57-67 (2010). Büntgen U. et al. Elevational range shifts in four mountain ungulate species from the Swiss Alps. Ecosphere 8 , e01761. https://doi.org/10.1002/ecs2.1761 (2017). Schai-Braun, S. C., Jenny, H., Ruf, T. & Hackländer, K. Temperature increase and frost decrease driving upslope elevational range shifts in Alpine grouse and hares. Glob. Change Biol. 27 , 6602–6614 (2021). Rehnus, M. Der Schneehase in den Alpen: Ein Überlebenskünstler mit ungewisser Zukunft . Bristol-Stiftung, Zürich; Haupt, Bern (2013). Stokes, A. W., Hofmeester, T. R., Thorsen, N. H., Odden, J., Linnell, J. D. C. & Pedersen, S. Altitude, latitude and climate zone as determinants of mountain hare ( Lepus timidus ) coat colour change. Ecol. Evol. 13 , e10548. https://doi.org/10.1002/ece3.10548 (2023). Sliwinski, K. et al. Habitat requirements of the European brown hare ( Lepus europaeus Pallas 1778) in an intensively used agriculture region (Lower Saxony, Germany). BMC Ecology , 19, 31. https://doi.org/10.1186/s12898-019-0247-7 (2019). Levänen, R., Thulin, C.-G., Spong, G., & Pohjoismäki J. L. O. Widespread introgression of L. timidus genes into Fennoscandian brown hare populations. PLoS ONE 13 , e0191790. https://doi.org/10.1371/journal.pone.0191790 (2018). Jansson, G. & Pehrson, Å. The recent expansion of the brown hare ( Lepus europaeus ) in Sweden with possible implications to the mountain hare ( L. timidus ). Eur. J. Wildl. Res. 53 , 125–130 (2007). La Morgia, V. et al. Global warming is promoting the rapid invasion of the L. timidus range by the European hare in the Alps. Biodivers. Conserv. 32 , 3875–3891 (2023). Edwards, J. E. et al. Domesticated equine species and their derived hybrids differ in their faecal microbiota. Anim. Microbiome 2 , Article 27. https://doi.org/10.1186/s42523-020-00027-7 (2020). Wei, L. et al. Hybridization altered the gut microbiota of pigs. Front. Microbiol. 14 , 1177947. https://doi.org/10.3389/fmicb.2023.1177947 (2023). Rowland, I. et al. Gut microbiota functions: metabolism of nutrients and other food components. Eur. J. Nutr. 57 , 1–24 (2018). Sittipo, P. et al. The function of gut microbiota in immune-related neurological disorders: a review. J. Neuroinflammation 19 , 154. https://doi.org/10.1186/s12974-022-02510-1 (2022). Kasarello, K., Cudnoch-Jędrzejewska, A. & Czarzasta, K. Communication of gut microbiota and brain via immune and neuroendocrine signaling. Front. Microbiol. 14 , 1118529. https://doi.org/10.3389/fmicb.2023.1118529 (2023). Moeller, A. H. et al. Experimental evidence for adaptation to species-specific gut microbiota in house mice. mSphere 4 , e00387-19. https://doi.org/10.1128/mSphere.00387-19 (2019). Hauffe, H. C. & Barelli, C. Conserve the germs: the gut microbiota and adaptive potential. Conserv. Genet. 20 , 19–27 (2019). Rojas, C. A., Ramírez-Barahona, S., Holekamp, K. E. & Theis, K. R. Host phylogeny and host ecology structure the mammalian gut microbiota at different taxonomic scales. Anim. Microbiome 3 , 33. https://doi.org/10.1186/s42523-021-00094-4 (2021). Neha, S. A., Hanson, J. D., Wilkinson, J. E., Bradley, R. D. & Phillips, C. D. Impacts of host phylogeny, diet, and geography on the gut microbiome of rodents. PLoS ONE 20 , e0316101. https://doi.org/10.1371/journal.pone.0316101 (2025). Barelli, C., Albanese, D., Stumpf, R. M., Asangba, A., Donati, C., Rovero, F. & Hauffe H. C. The gut microbiota communities of wild arboreal and ground-feeding tropical primates are affected differently by habitat disturbance. mSystems 5 , e00061-20. https://doi.org/10.1128/mSystems.00061-20 (2020). Amato, K. R. et al. The gut microbiota appears to compensate for seasonal diet variation in the wild black howler monkey ( Alouatta pigra ). Microb. Ecol. 69 , 434–443 (2015). Stalder, G. L. et al. Gut microbiota of the European brown hare ( Lepus europaeus ). Sci. Rep. 9 , 2738. https://doi.org/10.1038/s41598-019-39638-9 (2019). Zhang, L., Zhan, H., Xu, W., Yan, S. & Ng, S. C. The role of gut mycobiome in health and diseases. Therap. Adv. Gastroenterol. 14 , 17562848211047130. https://doi.org/10.1177/17562848211047130 (2021). Zandonai, A. et al. Six years of high-resolution climatic data collected along an elevation gradient in the Italian Alps. Sci. Data 11 , 751. https://doi.org/10.1038/s41597-024-03580-x (2024). Hilpold, A., Seeber, J., Fontana, V., Niedrist, G., Rief, A., Steinwandter, M., Tasser, E. & Tappeiner, U. Decline of rare and specialist species across multiple taxonomic groups after grassland intensification and abandonment . Biodivers. Conserv. 27 , 3729–3744 (2018). Praeg, N., Pauli, H. & Illmer, P. Microbial diversity in bulk and rhizosphere soil of Ranunculus glacialis along a high-alpine altitudinal gradient. Front. Microbiol. 10 , 1429. https://doi.org/10.3389/fmicb.2019.01429 (2019). Galla, G. et al. Mock community as an in situ positive control for amplicon sequencing of microbiotas from the same ecosystem. Sci. Rep. 13 , 4056. https://doi.org/10.1038/s41598-023-30916-1 (2023). Wilkinson, G. S. & Chapman, A. M. Length and sequence variation in evening bat D-loop mtDNA. Genetics 128 , 607–617 (1991). Melo-Ferreira, J., Boursot, P., Randi, E., Kryukov, A., Suchentrunk, F., Ferrand, N. & Alves, P. C. The rise and fall of Lepus timidus during Pleistocene glaciations: expansion and retreat with hybridization in the Iberian Peninsula. Mol. Ecol. 16 , 605–618 (2007). Pecchioli, E. Variabilità e struttura genetica in Lepus timidus : aspetti micro e macrogeografici, filogenetici e adattativi (unpublished PhD thesis, Università degli Studi di Ferrara) (2010). Zachos, F. E., Ben Slimen, H., Hackländer, K., Giacometti, M. & Suchentrunk, F. Regional genetic in situ differentiation despite phylogenetic heterogeneity in Alpine Lepus timidus . J. Zool. 282 , 47–53 (2010). Kumar S., Stecher G., Li M., Knyaz C, & Tamura K. MEGA X: Molecular Evolutionary Genetics Analysis across computing platforms. Mol. Biol. Evol. 35, 1547-1549 (2018). Klindworth, A. et al. Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. Nucleic Acids Res. 41 , e1 (2013). Walters, W. et al. Improved prokaryotic 16S rRNA gene (V4 and V4-5) and fungal internal transcribed spacer marker gene primers for microbial community surveys. mSystems 1 , e00009-15. https://doi.org/10.1128/mSystems.00009-15 (2015). Ihrmark, K. et al. New primers to amplify the fungal ITS2 region – evaluation by 454-sequencing of artificial and natural communities. FEMS Microbiol. Ecol. 82 , 666–677 (2012). White, T. J., Bruns, T., Lee, S. & Taylor, J. Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics. In: PCR Protocols: A Guide to Methods and Applications (eds. Innis, M. A., Gelfand, D. H., Sninsky, J. J. & White, T. J.) 315–322 (Academic Press, San Diego) (1990). Callahan, B. J. et al. DADA2: High-resolution sample inference from Illumina amplicon data. Nat. Methods 13 , 581–583 (2016). Quast, C. et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 41 , D590–D596 (2012). Nilsson, R. H. et al. The UNITE database for molecular identification of fungi: handling dark taxa and parallel taxonomic classifications. Nucleic Acids Res. 47 , D259–D264 (2019). Bodenhofer, U., Bonatesta, E., Horejs-Kainrath, C. & Hochreiter, S. msa: an R package for multiple sequence alignment. Bioinformatics 31 , 3997–3999 (2015). Schliep, K. P. phangorn: phylogenetic analysis in R. Bioinformatics 27 , 592–593 (2011). Davis, N. M., Proctor, D. M., Holmes, S. P., Relman, D. A. & Callahan, B. J. Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data. Microbiome 6 , 226. https://doi.org/10.1186/s40168-018-0605-2 (2018). McMurdie, P. J. & Holmes, S. phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One 8 , e61217. https://doi.org/10.1371/journal.pone.0061217 (2013). Liu, C., Cui, Y., Li, X. & Yao, M. microeco: An R package for data mining in microbial community ecology. FEMS Microbiol. Ecol. 97 , fiaa255 (2021). Wickham, H. ggplot2: Elegant Graphics for Data Analysis . Springer-Verlag, New York (2016). Neuwirth, E. RColorBrewer: ColorBrewer Palettes . R package version 1.1-3, available at https://CRAN.R-project.org/package=RColorBrewer (2022). Garnier, S. et al. viridis(Lite) – Colorblind-Friendly Color Maps for R. R package version 0.6.5. Available at: https://sjmgarnier.github.io/viridis/ (2024). Segata, N. et al. Metagenomic biomarker discovery and explanation. Genome Biol. 12 , R60. https://doi.org/10.1186/gb-2011-12-6-r60 (2011). Oksanen, J. et al. vegan: Community Ecology Package. R package version 2.6-4, available at https://CRAN.R-project.org/package=vegan (2022). Douglas, G. M. et al. PICRUSt2 for prediction of metagenome functions. Nat. Biotechnol. 38 , 685–688 (2020). Metsalu, T. & Vilo, J. ClustVis: a web tool for visualizing clustering of multivariate data using Principal Component Analysis and heatmap. Nucleic Acids Res. 43 , W566–W570 (2015). The GIMP Development Team. GIMP (Version 2.10.18) [Computer software]. https://www.gimp.org (2020). Tuddenham, S. & Sears, C. L. The intestinal microbiome and health. Curr. Opin. Infect. Dis. 28 , 464–470 (2015). Zhan, G. et al. Abnormal gut microbiota composition contributes to cognitive dysfunction in SAMP8 mice. Aging 10 , 1257–1267. https://doi.org/10.18632/aging.101464 (2018). Martínez-Haro, M., Triadó-Margarit, X., Mateo, R., Viñuela, J. A. & Casamayor, E. O. Taxonomic and functional fingerprints in the gut microbiota of Iberian hare populations ( Lepus granatensis ) inhabiting areas with organic and conventional farming. Environ. Pollut. 369 , 125862 (2025). Stothart, M. R. et al. Population dynamics and the microbiome in a wild boreal mammal: The snowshoe hare cycle and impacts of diet, season and predation risk. Mol. Ecol. 34 , e17629. https://doi.org/10.1111/mec.17629 (2025). Cao, Z. et al. Comparison and association of winter diets and gut microbiota using trnL and 16S rRNA gene sequencing for three herbivores in Taohongling, China. Glob. Ecol. Conserv. 53 , e03041. https://doi.org/10.1016/j.gecco.2024.e03041 (2024). Koskella, B. & Bergelson, J. The study of host–microbiome (co)evolution across levels of selection. Philos. Trans. R. Soc. B 375 , 20190604 (2020). Bock, A. Lepus europaeus (Lagomorpha: Leporidae). Mamm. Species 52 , 125–142. https://doi.org/10.1093/mspecies/seaa010 (2020). Angerbjörn, A. & Flux, J. E. C. Lepus timidus . Mamm. Species 495 , 1–11 (1995). Ley, R. E. et al. Evolution of mammals and their gut microbes. Science 320 , 1647–1651 (2008). Knowles, S. C. L., Eccles, R. M. & Baltrūnaitė, L. Species identity dominates over environment in shaping the microbiota of small mammals. Ecol. Lett. 22 , 826–837 (2019). Anders, J. L. et al. Comparing the gut microbiome along the gastrointestinal tract of three sympatric species of wild rodents. Sci. Rep. 11 , 19929. https://doi.org/10.1038/s41598-021-99379-6 (2021). Reichlin, T., Klansek, E. & Hackländer, K. Diet selection by hares ( Lepus europaeus ) in arable land and its implications for habitat management. Eur. J. Wildl. Res. 52 , 109–118 (2006). Schai-Braun, S. C. et al . The European hare ( Lepus europaeus ): a picky herbivore searching for plant parts rich in fat. PLoS ONE 10 , e0134278. https://doi.org/10.1371/journal.pone.0134278 (2015). Schai-Braun, S. C. et al. Diet niches of the Alpine Lepus timidus ( Lepus timidus varronis ) and European hare ( Lepus europaeus ) living in coexistence in the Alps. Mamm. Biol. 105 , 139–145 (2025). David, L. A. et al . Diet rapidly and reproducibly alters the human gut microbiome. Nature 505 , 559–563 (2014). Carmody, R. N. et al. Diet dominates host genotype in shaping the murine gut microbiota. Cell Host Microbe 17 , 72–84 (2015). Alessandri, G., Rizzo, S. M., Ossiprandi, M. C., van Sinderen, D. & Ventura, M. Creating an atlas to visualize the biodiversity of the mammalian gut microbiota. Curr. Opin. Biotechnol. 73 , 28–33 (2022). Zoetendal, E. G. & de Vos, W. M. Effect of diet on the intestinal microbiota and its activity. Curr. Opin. Gastroenterol. 30 , 189–195 (2014). Sommer, F. et al. The gut microbiota modulates energy metabolism in the hibernating brown bear Ursus arctos . Cell Rep. 14 , 1655–1661 (2016). Xu, Z. et al. Cold-induced lipid dynamics and transcriptional programs in white adipose tissue. BMC Biol. 17 , 74. https://doi.org/10.1186/s12915-019-0693-x (2019). Vázquez-Baeza, Y. et al. Dog and human inflammatory bowel disease rely on overlapping yet distinct dysbiosis networks. Nat. Microbiol. 1 , 16177 (2016). Lupp, C. et al . Host-mediated inflammation disrupts the intestinal microbiota and promotes the overgrowth of Enterobacteriaceae . Cell Host Microbe 2 , 119–129 (2007). Shin, N. R., Whon, T. W. & Bae, J.-W. Proteobacteria: microbial signature of dysbiosis in gut microbiota. Trends Biotechnol. 33 , 496–503 (2015). Praeg, N. et al. The role of land management and elevation in shaping soil microbial communities: Insights from the Central European Alps. Soil Biol. Biochem. 150 , 107951 (2020). Mouhamadou, B. et al. Differences in fungal communities associated to Festuca paniculata roots in subalpine grasslands. Fungal Divers. 47 , 55–63 (2011). Irwin, N. A., Twynstra, C. S., Mathur, V. & Keeling, P. J. The molecular phylogeny of Chionaster nivalis reveals a novel order of psychrophilic and globally distributed Tremellomycetes (Fungi, Basidiomycota). PLoS ONE 16 , e0247594. https://doi.org/10.1371/journal.pone.0247594 (2021). Stoppiello, G. A., Coleine, C., Moeller, R., Ripa, C., Billi, D., & Selbmann, L. Seasonality is the main determinant of microbial diversity associated to snow/ice around Concordia Station on the Antarctic Polar Plateau. Biology 12 , 1193. https://doi.org/10.3390/biology12091193 (2023). Kurtzman, C. P., Fell, J. W., & Boekhout, T. The Yeasts: A Taxonomic Study. Fifth Edition (2011). Olsen, L., Thum, E., & Rohner, N. Lipid metabolism in adaptation to extreme nutritional challenges. Develop. Cell. 56 (2021). Bauer, R., Begerow, D., & Oberwinkler, F. Ustilaginomycotina: The true smut fungi. Tree of Life Web Project available at: http://tolweb.org/Ustilaginomycotina/20530/2008.01.23 (2006). Murphy, M. A., Waits, L. E., & Kendall, K. C. Quantitative evaluation of fecal drying methods for brown bear DNA analysis. Wildl. Soc. Bull. 28 , 951–957 (2000). Elyasigorji, Z., Izadpanah, M., Hadi, F., & Zare, M. Mitochondrial genes as strong molecular markers for species identification. Nucleus 66 , 81–93. https://doi.org/10.1007/s13237-022-00393-4 (2022). Thulin, C.-G., Stone, J., Tegelström, H., & Walker, C. W. Species assignment and hybrid identification among Scandinavian hares Lepus europaeus and L. timidus . Wildl. Biol., 12 , 29–38 (2006). Schai-Braun, S. C., Schwienbacher, S., Smith, S., & Hackländer, K. Coexistence of European hares and Alpine mountain hares in the Alps: what drives the occurrence and frequency of their hybrids? J. Zool. 320 (1), 15–25 (2023). Wang, J. et al. Analysis of intestinal microbiota in hybrid house mice reveals evolutionary divergence in a vertebrate hologenome. Nat. Commun. 6:6440. https//:doi:10.1038/ncomms7440 (2015) Huggett, M. J. et al. Gut microbial communities of hybridising pygmy angelfishes reflect species boundaries. Commun. Biol. 6 , 542 (2023). Wei, L. et al. Hybridization alters red deer gut microbiome and metabolites. Front. Microbiol. 15 , 1387957 (2024). Suzuki, T. A., Martins, F. M., & Nachman, M. W. Altitudinal variation of the gut microbiota in wild house mice. Mol. Ecol. 28 (2019). Li, H., Zhou, R., Zhu, J., Huang, X., & Qu, J. Environmental filtering increases with elevation for the assembly of gut microbiota in wild pikas. Microb. Biotechnol. 12 (2019). Wang, X., et al. High-altitude drives the convergent evolution of alpha diversity and indicator microbiota in the gut microbiomes of ungulates. Front. Microbiol. 13 https://doi.org/10.3389/fmicb.2022.953234 (2022). Zhao, J. et al. Diet and high altitude strongly drive convergent adaptation of gut microbiota in wild macaques, humans, and dogs to high altitude environments. Front. Microbiol. 14 https://doi.org/10.3389/fmicb.2023.1067240 (2023). Perrotti, A. G., & van Asperen, E. Dung fungi as a proxy for megaherbivores: opportunities and limitations for archaeological applications. Veg. Hist. Archaeobot. 28 , 93–104 (2019). Tang, X., Zhang, L., Ren, S., Zhao, Y., Liu, K., & Zhang, Y. Stochastic processes derive gut fungi community assembly of plateau pikas ( Ochotona curzoniae ) along altitudinal gradients across warm and cold seasons. J. Fungi 9 https://doi.org/10.3390/jof9101032 (2023). Auchtung, T. A. et al. Investigating colonization of the healthy adult gastrointestinal tract by fungi. mSphere 3 https://doi.org/10.1128/mSphere.00092-18 (2018). Suhr, M. J., & Hallen-Adams, H. E. The human gut mycobiome: pitfalls and potentials - a mycologist’s perspective. Mycologia 107 , 1057-73 (2015). Li, J., Chen, L., Wang, H., Ouyang, S., Liu, X. & Lu, J. Pattern and drivers of soil fungal community along elevation gradient in the Abies georgei forests of Segila mountains, Southeast Tibet. Glob. Ecol. Conserv. 39 , e02446. https://doi.org/10.1016/j.gecco.2022.e02291 (2022). Vélez-Martínez, G. A. et al. Soil bacteria and fungi communities are shaped by elevation influences in Colombian forest and páramo natural ecosystems. Int. Microbiol. 27 , 377-391 (2024). Dean, R. et al. The Top 10 fungal pathogens in molecular plant pathology. Mol. Plant Pathol. 13 , 414–430 (2012). Rodriguez, R. J., White, J. F. Jr., Arnold, A. E., & Redman, R. S. Fungal endophytes: diversity and functional roles. New Phytol. 182 , 314–330 (2009). Semchenko, M. et al. Fungal diversity regulates plant-soil feedbacks in temperate grassland. Sci. Adv. 4 https://doi.org/10.1126/sciadv.aau4578 (2018). Grieneisen, L. E. et al. Genes, geology and germs: gut microbiota across a primate hybrid zone are explained by site soil properties, not host species. Proc. R. Soc. B 286 , 20190431 (2019). Banerjee, S., & van der Heijden, M. G. A. Soil microbiomes and one health. Nat. Rev. Microbiol. 21 , 6–20 (2023). Additional Declarations No competing interests reported. Supplementary Files MarinangelietalTIMIDUSmicrobiotaSUPPLTables04FEB2026.xlsx MarinangelietalFinalgraphsSUPPLrevised04Feb2026.pdf Cite Share Download PDF Status: Published Journal Publication published 04 Apr, 2026 Read the published version in Scientific Reports → Version 1 posted Editorial decision: Accepted 12 Mar, 2026 Editor assigned by journal 09 Mar, 2026 Editor invited by journal 04 Mar, 2026 Submission checks completed at journal 07 Feb, 2026 First submitted to journal 04 Feb, 2026 You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-7375012","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":594538137,"identity":"7f63135d-f269-45ff-ae54-c4663f64f108","order_by":0,"name":"Lara Marinangeli","email":"","orcid":"","institution":"Fondazione Edmund Mach, S. Michele all'Adige","correspondingAuthor":false,"prefix":"","firstName":"Lara","middleName":"","lastName":"Marinangeli","suffix":""},{"id":594538138,"identity":"dc70d30c-b9e5-4b44-983f-1d3ac07422d8","order_by":1,"name":"Barbara Crestanello","email":"","orcid":"","institution":"Fondazione Edmund Mach, S. Michele all'Adige","correspondingAuthor":false,"prefix":"","firstName":"Barbara","middleName":"","lastName":"Crestanello","suffix":""},{"id":594538139,"identity":"63907dc2-fd9a-4391-a003-f2fc57a6b930","order_by":2,"name":"Nadine Praeg","email":"","orcid":"","institution":"Universität Innsbruck","correspondingAuthor":false,"prefix":"","firstName":"Nadine","middleName":"","lastName":"Praeg","suffix":""},{"id":594538140,"identity":"18a74aab-9381-4b4b-8545-fbea8a58aa5f","order_by":3,"name":"Theresa Rzehak","email":"","orcid":"","institution":"Universität Innsbruck","correspondingAuthor":false,"prefix":"","firstName":"Theresa","middleName":"","lastName":"Rzehak","suffix":""},{"id":594538141,"identity":"915c4a5c-51b1-43e5-90d4-46d75c233e57","order_by":4,"name":"Julia Seeber","email":"","orcid":"","institution":"EURAC Research","correspondingAuthor":false,"prefix":"","firstName":"Julia","middleName":"","lastName":"Seeber","suffix":""},{"id":594538142,"identity":"d6c22b00-93a9-4c1c-b63c-9d7d4ce17e85","order_by":5,"name":"Paul Illmer","email":"","orcid":"","institution":"Universität Innsbruck","correspondingAuthor":false,"prefix":"","firstName":"Paul","middleName":"","lastName":"Illmer","suffix":""},{"id":594538143,"identity":"163a96f0-9bb4-40d9-b58b-238779634311","order_by":6,"name":"Francesco Nonnis Marzano","email":"","orcid":"","institution":"University of Parma (PR)","correspondingAuthor":false,"prefix":"","firstName":"Francesco","middleName":"Nonnis","lastName":"Marzano","suffix":""},{"id":594538144,"identity":"fe3e9a0b-50e0-4c74-b4b4-1489703b1c91","order_by":7,"name":"Heidi Christine Hauffe","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA6UlEQVRIie3QMQrCMBSA4RcC7SK4Vhy8QsRBhZZeJaHgpLtjRXASXSteopNzS6BdPEAlhVYE544FF6MUwSV1dMi/JfAleQHQ6f4wAsi3KEQAOC6hsuWWBa/lL8TwCApmDVEZ8jlWEtxpiPKasbleZSXkA5JgdLPpadA9+ohXCjLdxusJhfswTDAezakYBnnU8rCMbeQsHIXFPukvaoHC12hKUlzfxJW3mI8JFW47ydCbMEkMDFSwdnJmchbCvYOcpbelwgsy5kdnFUnT66Vecmcnf6yqqXD2gcerpYI08GuF/Fag0+l0OnVPYpFbz20oYd8AAAAASUVORK5CYII=","orcid":"","institution":"Fondazione Edmund Mach, S. Michele all'Adige","correspondingAuthor":true,"prefix":"","firstName":"Heidi","middleName":"Christine","lastName":"Hauffe","suffix":""},{"id":594538145,"identity":"52f201f3-0796-411c-b920-661eb2b813ad","order_by":8,"name":"Giulio Galla","email":"","orcid":"","institution":"Fondazione Edmund Mach, S. Michele all'Adige","correspondingAuthor":false,"prefix":"","firstName":"Giulio","middleName":"","lastName":"Galla","suffix":""}],"badges":[],"createdAt":"2025-08-14 14:53:27","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-7375012/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-7375012/v1","draftVersion":[],"editorialEvents":[{"content":"https://doi.org/10.1038/s41598-026-44592-4","type":"published","date":"2026-04-04T15:58:16+00:00"}],"editorialNote":"","failedWorkflow":false,"files":[{"id":103756477,"identity":"0fa7a680-beac-470e-a5b8-cbb136f37df3","added_by":"auto","created_at":"2026-03-02 14:11:39","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":1471905,"visible":true,"origin":"","legend":"\u003cp\u003eMap of the study area showing the 12 sampling plots and their respective elevations. Bottom right: large-scale overview of the sampling area. Top right: barplot showing the faecal pellet counts for the two study species across the four elevations (m a.s.l.). LEU: \u003cem\u003eLepus europaeus\u003c/em\u003e; LTI: \u003cem\u003eLepus timidus\u003c/em\u003e.\u003c/p\u003e","description":"","filename":"MarinangelietalFigure1.png","url":"https://assets-eu.researchsquare.com/files/rs-7375012/v1/ad40294b7c5fbf2fabae120d.png"},{"id":103756519,"identity":"cd0f1d95-0835-4333-874d-f7a14be0ff24","added_by":"auto","created_at":"2026-03-02 14:11:54","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":1330935,"visible":true,"origin":"","legend":"\u003cp\u003eRelative abundances of faecal prokaryote (16S) and fungal (ITS2) taxa. (\u003cstrong\u003ea\u003c/strong\u003e) Phyla based on 16S rRNA gene data and (\u003cstrong\u003eb\u003c/strong\u003e) classes based on ITS2 data. Each bar represents a sample, with colors indicating the most abundant taxa as detailed in the legend. Samples are arranged by species (LEU: \u003cem\u003eLepus europaeus\u003c/em\u003e; LTI: \u003cem\u003eLepus timidus\u003c/em\u003e) and elevation.\u003c/p\u003e","description":"","filename":"MarinangelietalFigure2.png","url":"https://assets-eu.researchsquare.com/files/rs-7375012/v1/2d16e88b661584992680c13d.png"},{"id":103756461,"identity":"81479680-8da0-420d-8e12-a06c14c17171","added_by":"auto","created_at":"2026-03-02 14:11:31","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":2425553,"visible":true,"origin":"","legend":"\u003cp\u003eDifferential abundance analysis of prokaryote and fungal communities between the two hare species. (\u003cstrong\u003ea\u003c/strong\u003e) Prokaryote phyla and (\u003cstrong\u003eb\u003c/strong\u003e) fungal classes with significant differences in relative abundance identified using LEfSE (microeco). The x-axis represents the Linear Discriminant Analysis (LDA) scores, indicating the magnitude of difference between group means. LEU: \u003cem\u003eLepus europaeus\u003c/em\u003e; LTI: \u003cem\u003eLepus timidus\u003c/em\u003e\u003c/p\u003e","description":"","filename":"MarinangelietalFigure3.png","url":"https://assets-eu.researchsquare.com/files/rs-7375012/v1/032b2efc79bb5eed4aa00355.png"},{"id":103756395,"identity":"5c83cfad-ba79-4e64-891a-e5e525c5d3af","added_by":"auto","created_at":"2026-03-02 14:11:10","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":3423941,"visible":true,"origin":"","legend":"\u003cp\u003eNon-Metric Multidimensional Scaling (NMDS) analyses comparing prokaryote and fungal communities of\u003cem\u003e L. europaeus\u003c/em\u003eand \u003cem\u003eL. timidus\u003c/em\u003e faecal pellets. Panels (\u003cstrong\u003ea, c\u003c/strong\u003e) show prokaryote communities, while (\u003cstrong\u003eb, d\u003c/strong\u003e) show fungal communities. Dissimilarities were calculated using Bray-Curtis (\u003cstrong\u003ea, b\u003c/strong\u003e), Unifrac (\u003cstrong\u003ec\u003c/strong\u003e), and Jaccard (\u003cstrong\u003ed\u003c/strong\u003e) indices. Each dot represents a single sample (N=95). LEU: \u003cem\u003eLepus europaeus\u003c/em\u003e; LTI: \u003cem\u003eLepus timidus\u003c/em\u003e\u003c/p\u003e","description":"","filename":"MarinangelietalFigure4.png","url":"https://assets-eu.researchsquare.com/files/rs-7375012/v1/d04b6f1dfdabba5156eacd5b.png"},{"id":103756381,"identity":"acb5dceb-a17d-449c-868b-05af432c4ae5","added_by":"auto","created_at":"2026-03-02 14:11:06","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":3546044,"visible":true,"origin":"","legend":"\u003cp\u003eNon-Metric Multidimensional Scaling (NMDS) and Spearman correlation analyses comparing prokaryote and fungal community structures across different elevations. (\u003cstrong\u003ea, b\u003c/strong\u003e) NMDS clustering of prokaryote and fungalcommunities in \u003cem\u003eL. europaeus \u003c/em\u003efaecal samples\u003cem\u003e;\u003c/em\u003e (\u003cstrong\u003ec, d\u003c/strong\u003e) NMDS clustering prokaryote and fungal communities in \u003cem\u003eL. timidus \u003c/em\u003efaecal samples. Each symbol represents a single faecal sample, with different shapes indicating elevation; (\u003cstrong\u003ee\u003c/strong\u003e) Spearman correlation analysis between environmental variables (\u003cem\u003ex\u003c/em\u003e-axis) and \u003cem\u003eL. europaeus\u003c/em\u003e faecal fungal ASVs (\u003cem\u003ey\u003c/em\u003e-axis). Correlation coefficients indicate the strength and direction of the relationships. Significant correlations are indicated (* p-value \u0026lt; 0.05; ** p-value \u0026lt; 0.005; *** p-value \u0026lt; 0.0005). LEU: \u003cem\u003eLepus europaeus\u003c/em\u003e; LTI: \u003cem\u003eLepus timidus\u003c/em\u003e.\u003c/p\u003e","description":"","filename":"MarinangelietalFigure5.png","url":"https://assets-eu.researchsquare.com/files/rs-7375012/v1/a1e94a511681469afbe9cfd3.png"},{"id":104400532,"identity":"e462ef30-254c-4e82-8663-8db583b2a2ed","added_by":"auto","created_at":"2026-03-11 12:10:18","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":5713,"visible":true,"origin":"","legend":"\u003cp\u003eFunctional diversity in \u003cem\u003eLepus \u003c/em\u003espp.predicted from 16S rRNA gene sequences. (\u003cstrong\u003ea, b\u003c/strong\u003e) Shannon diversity of metaCyc pathways in faecal samples from \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e collected at 2000 and 2500 m a.s.l., grouped by species (\u003cstrong\u003ea\u003c/strong\u003e) and elevation (\u003cstrong\u003eb\u003c/strong\u003e). Significant differences were assessed using the Wilcoxon rank-sum test. **: p-value ≤ 0.01; ns: p-value \u0026gt; 0.05. (\u003cstrong\u003ec\u003c/strong\u003e) Non-Metric Multidimensional Scaling (NMDS) based on Bray-Curtis dissimilarity estimates of metaCyc pathway abundances; (\u003cstrong\u003ed\u003c/strong\u003e) Heatmap of metaCyc pathways showing significant differences in abundance (e.g. Log2FC ≥ 1; FDR adj.p-value ≤ 0.01) between sympatric \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e across elevations. Rows (i.e. metaCyc pathways; N=75) and columns (i.e. \u003cem\u003eLepus\u003c/em\u003e spp. faecal samples; N=42) are clustered using correlation distance and average linkage. Data are row-centered and scaled to unit variance. Panels (a-c) were created using ggplot2 [43] (\u003cstrong\u003ea, c\u003c/strong\u003e) and ClustVis [49] (\u003cstrong\u003ed\u003c/strong\u003e) and formatted using GIMP v2.10.18 (The GIMP Development Team, 2019). LEU: \u003cem\u003eLepus europaeus\u003c/em\u003e; LTI: \u003cem\u003eLepus timidus\u003c/em\u003e.\u003c/p\u003e","description":"","filename":"placeholderimageCopy.png","url":"https://assets-eu.researchsquare.com/files/rs-7375012/v1/db80820c0f3f6a34632211b2.png"},{"id":106343373,"identity":"217f4e7e-5cef-439d-a582-aec5f274b793","added_by":"auto","created_at":"2026-04-07 16:03:36","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":14795006,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7375012/v1/6725d50f-7201-4679-8d14-b195739a3e03.pdf"},{"id":103756511,"identity":"4c535c6b-4d8a-47c0-9922-7c52ca1ff2bd","added_by":"auto","created_at":"2026-03-02 14:11:51","extension":"xlsx","order_by":0,"title":"","display":"","copyAsset":false,"role":"supplement","size":98351,"visible":true,"origin":"","legend":"","description":"","filename":"MarinangelietalTIMIDUSmicrobiotaSUPPLTables04FEB2026.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7375012/v1/4cd8d1e2f259dfdab520cd61.xlsx"},{"id":103756583,"identity":"a58a2d93-bed2-4834-a7e7-70daca8d3327","added_by":"auto","created_at":"2026-03-02 14:12:11","extension":"pdf","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":892556,"visible":true,"origin":"","legend":"","description":"","filename":"MarinangelietalFinalgraphsSUPPLrevised04Feb2026.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7375012/v1/3a589890eff10c20b7fdbc6a.pdf"}],"financialInterests":"No competing interests reported.","formattedTitle":"Sympatric Lepus spp. in the central Italian Alps host significantly different gut microbiotas","fulltext":[{"header":"Introduction","content":"\u003cp\u003eOne of the most dramatic and rapid impacts of climate warming in the Alps is the shifting of vegetational zones to higher elevations [\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e, \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e, \u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e], accompanied by upward distributional changes of animal species as they track their preferred habitats [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e, \u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e, \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e]. An emblematic example, estimated from 30 years of hunting bag data [\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e], is the mountain hare (\u003cem\u003eLepus timidus\u003c/em\u003e), a species generally associated with boreal zones above the treeline [\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e, \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e]. However, results in [\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e] also indicated that the European brown hare (\u003cem\u003eLepus europaeus\u003c/em\u003e), a steppe-adapted species [\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e], is also shifting upward even more rapidly, which may lead to a reduction in mountain hare habitat. In addition, sympatry of the two hare species could result in their hybridization, as already reported by several authors for northern European as well as Alpine populations (e.g. Fennoscandia: [\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e]; Swiss Alps: [\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e]; western Italian Alps: [\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e] and references therein). These risk factors could threaten the long-term persistence of the mountain hare and its unique genetic pool in Alpine habitats.\u003c/p\u003e \u003cp\u003eInterspecific hybridization and range shifts (associated with dietary changes) could also lead to the loss of local adaptation, for example by altering the gut flora, which in hybrids can be intermediate between the two parental lineages (e.g. equids: [\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e]; suids: [\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e]). It is now well-recognized that an intact gut microbiota (including prokaryotes, fungi, and viruses) is essential for immune and metabolic function, as well as for development of the neuroendocrine system [\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e, \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e, \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e], impacting host health and resilience. It has been suggested that the presence of specific microorganisms in the gut may be important for host fitness and should be considered for the conservation and management of wild animal species [\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e, \u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e]. While host phylogeny has a strong influence on gut flora for many mammals, with various authors reporting species-specific compositions [\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e, \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e], extrinsic factors such as habitat preference and diet are also expected to have an impact on gut flora alpha and beta diversities [\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e, \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eSince the European brown hare is rapidly colonizing higher elevations in the Alps due to climate warming, increasing sympatry with the mountain hare is inevitable. Thus, defining gut microbiota of the two hare species is useful for understanding the potential consequences of their interaction. Here, for the first time, we describe and compare the gut microbiota composition of the two \u003cem\u003eLepus\u003c/em\u003e species in an area of sympatry in Val Mazia/Matschertal, Province of Bolzano, Italy using MiSeq amplicon sequencing of the 16S rRNA and ITS2 genes. This research aims to understand how host species and elevation influences the microbiota of these lagomorphs. Although the prokaryote gut microbiota of the European brown hare has previously been reported by Stalder et al. [\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e], no similar data exist for that of the mountain hare. In addition, the mycobiota of both species has not been studied to date, despite its potential relevance for host health [\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e]. Importantly, no studies have yet examined the gut microbiota of sympatric populations.\u003c/p\u003e"},{"header":"Materials and methods","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eStudy area and sample collection\u003c/h2\u003e \u003cp\u003eSingle fresh pellets of \u003cem\u003eLepus\u003c/em\u003e spp. were collected at least 20 m apart, from 12 sampling areas, three from each of four elevations (1000, 1500, 2000 and 2500 m a.s.l.) in July 2019 and 2020 in the Vinschgau Valley LTSER, South Tyrol, Italy (site code LTER_EU_IT_097, 46.6928\u0026deg;, 1.6157\u0026deg;; \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://deims.org/11696de6-0ab9-4c94-a06b-7ce40f56c964\u003c/span\u003e\u003cspan address=\"https://deims.org/11696de6-0ab9-4c94-a06b-7ce40f56c964\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e; Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e) as part of the EUREGIO project MICROVALU. Freshly deposited faecal pellets were identified by their shiny (mucous smeared) external surface. Fieldworkers wore N95 masks and used sterile gloves and tweezers to collect the pellets, which were placed in certified sterile DNA/DNAse-free 15 ml tubes (Sarstedt, Germany) and stored at -20\u0026deg;C for up to 24h before being transferred to the Animal, Environmental and Antique DNA Platform at the Fondazione E. Mach (FEM), where they were archived at -80\u0026deg;C until further processing. Sample details are provided in the Supplementary Tab. S1.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eMeteorological data were collected on site using Campbell Scientific CR1000/CR1000X loggers and provided by Zandonai et al. [\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e], including daily average temperature in the three weeks prior to sampling and sum of precipitation of the entire vegetation period (May 1st to Sept 30th ). Temperature was recorded with three loggers per elevation, while precipitation was recorded by one meteorological station at each elevation. Species richness was estimated using a square plot with a grid of 100 cm\u003csup\u003e2\u003c/sup\u003e (10x10 cm) units to calculate the number of plant species in each elevational zone, while plant richness was obtained from Hilpold et al. [\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e]. Botanical data of the 2500 m a.s.l. sites was then added to the previously published elevation surveys.\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eDNA extraction\u003c/h3\u003e\n\u003cp\u003eUnder a sterile Class II biological safety cabinet (BSL2, Telstar, Spain) and using sterile DNA/DNase-free equipment and consumables, 50 mg of faecal material were taken from the centre of each pellet to avoid the outer layer that had potentially been in contact with soil. Total DNA was extracted from each 50 mg subsample of a single faecal pellet with the NucleoSpin Soil mini kit (Macherey-Nagel, Germany) using the lysis buffer SL1 in combination with 50 \u0026micro;l of Buffer SX, as suggested by Praeg, Pauli, and Illmer [\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e]. Each sample was co-extracted with a 0.15 dose of ZymoBIOMICS\u0026trade; Spike-in Control I (High Microbial Load; EuroClone, Milan, Italy), before the initial lysis step following Galla et al. [\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e], unless otherwise stated (Tab. S1). This mock community was added as an in situ positive control to validate the accuracy and reliability of sequencing and bioinformatics pipeline. Negative DNA extraction controls were included up to the sequencing step to monitor environmental and reagent contaminations. The purity and quantity of all DNA extracts were assessed by checking the UV/VIS spectra of each extract with a Spark multimode microplate reader (Tecan, Switzerland). Following quantification, all DNA extracts were diluted to a final concentration of 3 ng/\u0026micro;l using nuclease-free water and subsequently used for species identification and amplicon sequencing.\u003c/p\u003e\n\u003ch3\u003eHost species identification\u003c/h3\u003e\n\u003cp\u003eTo determine the origin of each pellet (\u003cem\u003eL. europaeus\u003c/em\u003e or \u003cem\u003eL. timidus\u003c/em\u003e), a\u0026thinsp;~\u0026thinsp;500 bp-fragment of the mitochondrial DNA D-loop region was amplified using GoTaq DNA polymerase by Promega with primers (H) GTTGCTGGTTTCACGGAGGTAG and (L) TCCTACCATCAGCACCCAAAGC\u0026rsquo; previously designed by Wilkinson \u0026amp; Chapman [\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e]. Cycling conditions consisted of an initial denaturation step at 94\u0026deg;C for 3 min followed by 40 cycles of 94\u0026deg;C for 30 s, 61\u0026deg;C for 30 s and 72\u0026deg;C for 45 s, with a final extension step at 72\u0026deg;C for 5 min. Negative PCR controls (PCR-grade water instead of DNA template) were included to check for contamination. Amplified DNA was purified using ExoProStarTM 1-Step (Cytiva, USA) and sequenced at the FEM Sequencing and Genotyping Platform with both primers H and L using the ABI Prism BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems). Sequencing reaction products were further purified with the Dynabeads Sequencing Clean-Up Kit (ThermoFisher Scientific) and subsequently analyzed on an ABI 3130XL DNA sequencer (Applied Biosystems, Foster City, CA, USA), 96 capillaries (50 cm) and POP07 Performance Optimized Polymer (Applied Biosystems). The sequences were processed with Sequencher software (Gene Codes Corporation-USA) to delete any automatic base assignment errors. BLASTn was used to identify sequences with \u0026ge;\u0026thinsp;99% identity as defined using the Standard Nucleotide method. Moreover, D-loop sequences were trimmed to a 334 bp region shared by all sequences as well as 147 D-loop sequences deposited in NCBI by Melo-Ferreira et al. [\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e], Pecchioli et al. (unpublished PhD thesis) [\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e] and Zachos et al. [\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e] and aligned with MEGA v. 10.0.5 (option Muscle; [\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e]). A Neighbor-Joining phylogenetic tree was constructed in MEGA, using the Phylogeny option, based on pairwise genetic distance.\u003c/p\u003e\n\u003ch3\u003eCharacterization of V3-V4 16S rRNA gene and ITS2 faecal communities\u003c/h3\u003e\n\u003cp\u003eThe amplification of the V3-V4 regions of the prokaryote 16S rRNA gene was performed using the KAPA HiFi HS ReadyMix (Roche) in a 25 \u0026micro;l reaction volume containing 1X KAPA HiFi HS ReadyMix Buffer, 0.3 \u0026micro;M each of primers 341F_ILL (CCTACGGGNGGCWGCAG) and 805R-2_ILL (GACTACNVGGGTWTCTAATCC) modified with Illumina overhang adapters [\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e, \u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e] (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://support.illumina.com/documents/documentation/chemistry_documentation/16s\u003c/span\u003e\u003cspan address=\"https://support.illumina.com/documents/documentation/chemistry_documentation/16s\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) and 9 ng of DNA. The thermal profile for 16S V3-V4 amplification reactions was: 3 min at 95\u0026deg;C, followed by 31 cycles of 30 sec at 95\u0026deg;C, 30 sec at 55\u0026deg;C, 30 sec at 72\u0026deg;C, and a single final extension step of 5 minutes at 72\u0026deg;C. The amplification of the fungal rRNA Internal Transcribed Spacer ITS2 was performed in a 25\u0026micro;l reaction volume containing 1X FastStart High Fidelity Reaction Buffer (Roche Applied Science), 0.4 \u0026micro;M each of primers gITS7 (Illumina_forward_overhang-GTGARTCATCGARTCTTTG) and ITS4 (Illumina_reverse_overhang-TCCTCCGCTTATTGATATGC) [\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e, \u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e], 200 \u0026micro;M each dNTPs, 9 ng of DNA, and 1.5 U of FastStart High Fidelity Enzyme Blend (Roche Applied Science, Germany). The thermal profile for ITS2 amplification reactions was 3 min at 95\u0026deg;C, followed by 35 cycles of 45 sec at 95\u0026deg;C, 45 sec at 50\u0026deg;C, 45 sec at 72\u0026deg;C, and a single final extension step of 5 minutes at 72\u0026deg;C. All amplifications were performed on a Veriti 96-Well Fast Thermal Cycler (Applied Biosystems, USA). Non-template controls were included in each amplification reaction. Amplicons were visualized by high-resolution capillary electrophoresis using the QIAxcel Advanced System (Qiagen, Hilden, Germany). Quantification of individual amplicon libraries, normalization at equimolar concentrations, pooling of indexed libraries and high throughput sequencing by Illumina technology were performed at the FEM Sequencing and Genotyping Platform. Libraries were sequenced on an Illumina MiSeq platform using Standard Flow Cells (PE300), targeting a minimum sequencing depth of 30,000 reads per amplicon.\u003c/p\u003e\n\u003ch3\u003eAmplicon sequencing data analysis\u003c/h3\u003e\n\u003cp\u003eBioinformatic analyses were performed in R v. 4.3.1. DADA2 v. 3.14 [\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e] was used for sequence quality inspection, primer trimming, chimera removal and ASVs taxonomy assignment following the standard operating procedure. Taxonomic classification of 16S rRNA gene and ITS2 ASVs was performed using Silva v. 138.1 [\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e] and UNITE v. 9.0 [\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e] reference databases, respectively. Multiple sequence alignment of 16S rRNA gene ASVs was done with msa v. 1.32.0 [\u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e] using ClustalW and default parameters. With the Phangorn package v. 2.11.1 [\u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e], a phylogenetic tree was constructed with 16S rRNA gene ASVs, which was then used for further analysis.\u003c/p\u003e \u003cp\u003eSequencing results of negative controls from both DNA extraction (N\u0026thinsp;=\u0026thinsp;12) and PCR amplification (N\u0026thinsp;=\u0026thinsp;13) were used to identify potential contaminants using decontam v. 1.20.0 [\u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e45\u003c/span\u003e] with the prevalence method. All mock community-related ASVs were removed from the dataset before performing subsequent analyses. Statistical analyses of 16S rRNA gene and ITS2 ASVs were performed using phyloseq v. 1.50.0 [\u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e] and microeco v. 4.3.1 [\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e]. Images were generated using ggplot2 v. 3.3.5 [\u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e48\u003c/span\u003e], using RColorBrewer v. 1.1-2 [\u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e49\u003c/span\u003e] and viridis v. 0.6.5 [\u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e50\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eThe impact of species and elevation on taxonomic composition (alpha and beta diversity) of faecal microbiota were investigated by focusing on elevations where the two species live in sympatry (Tab. S1). \u003cem\u003eL. timidus\u003c/em\u003e and \u003cem\u003eL. europaeus\u003c/em\u003e datasets were then analysed separately to assess intraspecific variation in faecal microbiota across all elevations where they were sampled.\u003c/p\u003e \u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003ch2\u003eTaxonomic profiling of 16S rRNA and ITS2 regions\u003c/h2\u003e \u003cp\u003eTaxonomic classification of 16S rRNA gene and ITS2 ASVs was estimated on non-rarefied data and expressed as relative abundances. Visualizations were generated using microeco, focusing on the 10 most abundant phyla and 16 most abundant classes and families. The relative abundance of prokaryote and fungal taxa was compared between the two species with LEfSE [\u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e51\u003c/span\u003e].\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eAlpha and beta diversity estimation\u003c/h3\u003e\n\u003cp\u003eTo avoid bias in richness estimates due to different sequencing depths, alpha diversity indices Chao1, inverse Simpson, Shannon and Faith\u0026rsquo;s Phylogenetic Diversity (PD; this latter index for prokaryotes only) were estimated using 16S rRNA gene and ITS2 datasets rarefied to 10,564 and 10,804 ASVs per sample, respectively. Alpha diversity estimates were compared between \u003cem\u003eLepus\u003c/em\u003e species and elevations using Wilcoxon rank-sum tests, with microeco. Venn diagrams representing the number of ASVs shared between \u003cem\u003eL. timidus\u003c/em\u003e and \u003cem\u003eL. europaeus\u003c/em\u003e were generated using the rarefied data, microeco. Beta diversity estimates were computed on non-rarefied ASV abundances normalized as relative abundances using Bray-Curtis, unweighted Unifrac (for 16S rRNA gene only) and Jaccard (for ITS2 only) indices. Ordination based on distance/dissimilarity matrices was done using non-metric multidimensional scaling (NMDS). Multivariate homogeneity of group dispersions, was tested using betadisper with Vegan v. 2.6.4 [\u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e52\u003c/span\u003e]. To compare beta-diversity estimates between groups, permutational multivariate ANOVA (PERMANOVA), and pairwise-PERMANOVA were performed with 999 permutations using the function adonis2 implemented in microeco. The same package was used for the analysis of similarities (ANOSIM). Pairwise-PERMANOVA p-values were adjusted using the false discovery rate (FDR) method to control for multiple testing.\u003c/p\u003e\n\u003ch3\u003eImpact of environmental variables on beta diversity\u003c/h3\u003e\n\u003cp\u003eSignificant differences in daily temperature, pre-sampling total precipitation and plant species richness between elevations were inspected using a Wilcox rank-sum test (data not shown). Autocorrelation among environmental variables was assessed using Spearman's Rank-Order Correlation using microeco. For both 16S rRNA gene and ITS2 datasets, Mantel tests employing Spearman rank order and Pearson correlation coefficients were used to examine associations between environmental parameters and beta diversity estimates. Spearman rank-order correlations were calculated to evaluate relationships between environmental parameters and the abundance of each 16S rRNA and ITS2 ASVs, also using microeco.\u003c/p\u003e \u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003eCharacterization of prokaryotic functional diversity\u003c/h2\u003e \u003cp\u003eProkaryote functional diversity was inspected using the PICRUSt2 pipeline v. 2.5.3 [\u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e53\u003c/span\u003e] implemented on Galaxy (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://usegalaxy.eu/\u003c/span\u003e\u003cspan address=\"https://usegalaxy.eu/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003e)\u003c/span\u003e with default settings. Estimated MetaCyc pathway abundances were normalized to relative abundances. Alpha diversity (Shannon index) and beta diversity (Bray-Curtis dissimilarity) were estimated from normalized MetaCyc pathway abundances using microeco. Differences in MetaCyc pathway abundances between \u003cem\u003eLepus\u003c/em\u003e species and elevations was tested using Wilcoxon rank-sum text with FDR p-value adjustment. Differential pathways were defined by a Log\u003csub\u003e2\u003c/sub\u003eFC cutoff of \u0026ge;\u0026thinsp;1 and an adjusted p-value of \u0026le;\u0026thinsp;0.01. Plots were generated with ggplot2. A heatmap clustering samples and pathways based on MetaCyc pathways abundances was generated using ClustVis [\u003cspan citationid=\"CR54\" class=\"CitationRef\"\u003e54\u003c/span\u003e] with default parameters and further formatted using GIMP v. 2.10.18 [\u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e55\u003c/span\u003e]\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003e)\u003c/span\u003e.\u003c/p\u003e \u003c/div\u003e"},{"header":"Results","content":"\u003cp\u003e \u003cb\u003eLepus\u003c/b\u003e \u003cb\u003especies identification\u003c/b\u003e\u003c/p\u003e \u003cp\u003eOf 108 total fresh faecal pellets, 95 provided reliable mtDNA D-loop results, identifying 72 \u003cem\u003eL. europaeus\u003c/em\u003e and 23 \u003cem\u003eL. timidus\u003c/em\u003e samples. As expected, \u003cem\u003eL. timidus\u003c/em\u003e was restricted to higher elevations, with nine and 14 samples found at 2000 m and 2500 m a.s.l., respectively. In contrast, \u003cem\u003eL. europaeus\u003c/em\u003e ranged across the study area, with 27, 26, 16 and 3 samples collected at 1000, 1500, 2000, and 2500 m a.s.l., respectively (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e).\u003c/p\u003e \u003cp\u003e \u003cb\u003eTaxonomic composition of\u003c/b\u003e \u003cb\u003eL. timidus\u003c/b\u003e \u003cb\u003eand\u003c/b\u003e \u003cb\u003eL. europaeus\u003c/b\u003e \u003cb\u003efaecal microbiota\u003c/b\u003e\u003c/p\u003e \u003cp\u003eThe profiling of \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e faecal microbiota identified 6145 ASVs for the prokaryote 16S rRNA gene and 2219 ASVs for the fungal ITS2 region. Of these, 22 and 16 ASVs were identified by decontam as potential contaminants and removed from the 16S rRNA gene and ITS2 datasets, respectively. Additionally, 24 ASVs were classified as belonging to the two taxa composing the mock community: \u003cem\u003eAllobacillus\u003c/em\u003e spp. (7 ASVs) and \u003cem\u003eImtechella\u003c/em\u003e spp. (17 ASVs); all of these ASVs were also removed from the dataset before statistical analyses. The overall median relative abundance of these mock community-related ASVs was 3.2% (ranging from 0.2% to 15%) and was comparable between \u003cem\u003eLepus\u003c/em\u003e species and elevations (Kruskal-Wallis, p-value: 0.3487). Following rarefaction, 16S rRNA gene and ITS2 datasets accounted for 4877 and 1514 ASVs, respectively (Fig. \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003ea-b). Rarefaction resulted in the loss of seven 16S rRNA gene (5 \u003cem\u003eL. europaeus\u003c/em\u003e and 2 \u003cem\u003eL. timidus\u003c/em\u003e) and five ITS2 libraries (three \u003cem\u003eL. europaeus\u003c/em\u003e and two \u003cem\u003eL. timidus\u003c/em\u003e).\u003c/p\u003e \u003cp\u003eBoth species were characterized by private ASVs accounting for approximately 82% of 16S rRNA and 12% ITS2 sequence reads. In contrast, only about 6% of prokaryote and 9% of fungal ASVs were shared by faecal samples of both species. Private \u003cem\u003eL. europaeus\u003c/em\u003e prokaryote ASVs prevailed in the gut community which consisted of about 70% of sequenced reads. In contrast, dominant fungal ASVs were shared, accounting for approximately 84% of generated sequence reads (Fig. \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003ea-b). Restricting the analysis to elevations where the two species were sympatric produced results consistent with those above, with approximately 45% of prokaryote ASVs found only in \u003cem\u003eL. europaeus\u003c/em\u003e samples, 35% of prokaryote ASVs uniquely found in \u003cem\u003eL. timidus\u003c/em\u003e samples (Fig. \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003ec) and only a small fraction of prokaryote ASVs found in both species with relatively high abundance (e.g. 10 ASVs detected in both species and elevations accounted for 10% of sequence reads). Regarding fungi, most ASVs were shared by both species at all elevations (Fig. \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003ed).\u003c/p\u003e \u003cp\u003eThe taxonomic classification of prokaryote phyla and fungal classes found in \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e faecal pellets is shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e. At all elevations, the faecal microbiota of \u003cem\u003eL. europaeus\u003c/em\u003e was dominated by the prokaryote phyla Firmicutes (median: 63.3%, range: 0.28% \u0026minus;\u0026thinsp;80.8%) and Bacteroidota (median: 16.4%, range: undetected\u0026thinsp;\u0026minus;\u0026thinsp;33.0%) followed, but with much lower percentages, by Actinobacteriota (median: 2.0%, range: undetected\u0026thinsp;\u0026minus;\u0026thinsp;70.3%), Spirochaetota (median: 1.2%, range: undetected\u0026thinsp;\u0026minus;\u0026thinsp;17.1%) and Patescibacteria (median: 1.0%, range: undetected\u0026thinsp;\u0026minus;\u0026thinsp;30.5%). In addition, about 15% of the 72 samples collected across all four elevations were characterised by a notably high abundance of Proteobacteria (median: 0.7%, range: 0.1% \u0026minus;\u0026thinsp;88.0%) (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea). In contrast, \u003cem\u003eL. timidus\u003c/em\u003e faecal microbiota was dominated by Proteobacteria (median: 62.9%, range: 4.6% \u0026minus;\u0026thinsp;78.9%), Bacteroiodota (median: 22.3%, range: undetected\u0026thinsp;\u0026minus;\u0026thinsp;29.2%) and, to a lesser extent, Actinobacteriota (median: 5.7%, range: 0.2% \u0026minus;\u0026thinsp;25.8%) and Acidobacteriota (median: 1.2%, range: undetected\u0026thinsp;\u0026minus;\u0026thinsp;4.2%) (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea). Firmicutes was detected in \u003cem\u003eL. timidus\u003c/em\u003e faecal samples as well, but in lower percentages (median: 0.3%, range: 0.1% \u0026minus;\u0026thinsp;26.2%).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eThe faecal mycobiota of both species was characterized by a high relative abundance of Dothideomycetes (\u003cem\u003eL. europaeus\u003c/em\u003e median: 51.8%, range: 6.3\u0026ndash;87.9%; \u003cem\u003eL. timidus\u003c/em\u003e median: 43.6%, range: 2.8\u0026ndash;75.1%), Leotiomicetes (\u003cem\u003eL. europaeus\u003c/em\u003e median: 15.5%, range: 1.0% \u0026minus;\u0026thinsp;74.9%; \u003cem\u003eL. timidus\u003c/em\u003e median: 22.1%, range: 0.4% \u0026minus;\u0026thinsp;81.3%) and Sordariomycetes (\u003cem\u003eL. europaeus\u003c/em\u003e median: 6.0%, range: undetected\u0026thinsp;\u0026minus;\u0026thinsp;49.1%; \u003cem\u003eL. timidus\u003c/em\u003e median: 0.9%, range: 0.1\u0026ndash;20.3%). \u003cem\u003eL. europaeus\u003c/em\u003e faecal samples also had non negligible occurrences of Pezizomycetes (median: 1.2%, range: undetected\u0026thinsp;\u0026minus;\u0026thinsp;66.6%) and Microbotryomycetes (median: 0.1%, range: undetected\u0026thinsp;\u0026minus;\u0026thinsp;40.9%) (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eb), whereas \u003cem\u003eL. timidus\u003c/em\u003e faecal samples were notable for the abundance of Tremellomycetes (median: 14.9%, range: 0.1% \u0026minus;\u0026thinsp;53.4%), this latter class being particularly abundant in samples collected at 2500 m a.s.l. (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eb).\u003c/p\u003e \u003cp\u003eThe taxonomic composition of faecal microbiota was significantly different between the two study species where they were found to be sympatric (2000 and 2500 m a.s.l.). As shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e and Tab. S2, \u003cem\u003eL. europaeus\u003c/em\u003e faecal microbiota was significantly enriched with the following prokaryote phyla: Firmicutes (LEfSE; LDA\u0026thinsp;=\u0026thinsp;5.47, adjusted p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.0001), Spirochaetota (LEfSE; LDA\u0026thinsp;=\u0026thinsp;4.11, adjusted p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.0001), Verrucomicrobiota (LEfSE; LDA\u0026thinsp;=\u0026thinsp;3.57, adjusted p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.0001). Instead, \u003cem\u003eL. timidus\u003c/em\u003e faecal microbiota had significantly higher abundances of Proteobacteria (LEfSE; LDA\u0026thinsp;=\u0026thinsp;5.42, adjusted p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.0001), Actinobacteriota (LEfSE; LDA\u0026thinsp;=\u0026thinsp;4.05, adjusted p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.0001), Acidobacteriota (LEfSE; LDA\u0026thinsp;=\u0026thinsp;3.93, adjusted p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.0001) and Myxococcota (LEfSE; LDA\u0026thinsp;=\u0026thinsp;3.38, adjusted p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.001), among others (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ea, Fig. \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e, Tab. S2).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eWe also found differences in relative abundance of fungal classes between the two host species at 2000 and 2500 m a.s.l., mainly related to taxa belonging to Ustilaginomycetes (phylum Basidiomycota; LEfSE; LDA\u0026thinsp;=\u0026thinsp;4.00, adjusted p-value\u0026thinsp;=\u0026thinsp;2.57E-3) and Pezizomycetes (phylum Ascomycota; LEfSE; LDA\u0026thinsp;=\u0026thinsp;3.86, adjusted p-value\u0026thinsp;=\u0026thinsp;4.13E-3), that were all more abundant in \u003cem\u003eL. europaeus\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e, Tab. S2). Instead, the two Basidiomycota fungal classes Tremellomycetes (LEfSE; LDA\u0026thinsp;=\u0026thinsp;4.93, adjusted p-value\u0026thinsp;=\u0026thinsp;9.34E-4) and Cystobasidiomycetes (LEfSE; LDA\u0026thinsp;=\u0026thinsp;3.67, adjusted p-value\u0026thinsp;=\u0026thinsp;5.57E-3) had a higher relative abundance in \u003cem\u003eL. timidus\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eb, Fig. \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e, Tab. S2).\u003c/p\u003e \u003cp\u003e \u003cb\u003eInterspecific variation in diversity of\u003c/b\u003e \u003cb\u003eL. timidus\u003c/b\u003e \u003cb\u003eand\u003c/b\u003e \u003cb\u003eL. europaeus\u003c/b\u003e \u003cb\u003efaecal microbiota\u003c/b\u003e\u003c/p\u003e \u003cp\u003eWilcoxon rank-sum tests indicated no significant differences in faecal microbiota richness (Chao1), diversity (Shannon, Inverse Simpson), or Faith\u0026rsquo;s phylogenetic diversity (PD) between the two \u003cem\u003eLepus\u003c/em\u003e spp. species (Fig. S3), regardless of whether all four elevations were included or only the two where both species were found in sympatry (i.e. 2000m, 2500 m a.s.l.). However, NMDS based on beta diversity estimates highlighted a clear separation between the microbial communities of \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea, c). Consistently, the PERMANOVA analysis performed on samples collected from both species at 2000 and 2500 m a.s.l. which indicated significant differences in prokaryote community composition (PERMANOVA: unweighted Unifrac: R\u003csup\u003e2\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.20, p-value\u0026thinsp;=\u0026thinsp;0.001; Bray-Curtis: R\u003csup\u003e2\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.27, p-value\u0026thinsp;=\u0026thinsp;0.001; Table\u0026nbsp;\u003cspan refid=\"Tab1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). These findings were further corroborated by the lack of differences in group dispersions across species and elevations (betadisper; p-value\u0026thinsp;\u0026gt;\u0026thinsp;0.05 in all cases except \u003cem\u003eL. timidus\u003c/em\u003e at 2000 m a.s.l. vs \u003cem\u003eL. timidus\u003c/em\u003e at 2500 m a.s.l.). The clear differentiation in faecal prokaryotic communities between \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e at 2000 and 2500 m a.s.l. was further confirmed by ANOSIM (Bray-Curtis: R\u0026thinsp;=\u0026thinsp;0.886, p-value\u0026thinsp;=\u0026thinsp;0.001; unweighted Unifrac R\u0026thinsp;=\u0026thinsp;0.675, p-value\u0026thinsp;=\u0026thinsp;0.001).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003cdiv class=\"gridtable\"\u003e\u003ctable float=\"Yes\" id=\"Tab1\" border=\"1\"\u003e \u003ccaption language=\"En\"\u003e \u003cdiv class=\"CaptionNumber\"\u003eTable 1\u003c/div\u003e \u003cdiv class=\"CaptionContent\"\u003e \u003cp\u003eResults of PERMANOVA analysis comparing the faecal microbiota of \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e at sympatric elevations. Microbial diversity was estimated using Bray-Curtis, Unweighted Unifrac (for 16S rRNA gene only) and Jaccard index (ITS2 locus only). For each microbial community and diversity index, statistical tests were performed by considering species, elevation and their interactions.\u003c/p\u003e \u003c/div\u003e \u003c/caption\u003e \u003ccolgroup cols=\"7\"\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c4\" colnum=\"4\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c5\" colnum=\"5\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c6\" colnum=\"6\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c7\" colnum=\"7\"\u003e\u003c/div\u003e \u003cthead\u003e \u003ctr\u003e \u003cth align=\"left\" colname=\"c1\"\u003e \u003cp\u003eMicrobial community (marker)\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c2\"\u003e \u003cp\u003eDiversity Index\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c3\"\u003e \u003cp\u003eVariable\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c4\"\u003e \u003cp\u003eR\u003csup\u003e2\u003c/sup\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c5\"\u003e \u003cp\u003eF\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c6\"\u003e \u003cp\u003ep-value\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c7\"\u003e \u003cp\u003eSig.\u003c/p\u003e \u003c/th\u003e \u003c/tr\u003e \u003c/thead\u003e \u003ctbody\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\" morerows=\"7\" rowspan=\"8\"\u003e \u003cp\u003e\u003cb\u003eProkaryotes\u003c/b\u003e\u003c/p\u003e \u003cp\u003e\u003cb\u003e(16S rRNA gene)\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\" morerows=\"3\" rowspan=\"4\"\u003e \u003cp\u003e\u003cb\u003eBray-Curtis\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eSpecies\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e0.273\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e15.207\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e0.001\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003e***\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eElevation\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e0.022\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e1.252\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e0.203\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003ens\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eSpecies:Elevation\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e0.024\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e1.336\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e0.171\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003ens\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eResidual\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e0.681\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e\u0026nbsp;\u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e\u0026nbsp;\u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003ens\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c2\" morerows=\"3\" rowspan=\"4\"\u003e \u003cp\u003e\u003cb\u003eUnweighted Unifrac\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eSpecies\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e0.200\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e9.903\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e0.001\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003e***\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eElevation\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e0.017\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e0.855\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e0.528\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003ens\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eSpecies:Elevation\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e0.016\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e0.774\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e0.655\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003ens\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eResidual\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e0.767\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e\u0026nbsp;\u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e\u0026nbsp;\u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003ens\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\" morerows=\"7\" rowspan=\"8\"\u003e \u003cp\u003e\u003cb\u003eFungi\u003c/b\u003e\u003c/p\u003e \u003cp\u003e\u003cb\u003e(ITS2 locus)\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\" morerows=\"3\" rowspan=\"4\"\u003e \u003cp\u003e\u003cb\u003eBray-Curtis\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eSpecies\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e0.240\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e13.659\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e0.001\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003e***\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eElevation\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e0.052\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e2.967\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e0.007\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003e**\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eSpecies:Elevation\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e0.041\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e2.332\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e0.022\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003e*\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eResidual\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e0.667\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e\u0026nbsp;\u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e\u0026nbsp;\u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003ens\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c2\" morerows=\"3\" rowspan=\"4\"\u003e \u003cp\u003e\u003cb\u003eJaccard\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eSpecies\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e0.133\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e6.266\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e0.001\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003e***\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eElevation\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e0.033\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e1.534\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e0.046\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003e*\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eSpecies:Elevation\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e0.029\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e1.364\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e0.085\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003ens\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eResidual\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e0.806\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e\u0026nbsp;\u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e\u0026nbsp;\u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003ens\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003c/tbody\u003e \u003c/colgroup\u003e \u003ctfoot\u003e \u003ctr\u003e\u003ctd colspan=\"7\"\u003ep-value\u0026thinsp;\u0026le;\u0026thinsp;0.001: ***; p-value\u0026thinsp;\u0026le;\u0026thinsp;0.01: **; p-value\u0026thinsp;\u0026le;\u0026thinsp;0.05: *; ns: not significant.\u003c/td\u003e\u003c/tr\u003e \u003c/tfoot\u003e \u003c/table\u003e\u003c/div\u003e \u003c/p\u003e \u003cp\u003eFor the fungal community, NMDS clustering of \u003cem\u003eL. timidus\u003c/em\u003e and \u003cem\u003eL. europaeus\u003c/em\u003e faecal samples collected at sympatric elevations revealed distinct clusters associated with species (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eb, d). This was supported by PERMANOVA analyses, which highlighted significant differences in the fungal community composition between \u003cem\u003eL. timidus\u003c/em\u003e and \u003cem\u003eL. europaeus\u003c/em\u003e (PERMANOVA: Jaccard: R\u003csup\u003e2\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.133, p-value\u0026thinsp;=\u0026thinsp;0.001; Bray-Curtis R\u003csup\u003e2\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.240, p-value\u0026thinsp;=\u0026thinsp;0.001), as well as across the two considered elevations (PERMANOVA: Jaccard: R\u003csup\u003e2\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.033, p-value\u0026thinsp;=\u0026thinsp;0.046; Bray-Curtis R\u003csup\u003e2\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.052, p-value\u0026thinsp;=\u0026thinsp;0.007). A weaker but still significant species-by-elevation interaction was observed (Jaccard: R\u003csup\u003e2\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.029, p-value\u0026thinsp;=\u0026thinsp;0.085; Bray-Curtis R\u003csup\u003e2\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.041, p-value\u0026thinsp;=\u0026thinsp;0.022; Table\u0026nbsp;\u003cspan refid=\"Tab1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). Additional support for species-specific differences in mycobiota at sympatric elevations was provided by ANOSIM (Bray-Curtis: R\u0026thinsp;=\u0026thinsp;0.438, p-value\u0026thinsp;=\u0026thinsp;0.001; Jaccard R\u0026thinsp;=\u0026thinsp;0.517, p-value\u0026thinsp;=\u0026thinsp;0.001).\u003c/p\u003e \u003cdiv id=\"Sec13\" class=\"Section2\"\u003e \u003ch2\u003eIntraspecific variation in diversity along an elevational gradient\u003c/h2\u003e \u003cp\u003eNMDS clustering and PERMANOVA analysis of \u003cem\u003eL. europaeus\u003c/em\u003e samples across the entire elevational gradient (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ea, b; Tab. S3) provided limited evidence for an association between elevation and prokaryote community composition (PERMANOVA: Bray-Curtis: R\u003csup\u003e2\u003c/sup\u003e: 0.064, p-value\u0026thinsp;=\u0026thinsp;0.011; Unifrac: R\u003csup\u003e2\u003c/sup\u003e: 0.053, p-value\u0026thinsp;=\u0026thinsp;0.082). This result was also confirmed by the ANOSIM results which showed no significant differences across elevations (Bray-Curtis: R\u0026thinsp;=\u0026thinsp;0.046, p-value\u0026thinsp;=\u0026thinsp;0.075; unweighted Unifrac R\u0026thinsp;=\u0026thinsp;0.032, p-value\u0026thinsp;=\u0026thinsp;0.129). Conversely, the same ordination and statistical analysis of \u003cem\u003eL. europaeus\u003c/em\u003e fungal communities revealed a clear association between elevation and composition of faecal mycobiota (PERMANOVA: Bray-Curtis R\u003csup\u003e2\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.153, p-value\u0026thinsp;=\u0026thinsp;0.001. Jaccard R\u003csup\u003e2\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.098, p-value\u0026thinsp;=\u0026thinsp;0.001; ANOSIM; Bray-Curtis: R\u0026thinsp;=\u0026thinsp;0.244, p-value\u0026thinsp;=\u0026thinsp;0.075; Jaccard: R\u0026thinsp;=\u0026thinsp;0.3075, p-value\u0026thinsp;=\u0026thinsp;0.129) (Tab. S3, Fig. S4b). Similarly, no association with elevation was detected for \u003cem\u003eL. timidus\u003c/em\u003e prokaryote beta diversity estimates (PERMANOVA and ANOSIM: p-value\u0026thinsp;\u0026gt;\u0026thinsp;0.05; Tab. S3; Fig. S4c). In contrast, mycobiota composition showed a significant association with elevation (PERMANOVA: Bray-Curtis R\u003csup\u003e2\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.118, p-value\u0026thinsp;=\u0026thinsp;0.018. Jaccard R\u003csup\u003e2\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.079, p-value\u0026thinsp;=\u0026thinsp;0.014; ANOSIM; Bray-Curtis: R\u0026thinsp;=\u0026thinsp;0.339, p-value\u0026thinsp;=\u0026thinsp;0.001; Jaccard: R\u0026thinsp;=\u0026thinsp;0.360, p-value\u0026thinsp;=\u0026thinsp;0.005; Tab. S3, Fig. S4d). Consistently, pairwise PERMANOVAs performed separately for \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e revealed significant differences in fungal beta diversity between faecal samples collected at all elevation pairs, including those from neighbouring elevations (pairwise PERMANOVA: adjusted p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.05 for both species in all pairwise comparisons; Tab. S4). Conversely, no significant differences were observed in prokaryote communities across elevations for either species, including those between the most distant sites (pairwise PERMANOVA: adjusted p-value\u0026thinsp;\u0026gt;\u0026thinsp;0.05 for both species in all pairwise comparisons; Tab. S4).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003cb\u003eInfluence of environmental factors on\u003c/b\u003e \u003cb\u003eLepus\u003c/b\u003e \u003cb\u003espp. microbial communities\u003c/b\u003e\u003c/p\u003e \u003cp\u003eNo significant correlations were detected between prokaryote diversity and temperature, pre-sampling precipitation, or plant richness for either \u003cem\u003eL. europaeus\u003c/em\u003e or \u003cem\u003eL. timidus\u003c/em\u003e (Mantel test: Spearman's ρ adjusted p-values\u0026thinsp;\u0026gt;\u0026thinsp;0.01; Tab. S5). Instead, fungal diversity in faecal samples of \u003cem\u003eL. europaeus\u003c/em\u003e, but not \u003cem\u003eL. timidus\u003c/em\u003e, was significantly associated with all three environmental variables: temperature (Mantel test; Bray-Curtis: Spearman's ρ\u0026thinsp;=\u0026thinsp;0.288, adjusted p-values\u0026thinsp;=\u0026thinsp;0.002; Jaccard: Spearman's ρ\u0026thinsp;=\u0026thinsp;0.262, adjusted p-values\u0026thinsp;=\u0026thinsp;0.001), pre-sampling total precipitation (Mantel test; Bray-Curtis: Spearman's ρ\u0026thinsp;=\u0026thinsp;0.118, adjusted p-values\u0026thinsp;=\u0026thinsp;0.006; Jaccard: Spearman's ρ\u0026thinsp;=\u0026thinsp;0.255, adjusted p-values\u0026thinsp;=\u0026thinsp;0.001) and plant richness (Mantel test; Bray-Curtis: Spearman's ρ\u0026thinsp;=\u0026thinsp;0.123, adjusted p-values\u0026thinsp;=\u0026thinsp;0.002; Jaccard: Spearman's ρ\u0026thinsp;=\u0026thinsp;0.279, adjusted p-values\u0026thinsp;=\u0026thinsp;0.001) (Tab. S5). Furthermore, several fungal ASVs were found to be correlated with pre-sampling precipitation (four ASVs, positively correlated), temperature (15 ASVs positively correlated; seven ASVs negatively correlated), and plant richness (five ASVs positively correlated; 17 ASVs negatively correlated) (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ee).\u003c/p\u003e \u003cp\u003e \u003cb\u003eFunctional diversity of prokaryote microbiota in\u003c/b\u003e \u003cb\u003eLepus\u003c/b\u003e \u003cb\u003espp.\u003c/b\u003e\u003c/p\u003e \u003cp\u003eThe Wilcoxon Rank-Sum on predicted metaCyc pathway abundances suggested there were significant differences in Shannon (functional) diversity between faecal samples of sympatric \u003cem\u003eL. timidus\u003c/em\u003e and \u003cem\u003eL. europaeus\u003c/em\u003e, with \u003cem\u003eL. timidus\u003c/em\u003e showing higher Shannon estimates than \u003cem\u003eL. europaeus\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ea). However, no intraspecific difference in functional diversity estimates was found between elevations (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eb). The clustering of faecal prokaryote communities using NMDS with beta diversity estimates (i.e. Bray-Curtis dissimilarity) on predicted metaCyc pathway abundances highlighted a clear separation between the two \u003cem\u003eLepus\u003c/em\u003e spp. species (PERMANOVA: Bray-Curtis R\u003csup\u003e2\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.498, p-value\u0026thinsp;=\u0026thinsp;0.001; ANOSIM: Bray-Curtis R\u0026thinsp;=\u0026thinsp;0.327, p-value\u0026thinsp;=\u0026thinsp;0.001; Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ec). The clustering of \u003cem\u003eLepus\u003c/em\u003e spp. samples based on pathway abundances (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ed) identified the species as the main clustering variable. Consistently, differential abundance testing comparing the predicted metaCyc pathway abundances across \u003cem\u003eLepus\u003c/em\u003e spp. species and elevations identified 75 pathways, all showing significant differences between \u003cem\u003eL. europaeus\u003c/em\u003e at 2000 m a.s.l. and \u003cem\u003eL. timidus\u003c/em\u003e at 2500 m a.s.l. (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ed, Tab. S7). Of note, 61 out of these 75 pathways displayed significant differences in abundance between \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e at 2000 m a.s.l. as well. Additionally, the contrast between \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e at 2500 m a.s.l. identified 24 metaCyc pathways with significant differences between the two species, the majority of which (70.8%) were also included in the contrast between \u003cem\u003eL. europaeus\u003c/em\u003e at 2000 m a.s.l. and \u003cem\u003eL. timidus\u003c/em\u003e at 2500 m a.s.l. Differential pathways were involved in different classes of biological processes (e.g. superclasses), with Cofactor, Carrier, and Vitamin Biosynthesis (18 pathways, 24% of differential pathways), Fatty Acid and Lipid Biosynthesis (7 pathways, 9% of differential pathways), Carbohydrate Degradation (7 pathways, 9% of differential pathways), Amino Acid Degradation (5 pathways, 7% of differential pathways), and Nucleoside and Nucleotide Degradation (5 pathways, 7% of differential pathways) being the most represented (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ed, Tab. S7). Of note, while most of these superclasses displayed a relatively balanced proportion of pathway abundances between the two species, all differential pathways within the superclasses Fatty Acid and Lipid Biosynthesis and Amino Acid Degradation were significantly enriched in \u003cem\u003eL. timidus\u003c/em\u003e compared to \u003cem\u003eL. europaeus\u003c/em\u003e.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e"},{"header":"Discussion","content":"\u003cp\u003eGut microbiota diversity and composition are known to significantly impact animal health and survival [\u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e56\u003c/span\u003e, \u003cspan citationid=\"CR57\" class=\"CitationRef\"\u003e57\u003c/span\u003e]; thus, investigating gut microbial communities in species potentially threatened by rapid biotic and abiotic environmental changes may enhance our understanding of their conservation status and extinction risk. In this study, we characterized the gut microbiota of the mountain (\u003cem\u003eL. timidus\u003c/em\u003e) and European brown (\u003cem\u003eL. europaeus\u003c/em\u003e) hares using field-collected faecal samples to investigate the consequences of climate change-driven sympatry on their gut prokaryote and fungal communities. Our data confirms the effectiveness of non-invasive samples for the monitoring of gut microbial diversity and composition in elusive species such as \u003cem\u003eLepus\u003c/em\u003e spp. Unexpectedly, despite overlapping in the upper 1000 m of their vertical distribution, our results demonstrated that \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e have distinct faecal microbiota and mycobiota profiles. For the first time, we also identified a significant association between elevation and fungal diversity in both species.\u003c/p\u003e \u003cp\u003e \u003cb\u003eThe composition of\u003c/b\u003e \u003cb\u003eL. timidus\u003c/b\u003e \u003cb\u003eand\u003c/b\u003e \u003cb\u003eL. europaeus\u003c/b\u003e \u003cb\u003efaecal microbiotas are distinct\u003c/b\u003e\u003c/p\u003e \u003cp\u003eThe prokaryote faecal microbiota of the two \u003cem\u003eLepus\u003c/em\u003e spp. were characterized by comparable alpha diversity estimates (Fig. S3), but marked differences in community structure (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e), with up to 27% of observed variation associated with host species (Table\u0026nbsp;\u003cspan refid=\"Tab1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). Accordingly, prokaryote communities displayed distinct taxonomic profiles (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e, \u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e), with only a limited number of shared abundant ASVs (i.e. ~6% ASVs accounting for ~\u0026thinsp;21% sequence reads; Fig. \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e). \u003cem\u003eL. europaeus\u003c/em\u003e faecal samples were primarily dominated by the phylum Firmicutes, and rich in Bacteroidota and Spirochaetota. Instead, the most abundant phylum detected in \u003cem\u003eL. timidus\u003c/em\u003e was Proteobacteria, followed by Bacteroidota and Acidobacteriota. A dominance of Firmicutes and Bacteroidota has already been reported for \u003cem\u003eL. europaeus\u003c/em\u003e faecal samples [\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e] as well as \u003cem\u003eL. granatensis\u003c/em\u003e [\u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e58\u003c/span\u003e] and \u003cem\u003eL. americanus\u003c/em\u003e [\u003cspan citationid=\"CR59\" class=\"CitationRef\"\u003e59\u003c/span\u003e]. Instead, a high relative abundance of Proteobacteria, as seen in \u003cem\u003eL. timidus\u003c/em\u003e, was previously reported for \u003cem\u003eL. sinensis\u003c/em\u003e [\u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e60\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eThe clear differentiation observed between the microbial communities of the two species, particularly at higher elevations where both species were present, suggests that any contaminants from environmental sources may have had only marginal effects on the diversity or composition of faecal bacterial communities. While Firmicutes, Bacteroidota and Proteobacteria are both very commonly found in the mammalian gut, the observed differences in prokaryote composition may reflect co-evolution of host species and associated microbiota after speciation and/or adaptation to different ecological niches [\u003cspan citationid=\"CR61\" class=\"CitationRef\"\u003e61\u003c/span\u003e]. For example, \u003cem\u003eL. europaeus\u003c/em\u003e is primarily a steppe species (with the earliest fossil evidence in the Alpine region dating to approximately 8,000 years ago [\u003cspan citationid=\"CR62\" class=\"CitationRef\"\u003e62\u003c/span\u003e]), while \u003cem\u003eL. timidus\u003c/em\u003e is a boreal-adapted species [\u003cspan citationid=\"CR63\" class=\"CitationRef\"\u003e63\u003c/span\u003e]. While there are no other studies addressing microbial diversity in sympatric Leporidae, the greater impact of phylogeny compared to other environmental factors, including diet, has been already documented in other mammals [\u003cspan citationid=\"CR64\" class=\"CitationRef\"\u003e64\u003c/span\u003e, \u003cspan citationid=\"CR65\" class=\"CitationRef\"\u003e65\u003c/span\u003e, \u003cspan citationid=\"CR66\" class=\"CitationRef\"\u003e66\u003c/span\u003e]. Therefore, despite their coexistence at higher elevations in our Alpine study area, in keeping with their evolutionary history, these species probably exploit distinct resources, with \u003cem\u003eL. europaeus\u003c/em\u003e grazing on the high fat parts of weeds/grasses and various crop types richer in fats and proteins, as noted by Reichlin et al. [\u003cspan citationid=\"CR67\" class=\"CitationRef\"\u003e67\u003c/span\u003e] and Schai-Braun et al. [\u003cspan citationid=\"CR68\" class=\"CitationRef\"\u003e68\u003c/span\u003e]), while \u003cem\u003eL. timidus\u003c/em\u003e may exploit significantly more ligneous plants and Ericaceae [\u003cspan citationid=\"CR69\" class=\"CitationRef\"\u003e69\u003c/span\u003e]. Future studies combining longitudinal vegetation surveys, diet metabarcoding, and metataxonomic investigations are needed to clarify the links existing between seasonal vegetation variation, diet and gut microbial diversity in the two species.\u003c/p\u003e \u003cp\u003eAlthough the effects of seasonal dietary variation over the course of the year on gut microbiota are still largely unknown for these species, the impact of dietary niche on the gut microbiota has been well-documented for many mammal species including humans [\u003cspan citationid=\"CR70\" class=\"CitationRef\"\u003e70\u003c/span\u003e], laboratory mice [\u003cspan citationid=\"CR71\" class=\"CitationRef\"\u003e71\u003c/span\u003e] and wild mammals (see review in: Ley et al. [\u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e58\u003c/span\u003e]; Alessandri et al. [\u003cspan citationid=\"CR72\" class=\"CitationRef\"\u003e72\u003c/span\u003e]). Additionally, the contrasting taxonomic composition of \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e gut microbiota implies different strategies for energy resorption adopted by the two microbial communities [\u003cspan citationid=\"CR64\" class=\"CitationRef\"\u003e64\u003c/span\u003e, \u003cspan citationid=\"CR73\" class=\"CitationRef\"\u003e73\u003c/span\u003e, \u003cspan citationid=\"CR74\" class=\"CitationRef\"\u003e74\u003c/span\u003e]. In fact, our PICRUSt2 analysis supports the hypothesis that the two host-associated microbiotas have alternative biosynthetic and metabolic potential, as indicated by the differences in the number and relative abundance of predicted functions (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e). Although a high abundance of Ericaceae has been reported in the \u003cem\u003eL. timidus\u003c/em\u003e diet [\u003cspan citationid=\"CR69\" class=\"CitationRef\"\u003e69\u003c/span\u003e], we found no differences in metabolic pathways associated with the degradation of phenolic compounds which are abundant in these plants. Interestingly, all differential pathways involved in fatty acid / lipid biosynthesis and amino acid degradation were significantly enriched in \u003cem\u003eL. timidus\u003c/em\u003e, which may again reflect differences in diet between the two species [\u003cspan citationid=\"CR68\" class=\"CitationRef\"\u003e68\u003c/span\u003e, \u003cspan citationid=\"CR69\" class=\"CitationRef\"\u003e69\u003c/span\u003e]. We also speculate that since host-mediated biosynthesis of fatty acids and lipids support metabolic processes that generate heat, thereby aiding in thermoregulation and overall survival in colder climates [\u003cspan citationid=\"CR75\" class=\"CitationRef\"\u003e75\u003c/span\u003e], microbial mediated biosynthesis of these molecules may help \u003cem\u003eL. timidus\u003c/em\u003e to cope with the cold temperatures typically found in its natural home range. Overall, these results suggest that the composition of the gut microbiota may be associated with the adaptation of \u003cem\u003eL. timidus\u003c/em\u003e to its boreal habitat. However, again, additional studies on the actual diet of the individuals sampled in this study are needed to confirm this hypothesis.\u003c/p\u003e \u003cp\u003eWhile the prevalence and relative abundance of Proteobacteria in \u003cem\u003eL. timidus\u003c/em\u003e were consistent across faecal pellets, most \u003cem\u003eL. europaeus\u003c/em\u003e samples were characterized by relatively low abundance of this phylum. However, 15% of investigated faecal samples displayed unexpectedly high abundances (e.g. \u0026gt; 80% of sequence reads) of Proteobacteria, mostly belonging to Enterobacteraceae, Erwiniaceae and Oxalobacteraceae (Fig. \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e). Enterobacteriaceae is the most studied family of these three and includes both commensal bacteria contributing to the maintenance of the gut anaerobic environment, the production of secondary metabolites (e.g. vitamins) and the protection against gut pathogens, as well as opportunistic disease-causing pathogens [\u003cspan citationid=\"CR76\" class=\"CitationRef\"\u003e76\u003c/span\u003e, \u003cspan citationid=\"CR77\" class=\"CitationRef\"\u003e77\u003c/span\u003e]. Although Enterobacteriaceae are commonly found in low abundance in the gut microbiota of healthy mammals, their proliferation in humans is also considered to be a biomarker of gut dysbiosis, which has been associated with several inflammatory bowel diseases [\u003cspan citationid=\"CR78\" class=\"CitationRef\"\u003e78\u003c/span\u003e]. Because only non-invasive faecal samples were collected, we do not know if the observed abundance of Proteobacteria in a limited, but not negligible, number of \u003cem\u003eL. europaeus\u003c/em\u003e samples could be attributable to a state of inflammation of the host. Changes in microbiota composition or diet may also lead to favourable environmental conditions for Proteobacteria (e.g. variation in carbon sources or decreased hypoxia, among others [\u003cspan citationid=\"CR78\" class=\"CitationRef\"\u003e78\u003c/span\u003e]). This unusually high abundance of Proteobacteria should be further investigated, since a disrupted gut microbiota, or microbiota consisting of a number of potentially pathogenic taxa could impact the health of \u003cem\u003eLepus\u003c/em\u003e individuals. Furthermore, future studies should also plan detailed investigations targeting specific Enterobacteriaceae taxa as biomarkers, which could be useful for health monitoring of these species.\u003c/p\u003e \u003cp\u003eSimilarly, both diversity and taxonomic composition of fungal taxa between \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e diverged significantly, (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e), with up to 24% of observed variation, associated with host species (Table\u0026nbsp;\u003cspan refid=\"Tab1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). The fungal classes Tremellomycetes and Cystobasidiomycetes were significantly enriched in mountain hare gut microbiota, while higher abundances of Ustilaginomycetes and Pezizomycetes were detected in the European brown hare. The class Tremellomycetes is a nutritionally heterogeneous group comprising saprotrophs, animal parasites, and fungicolous species previously isolated from various habitats including alpine soil [\u003cspan citationid=\"CR79\" class=\"CitationRef\"\u003e79\u003c/span\u003e], subalpine grasses [\u003cspan citationid=\"CR80\" class=\"CitationRef\"\u003e80\u003c/span\u003e], snow collected at high elevations [\u003cspan citationid=\"CR81\" class=\"CitationRef\"\u003e81\u003c/span\u003e] and the Antarctic Polar Plateau [\u003cspan citationid=\"CR82\" class=\"CitationRef\"\u003e82\u003c/span\u003e]. Due to the ability of some Tremellomycetes to survive in mountainous or cold habitats, we hypothesize that the high abundance of Tremellomycetes found in \u003cem\u003eL. timidus\u003c/em\u003e resulted from these animals having more opportunities to feed on these fungi or have contact with them at higher elevations. However, this fungal class has also been linked to high lipid production [\u003cspan citationid=\"CR83\" class=\"CitationRef\"\u003e83\u003c/span\u003e], which may contribute to energy homeostasis and thermoregulation [\u003cspan citationid=\"CR84\" class=\"CitationRef\"\u003e84\u003c/span\u003e] in \u003cem\u003eL. timidus\u003c/em\u003e. On the other hand, Ustilaginomycetes, found here to be more abundant in \u003cem\u003eL. europaeus\u003c/em\u003e, are well-known to infect vascular plants, especially grass families like Poaceae [\u003cspan citationid=\"CR85\" class=\"CitationRef\"\u003e85\u003c/span\u003e]. Indeed, \u003cem\u003eL. europaeus\u003c/em\u003e are primarily grazers and have a diet that heavily relies on grasses [\u003cspan citationid=\"CR69\" class=\"CitationRef\"\u003e69\u003c/span\u003e], possibly explaining the higher abundance of this fungal class, as well as Pezizomycetes (phylum Ascomycota) in their gut compared to \u003cem\u003eL. timidus\u003c/em\u003e. Again, molecular diet analysis of the faecal pellets would allow us to confirm that the patterns in fungal diversity and composition are mainly due to differences in feeding preferences.\u003c/p\u003e \u003cp\u003eIn this study, host species identification was conducted using mitochondrial DNA (mtDNA) barcoding, which is currently regarded as the most reliable approach for species assignment from faecal samples, owing to the high copy number of mitochondrial genomes per cell and their relatively high mutation rate which allow for good taxonomic discrimination, even between closely related species [\u003cspan citationid=\"CR86\" class=\"CitationRef\"\u003e86\u003c/span\u003e, \u003cspan citationid=\"CR87\" class=\"CitationRef\"\u003e87\u003c/span\u003e]. However, because mtDNA barcoding does not allow discrimination between hybrids and their maternal parental species, we cannot exclude the possibility that a small number of the analysed samples may correspond to hybrid individuals. However, the occurrence in northern European and Alpine regions has been estimated to occur at a frequency of only 2\u0026ndash;5% [\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e, \u003cspan citationid=\"CR88\" class=\"CitationRef\"\u003e88\u003c/span\u003e, \u003cspan citationid=\"CR89\" class=\"CitationRef\"\u003e89\u003c/span\u003e]. Therefore, while hybridization is frequently associated with substantial shifts in the composition and diversity of gut microbial communities [\u003cspan citationid=\"CR90\" class=\"CitationRef\"\u003e90\u003c/span\u003e, \u003cspan citationid=\"CR91\" class=\"CitationRef\"\u003e91\u003c/span\u003e, \u003cspan citationid=\"CR92\" class=\"CitationRef\"\u003e92\u003c/span\u003e], we do not know yet how the observed differentiation in composition and predicted functional potential between \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e could be translated into the gut microbiota of their hybrids. \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e hybrids might exhibit transgressive microbiome configurations, not directly resembling the composition and the predicted functional potential between parental species as seen for hybrids between subspecies \u003cem\u003eMus musculus musculus\u003c/em\u003e \u0026times; \u003cem\u003eM. m. domesticus\u003c/em\u003e [\u003cspan citationid=\"CR90\" class=\"CitationRef\"\u003e90\u003c/span\u003e], intermediate microbiome configurations (e.g. \u003cem\u003eCentropyge\u003c/em\u003e spp.; [\u003cspan citationid=\"CR91\" class=\"CitationRef\"\u003e91\u003c/span\u003e]), or asymmetric, parent-biased microbiomes as seen for \u003cem\u003eCervus elaphus\u003c/em\u003e [\u003cspan citationid=\"CR92\" class=\"CitationRef\"\u003e92\u003c/span\u003e]. Therefore, at present we are unable to predict the most likely microbiota configuration of \u003cem\u003eLepus\u003c/em\u003e hybrids, or its influence on host nutritional behavior, disease susceptibility, or adaptive potential to changing environments. Genotyping of faecal samples could resolve the issue by identifying hybrid individuals; however, many more samples would be needed since the success rate of genotyping from hare faecal samples is only about 50%, and the occurrence of hybrids is so low and unpredictable.\u003c/p\u003e \u003cdiv id=\"Sec15\" class=\"Section2\"\u003e \u003ch2\u003eWithin-species variation in microbial diversity across an elevational gradient\u003c/h2\u003e \u003cp\u003eElevation did not appear to affect the alpha diversity of prokaryote and fungal communities hosted by the two \u003cem\u003eLepus\u003c/em\u003e species, and only marginally impacted (Tab. S3) beta diversity estimates of \u003cem\u003eL. europaeus\u003c/em\u003e prokaryote gut communities. The resilience of the hare gut prokaryote microbiota to variation in elevation and associated environmental variables was unexpected based on studies performed on other mammal species like wild house mice [\u003cspan citationid=\"CR93\" class=\"CitationRef\"\u003e93\u003c/span\u003e], pikas [\u003cspan citationid=\"CR94\" class=\"CitationRef\"\u003e94\u003c/span\u003e], ungulates [\u003cspan citationid=\"CR95\" class=\"CitationRef\"\u003e95\u003c/span\u003e]; macaques, humans and domestic dogs [\u003cspan citationid=\"CR96\" class=\"CitationRef\"\u003e96\u003c/span\u003e]. However, the beta diversity of the prokaryote community structure was in striking contrast to that observed in the gut fungal communities of both species; in fact, marked differences in fungal microbiota composition were detected using both Bray-Curtis (e.g. accounting for ASV abundances) and Jaccard (e.g. based on ASV presence/absence) indices, indicating that variation in fungal composition across elevations reflected changes in both the shared abundant and rare components of the environmental fungal community. Furthermore, we found elevation to be a significant driver of fungal diversity even in pairwise comparisons between neighbouring elevations (e.g. 1000m vs 1500m a.s.l. for \u003cem\u003eL. europaeus\u003c/em\u003e and 2000 m vs 2500 m for \u003cem\u003eL. timidus\u003c/em\u003e). Consistently, Mantel tests highlighted low but significant correlations between the composition of overall fungal communities and pre-sampling precipitation, temperature and plant richness. Moreover, several fungal ASVs in \u003cem\u003eL. europaeus\u003c/em\u003e, including seven matching the coprophilous genus \u003cem\u003eSporormiella\u003c/em\u003e (phylum Ascomycota), that uses herbivore dung as a primary substrate [\u003cspan citationid=\"CR97\" class=\"CitationRef\"\u003e97\u003c/span\u003e], was significantly associated with these environmental variables. Interestingly, variation in the relative abundance of \u003cem\u003eSporomiella\u003c/em\u003e in the gut of the lagomorph \u003cem\u003eOchotona curzoniae\u003c/em\u003e (i.e. plateau pikas) captured at different elevations was also reported by Tang et al. [\u003cspan citationid=\"CR98\" class=\"CitationRef\"\u003e98\u003c/span\u003e]. However, only a few fungal taxa are considered gut residents [\u003cspan citationid=\"CR99\" class=\"CitationRef\"\u003e99\u003c/span\u003e, \u003cspan citationid=\"CR100\" class=\"CitationRef\"\u003e100\u003c/span\u003e] and environmental fungi ingested with diet are almost certainly more exposed to abiotic variation. Consistently, differences in soil fungal communities driven by elevation have been reported frequently [\u003cspan citationid=\"CR79\" class=\"CitationRef\"\u003e79\u003c/span\u003e, \u003cspan citationid=\"CR101\" class=\"CitationRef\"\u003e101\u003c/span\u003e, \u003cspan citationid=\"CR102\" class=\"CitationRef\"\u003e102\u003c/span\u003e]. Additionally, the survival of many fungal taxa relies on interactions with plants and soil [\u003cspan citationid=\"CR103\" class=\"CitationRef\"\u003e103\u003c/span\u003e, \u003cspan citationid=\"CR104\" class=\"CitationRef\"\u003e104\u003c/span\u003e, \u003cspan citationid=\"CR105\" class=\"CitationRef\"\u003e105\u003c/span\u003e]. Therefore, the high number of taxa common to both species (80% of generated sequence reads, Fig. \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e), together with the marked association with elevation and related environmental variables (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e, Tab. S5), indicate a strong contribution of diet (and possibly geophagy) to the gut fungal assembly in these two \u003cem\u003eLepus\u003c/em\u003e species, as suggested by Grieneisen et al. [\u003cspan citationid=\"CR106\" class=\"CitationRef\"\u003e106\u003c/span\u003e] for a primate hybrid zone (see also [\u003cspan citationid=\"CR107\" class=\"CitationRef\"\u003e107\u003c/span\u003e]).\u003c/p\u003e \u003cp\u003eThe clear differentiation between prokaryote and fungal communities in their resilience to environmental parameters could have interesting implications in relation to the observed elevation shifts in the distribution of \u003cem\u003eL. europaeus\u003c/em\u003e and the associated potential reduction in habitat availability for \u003cem\u003eL. timidus\u003c/em\u003e. The lack of differentiation in prokaryote communities across elevations suggests that gut functionality was being maintained along the entire gradient. Our data also suggested that \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e had different feeding strategies (e.g. selecting different plant species) maintained along the same gradient, allowing them to co-exist in the same areas. However, we expect a selective feeding strategy to represent a potential vulnerability for \u003cem\u003eL. timidus\u003c/em\u003e, particularly if climate change continues to alter the distribution of plant taxa in the \u003cem\u003eL. timidus\u003c/em\u003e diet. Furthermore, the unique components of the \u003cem\u003eL. timidus\u003c/em\u003e gut microbiota biodiversity also risk disappearing, with unknown impacts on the Alpine ecosystem. Future studies addressing the health impacts of diet, geophagy and other environmental factors across elevations are needed to disentangle the associations of these factors with the taxonomic and functional diversity of \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e in the Alps, and consequently, of the implications for their survival during climate warming.\u003c/p\u003e \u003c/div\u003e"},{"header":"Conclusions","content":"\u003cp\u003eOur findings demonstrate that sympatric populations of \u003cem\u003eL. europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e maintained distinct, species-specific prokaryote and fungal gut microbiota in terms of taxonomy and structure, but not diversity, at least within our study area. Our data act as a baseline for future monitoring of the process of adaptation of these species to climate warming and associated environmental changes. Future studies should prioritize the characterization of gut prokaryote and fungal composition in hybrids of the two lagomorph species across multiple populations across the Alps. Moreover, future research should explore how climate-driven changes in soil microbiota and food plant composition may shape the future diversity and structure of host-associated gut microbiota in these Lagomorph species.\u003c/p\u003e "},{"header":"Declarations","content":"\u003ch2\u003eEthics approval and consent to participate\u003c/h2\u003e\n\u003cp\u003eNot applicable.\u003c/p\u003e\n\u003ch2\u003eCompeting interests\u003c/h2\u003e\n\u003cp\u003eThe authors declare no competing interests.\u003c/p\u003e\n\u003ch2\u003eConsent for publication\u003c/h2\u003e\n\u003cp\u003eNot applicable.\u003c/p\u003e\n\u003ch2\u003eAdditional information\u003c/h2\u003e\n\u003cp\u003eThe authors declare no competing interests\u003c/p\u003e\n\u003ch2\u003eFunding\u003c/h2\u003e\n\u003cp\u003eThe EUREGIO project: MICROVALU - Evaluating microbiodiversity in alpine pastures (Project ID: IPN94), awarded to PI, JS and HCH, was funded by the \u0026quot;Euregio Tirolo-Alto Adige-Trentino\u0026quot; Interregional Project Network. This study was also partially carried out with funding to HCH at the Fondazione E. Mach (Project BIOALPEC) under the National Biodiversity Future Centre (NBFC) Project (code CN_00000033, Concession Decree No. 1034 of 17 June 2022 adopted by the Italian Ministry of University and Research, CUPD43C22001280006), funded under the National Recovery and Resilience Plan (NRRP), Mission 4 Component 2 Investment 1.4 - Call for tender No. 3138 of 16 December 2021, rectified by Decree n.3175 of 18 December 2021 of Italian Ministry of University and Research funded by the European Union \u0026ndash; NextGenerationEU. This manuscript reflects only the authors\u0026rsquo; views and opinions, neither the European Union nor the European Commission can be considered responsible for them.\u003c/p\u003e\n\u003ch2\u003eAuthor Contribution\u003c/h2\u003e\n\u003cp\u003eB.C., H.C.H., and G.G. conducted the sampling. L.M., B.C. and G.G. completed the laboratory analyses. Data analyses were provided by L.M. and G.G., with support from N.P., T.R., J.S., P.I. and F.N.M.; L.M., H.C.H. and G.G. drafted the manuscript. All authors edited and approved the final manuscript.\u003c/p\u003e\n\u003ch2\u003eAcknowledgement\u003c/h2\u003e\n\u003cp\u003eThe authors wish to thank the Fondazione E. Mach for access to facilities, and the staff of the Sequencing and Genotyping Platform for their outstanding support.\u003c/p\u003e\n\u003ch2\u003eData Availability\u003c/h2\u003e\n\u003cp\u003eSanger sequences have been deposited at NCBI GenBank with accession numbers PX122685 - PX122779. Sanger sequences of the mitochondrial D-loop generated in our study were made available to the editor and reviewers with the uploaded file: Submission2992130.txt.gz. The raw amplicon-sequencing data has been deposited at NCBI Sequence Read Archive (SRA) under the BioProject ID PRJNA1304890. Reviewers can access BioProject and associated SRA metadata at [https://dataview.ncbi.nlm.nih.gov/object/PRJNA1304890?reviewer=rsuoanlfj6uo274kij1dik96rn] .\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eLenoir, J., G\u0026eacute;gout, J. C., Marquet, P. A., de Ruffray, P. \u0026amp; Brisse, H. A. significant upward shift in plant species optimum elevation during the 20th century. \u003cem\u003eScience\u003c/em\u003e \u003cstrong\u003e320\u003c/strong\u003e, 1768\u0026ndash;1771 (2008).\u003c/li\u003e\n\u003cli\u003eHarsch, M. A., Hulme, P. E., McGlone, M. S. \u0026amp; Duncan, R. P. Are treelines advancing? A global meta-analysis of treeline response to climate warming. \u003cem\u003eEcol. Lett.\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, 1040\u0026ndash;1049 (2009).\u003c/li\u003e\n\u003cli\u003eNoce, S., Cipriano, C. \u0026amp; Santini, M. Altitudinal shifting of major forest tree species in Italian mountains under climate change. \u003cem\u003eFront. For. Glob. Change\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, 1250651; https://doi.org/10.3389/ffgc.2023.1250651 (2023).\u003c/li\u003e\n\u003cli\u003eHickling, R., Roy, D. B., Hill, J. K., Fox, R. \u0026amp; Thomas C. D. The distributions of a wide range of taxonomic groups are expanding polewards. \u003cem\u003eGlob. Change Biol.\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, 450-455 (2006).\u003c/li\u003e\n\u003cli\u003ePopy, S., Bordignon L. \u0026amp; Prodon, R. A. weak upward elevational shift in the distributions of breeding birds in the Italian Alps. \u003cem\u003eJ. Biogeogr.\u003c/em\u003e \u003cstrong\u003e37\u003c/strong\u003e, 57-67 (2010).\u003c/li\u003e\n\u003cli\u003eB\u0026uuml;ntgen U. et al. Elevational range shifts in four mountain ungulate species from the Swiss Alps. \u003cem\u003eEcosphere\u003c/em\u003e \u003cstrong\u003e8\u003c/strong\u003e, e01761. https://doi.org/10.1002/ecs2.1761 (2017).\u003c/li\u003e\n\u003cli\u003eSchai-Braun, S. C., Jenny, H., Ruf, T. \u0026amp; Hackl\u0026auml;nder, K. Temperature increase and frost decrease driving upslope elevational range shifts in Alpine grouse and hares. \u003cem\u003eGlob. Change Biol.\u003c/em\u003e \u003cstrong\u003e27\u003c/strong\u003e, 6602\u0026ndash;6614 (2021).\u003c/li\u003e\n\u003cli\u003eRehnus, M. \u003cem\u003eDer Schneehase in den Alpen: Ein \u0026Uuml;berlebensk\u0026uuml;nstler mit ungewisser Zukunft\u003c/em\u003e. Bristol-Stiftung, Z\u0026uuml;rich; Haupt, Bern (2013).\u003c/li\u003e\n\u003cli\u003eStokes, A. W., Hofmeester, T. R., Thorsen, N. H., Odden, J., Linnell, J. D. C. \u0026amp; Pedersen, S. Altitude, latitude and climate zone as determinants of mountain hare (\u003cem\u003eLepus timidus\u003c/em\u003e) coat colour change. \u003cem\u003eEcol. Evol.\u003c/em\u003e \u003cstrong\u003e13\u003c/strong\u003e, e10548. https://doi.org/10.1002/ece3.10548 (2023).\u003c/li\u003e\n\u003cli\u003eSliwinski, K. et al. Habitat requirements of the European brown hare (\u003cem\u003eLepus europaeus\u003c/em\u003e Pallas 1778) in an intensively used agriculture region (Lower Saxony, Germany). \u003cem\u003eBMC Ecology\u003c/em\u003e, 19, 31. https://doi.org/10.1186/s12898-019-0247-7 (2019).\u003c/li\u003e\n\u003cli\u003eLev\u0026auml;nen, R., Thulin, C.-G., Spong, G., \u0026amp; Pohjoism\u0026auml;ki J. L. O. Widespread introgression of \u003cem\u003eL. timidus\u003c/em\u003e genes into Fennoscandian brown hare populations. \u003cem\u003ePLoS ONE\u003c/em\u003e \u003cstrong\u003e13\u003c/strong\u003e, e0191790. https://doi.org/10.1371/journal.pone.0191790 (2018).\u003c/li\u003e\n\u003cli\u003eJansson, G. \u0026amp; Pehrson, \u0026Aring;. The recent expansion of the brown hare (\u003cem\u003eLepus europaeus\u003c/em\u003e) in Sweden with possible implications to the mountain hare (\u003cem\u003eL. timidus\u003c/em\u003e). \u003cem\u003eEur. J. Wildl. Res.\u003c/em\u003e \u003cstrong\u003e53\u003c/strong\u003e, 125\u0026ndash;130 (2007).\u003c/li\u003e\n\u003cli\u003eLa Morgia, V. et al. Global warming is promoting the rapid invasion of the \u003cem\u003eL. timidus\u003c/em\u003e range by the European hare in the Alps. \u003cem\u003eBiodivers. Conserv.\u003c/em\u003e \u003cstrong\u003e32\u003c/strong\u003e, 3875\u0026ndash;3891 (2023).\u003c/li\u003e\n\u003cli\u003eEdwards, J. E. et al. Domesticated equine species and their derived hybrids differ in their faecal microbiota. \u003cem\u003eAnim. Microbiome\u003c/em\u003e \u003cstrong\u003e2\u003c/strong\u003e, Article 27. https://doi.org/10.1186/s42523-020-00027-7 (2020).\u003c/li\u003e\n\u003cli\u003eWei, L. et al. Hybridization altered the gut microbiota of pigs. \u003cem\u003eFront. Microbiol.\u003c/em\u003e \u003cstrong\u003e14\u003c/strong\u003e, 1177947. https://doi.org/10.3389/fmicb.2023.1177947 (2023).\u003c/li\u003e\n\u003cli\u003eRowland, I. et al. Gut microbiota functions: metabolism of nutrients and other food components. \u003cem\u003eEur. J. Nutr.\u003c/em\u003e \u003cstrong\u003e57\u003c/strong\u003e, 1\u0026ndash;24 (2018).\u003c/li\u003e\n\u003cli\u003eSittipo, P. et al. The function of gut microbiota in immune-related neurological disorders: a review. \u003cem\u003eJ. Neuroinflammation\u003c/em\u003e \u003cstrong\u003e19\u003c/strong\u003e, 154. https://doi.org/10.1186/s12974-022-02510-1 (2022).\u003c/li\u003e\n\u003cli\u003eKasarello, K., Cudnoch-Jędrzejewska, A. \u0026amp; Czarzasta, K. Communication of gut microbiota and brain via immune and neuroendocrine signaling. \u003cem\u003eFront. Microbiol.\u003c/em\u003e \u003cstrong\u003e14\u003c/strong\u003e, 1118529. https://doi.org/10.3389/fmicb.2023.1118529 (2023).\u003c/li\u003e\n\u003cli\u003eMoeller, A. H. et al. Experimental evidence for adaptation to species-specific gut microbiota in house mice. \u003cem\u003emSphere\u003c/em\u003e \u003cstrong\u003e4\u003c/strong\u003e, e00387-19. https://doi.org/10.1128/mSphere.00387-19 (2019).\u003c/li\u003e\n\u003cli\u003eHauffe, H. C. \u0026amp; Barelli, C. Conserve the germs: the gut microbiota and adaptive potential. \u003cem\u003eConserv. Genet.\u003c/em\u003e \u003cstrong\u003e20\u003c/strong\u003e, 19\u0026ndash;27 (2019).\u003c/li\u003e\n\u003cli\u003eRojas, C. A., Ram\u0026iacute;rez-Barahona, S., Holekamp, K. E. \u0026amp; Theis, K. R. Host phylogeny and host ecology structure the mammalian gut microbiota at different taxonomic scales. \u003cem\u003eAnim. Microbiome\u003c/em\u003e \u003cstrong\u003e3\u003c/strong\u003e, 33. https://doi.org/10.1186/s42523-021-00094-4 (2021).\u003c/li\u003e\n\u003cli\u003eNeha, S. A., Hanson, J. D., Wilkinson, J. E., Bradley, R. D. \u0026amp; Phillips, C. D. Impacts of host phylogeny, diet, and geography on the gut microbiome of rodents. \u003cem\u003ePLoS ONE\u003c/em\u003e \u003cstrong\u003e20\u003c/strong\u003e, e0316101. https://doi.org/10.1371/journal.pone.0316101 (2025).\u003c/li\u003e\n\u003cli\u003eBarelli, C., Albanese, D., Stumpf, R. M., Asangba, A., Donati, C., Rovero, F. \u0026amp; Hauffe H. C. The gut microbiota communities of wild arboreal and ground-feeding tropical primates are affected differently by habitat disturbance. \u003cem\u003emSystems\u003c/em\u003e \u003cstrong\u003e5\u003c/strong\u003e, e00061-20. https://doi.org/10.1128/mSystems.00061-20 (2020).\u003c/li\u003e\n\u003cli\u003eAmato, K. R. et al. The gut microbiota appears to compensate for seasonal diet variation in the wild black howler monkey (\u003cem\u003eAlouatta pigra\u003c/em\u003e). \u003cem\u003eMicrob. Ecol.\u003c/em\u003e \u003cstrong\u003e69\u003c/strong\u003e, 434\u0026ndash;443 (2015).\u003c/li\u003e\n\u003cli\u003eStalder, G. L. et al. Gut microbiota of the European brown hare (\u003cem\u003eLepus europaeus\u003c/em\u003e). \u003cem\u003eSci. Rep.\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, 2738. https://doi.org/10.1038/s41598-019-39638-9 (2019).\u003c/li\u003e\n\u003cli\u003eZhang, L., Zhan, H., Xu, W., Yan, S. \u0026amp; Ng, S. C. The role of gut mycobiome in health and diseases. \u003cem\u003eTherap. Adv. Gastroenterol.\u003c/em\u003e \u003cstrong\u003e14\u003c/strong\u003e, 17562848211047130. https://doi.org/10.1177/17562848211047130 (2021).\u003c/li\u003e\n\u003cli\u003eZandonai, A. et al. Six years of high-resolution climatic data collected along an elevation gradient in the Italian Alps. \u003cem\u003eSci. Data\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, 751. https://doi.org/10.1038/s41597-024-03580-x (2024).\u003c/li\u003e\n\u003cli\u003eHilpold, A., Seeber, J., Fontana, V., Niedrist, G., Rief, A., Steinwandter, M., Tasser, E. \u0026amp; Tappeiner, U. Decline of rare and specialist species across multiple taxonomic groups after grassland intensification and abandonment\u003cem\u003e. Biodivers. Conserv.\u003c/em\u003e \u003cstrong\u003e27\u003c/strong\u003e, 3729\u0026ndash;3744 (2018).\u003c/li\u003e\n\u003cli\u003ePraeg, N., Pauli, H. \u0026amp; Illmer, P. Microbial diversity in bulk and rhizosphere soil of \u003cem\u003eRanunculus glacialis\u003c/em\u003e along a high-alpine altitudinal gradient. \u003cem\u003eFront. Microbiol.\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, 1429. https://doi.org/10.3389/fmicb.2019.01429 (2019).\u003c/li\u003e\n\u003cli\u003eGalla, G. et al. Mock community as an in situ positive control for amplicon sequencing of microbiotas from the same ecosystem. \u003cem\u003eSci. Rep.\u003c/em\u003e \u003cstrong\u003e13\u003c/strong\u003e, 4056. https://doi.org/10.1038/s41598-023-30916-1 (2023).\u003c/li\u003e\n\u003cli\u003eWilkinson, G. S. \u0026amp; Chapman, A. M. Length and sequence variation in evening bat D-loop mtDNA. \u003cem\u003eGenetics\u003c/em\u003e \u003cstrong\u003e128\u003c/strong\u003e, 607\u0026ndash;617 (1991).\u003c/li\u003e\n\u003cli\u003eMelo-Ferreira, J., Boursot, P., Randi, E., Kryukov, A., Suchentrunk, F., Ferrand, N. \u0026amp; Alves, P. C. The rise and fall of \u003cem\u003eLepus timidus\u003c/em\u003e during Pleistocene glaciations: expansion and retreat with hybridization in the Iberian Peninsula. \u003cem\u003eMol. Ecol.\u003c/em\u003e \u003cstrong\u003e16\u003c/strong\u003e, 605\u0026ndash;618 (2007).\u003c/li\u003e\n\u003cli\u003ePecchioli, E. \u003cem\u003eVariabilit\u0026agrave; e struttura genetica in \u003c/em\u003eLepus timidus\u003cem\u003e: aspetti micro e macrogeografici, filogenetici e adattativi\u003c/em\u003e (unpublished PhD thesis, Universit\u0026agrave; degli Studi di Ferrara) (2010).\u003c/li\u003e\n\u003cli\u003eZachos, F. E., Ben Slimen, H., Hackl\u0026auml;nder, K., Giacometti, M. \u0026amp; Suchentrunk, F. Regional genetic \u003cem\u003ein situ\u003c/em\u003e differentiation despite phylogenetic heterogeneity in Alpine \u003cem\u003eLepus timidus\u003c/em\u003e. \u003cem\u003eJ. Zool.\u003c/em\u003e \u003cstrong\u003e282\u003c/strong\u003e, 47\u0026ndash;53 (2010).\u003c/li\u003e\n\u003cli\u003eKumar S., Stecher G., Li M., Knyaz C, \u0026amp; Tamura K. MEGA X: Molecular Evolutionary Genetics Analysis across computing platforms. \u003cem\u003eMol. Biol. Evol.\u003c/em\u003e \u003cstrong\u003e35, \u003c/strong\u003e1547-1549 (2018).\u003c/li\u003e\n\u003cli\u003eKlindworth, A. et al. Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. \u003cem\u003eNucleic Acids Res.\u003c/em\u003e \u003cstrong\u003e41\u003c/strong\u003e, e1 (2013).\u003c/li\u003e\n\u003cli\u003eWalters, W. et al. Improved prokaryotic 16S rRNA gene (V4 and V4-5) and fungal internal transcribed spacer marker gene primers for microbial community surveys. \u003cem\u003emSystems\u003c/em\u003e \u003cstrong\u003e1\u003c/strong\u003e, e00009-15. https://doi.org/10.1128/mSystems.00009-15 (2015).\u003c/li\u003e\n\u003cli\u003eIhrmark, K. et al. New primers to amplify the fungal ITS2 region \u0026ndash; evaluation by 454-sequencing of artificial and natural communities. \u003cem\u003eFEMS Microbiol. Ecol.\u003c/em\u003e \u003cstrong\u003e82\u003c/strong\u003e, 666\u0026ndash;677 (2012).\u003c/li\u003e\n\u003cli\u003eWhite, T. J., Bruns, T., Lee, S. \u0026amp; Taylor, J. Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics. In: \u003cem\u003ePCR Protocols: A Guide to Methods and Applications\u003c/em\u003e (eds. Innis, M. A., Gelfand, D. H., Sninsky, J. J. \u0026amp; White, T. J.) 315\u0026ndash;322 (Academic Press, San Diego) (1990).\u003c/li\u003e\n\u003cli\u003eCallahan, B. J. et al. DADA2: High-resolution sample inference from Illumina amplicon data. \u003cem\u003eNat. Methods\u003c/em\u003e \u003cstrong\u003e13\u003c/strong\u003e, 581\u0026ndash;583 (2016).\u003c/li\u003e\n\u003cli\u003eQuast, C. et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. \u003cem\u003eNucleic Acids Res.\u003c/em\u003e \u003cstrong\u003e41\u003c/strong\u003e, D590\u0026ndash;D596 (2012).\u003c/li\u003e\n\u003cli\u003eNilsson, R. H. et al. The UNITE database for molecular identification of fungi: handling dark taxa and parallel taxonomic classifications. \u003cem\u003eNucleic Acids Res.\u003c/em\u003e \u003cstrong\u003e47\u003c/strong\u003e, D259\u0026ndash;D264 (2019).\u003c/li\u003e\n\u003cli\u003eBodenhofer, U., Bonatesta, E., Horejs-Kainrath, C. \u0026amp; Hochreiter, S. msa: an R package for multiple sequence alignment. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e31\u003c/strong\u003e, 3997\u0026ndash;3999 (2015).\u003c/li\u003e\n\u003cli\u003eSchliep, K. P. phangorn: phylogenetic analysis in R. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e27\u003c/strong\u003e, 592\u0026ndash;593 (2011).\u003c/li\u003e\n\u003cli\u003eDavis, N. M., Proctor, D. M., Holmes, S. P., Relman, D. A. \u0026amp; Callahan, B. J. Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data. \u003cem\u003eMicrobiome\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, 226. https://doi.org/10.1186/s40168-018-0605-2 (2018).\u003c/li\u003e\n\u003cli\u003eMcMurdie, P. J. \u0026amp; Holmes, S. phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data. \u003cem\u003ePLoS One\u003c/em\u003e \u003cstrong\u003e8\u003c/strong\u003e, e61217. https://doi.org/10.1371/journal.pone.0061217 (2013).\u003c/li\u003e\n\u003cli\u003eLiu, C., Cui, Y., Li, X. \u0026amp; Yao, M. microeco: An R package for data mining in microbial community ecology. \u003cem\u003eFEMS Microbiol. Ecol.\u003c/em\u003e \u003cstrong\u003e97\u003c/strong\u003e, fiaa255 (2021).\u003c/li\u003e\n\u003cli\u003eWickham, H. \u003cem\u003eggplot2: Elegant Graphics for Data Analysis\u003c/em\u003e. Springer-Verlag, New York (2016).\u003c/li\u003e\n\u003cli\u003eNeuwirth, E. \u003cem\u003eRColorBrewer: ColorBrewer Palettes\u003c/em\u003e. R package version 1.1-3, available at https://CRAN.R-project.org/package=RColorBrewer (2022).\u003c/li\u003e\n\u003cli\u003eGarnier, S. et al. viridis(Lite) \u0026ndash; Colorblind-Friendly Color Maps for R. R package version 0.6.5. Available at: https://sjmgarnier.github.io/viridis/ (2024).\u003c/li\u003e\n\u003cli\u003eSegata, N. et al. Metagenomic biomarker discovery and explanation. \u003cem\u003eGenome Biol.\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, R60. https://doi.org/10.1186/gb-2011-12-6-r60 (2011).\u003c/li\u003e\n\u003cli\u003eOksanen, J. et al. vegan: Community Ecology Package. R package version 2.6-4, available at https://CRAN.R-project.org/package=vegan (2022).\u003c/li\u003e\n\u003cli\u003eDouglas, G. M. et al. PICRUSt2 for prediction of metagenome functions. \u003cem\u003eNat. Biotechnol.\u003c/em\u003e \u003cstrong\u003e38\u003c/strong\u003e, 685\u0026ndash;688 (2020).\u003c/li\u003e\n\u003cli\u003eMetsalu, T. \u0026amp; Vilo, J. ClustVis: a web tool for visualizing clustering of multivariate data using Principal Component Analysis and heatmap. \u003cem\u003eNucleic Acids Res.\u003c/em\u003e \u003cstrong\u003e43\u003c/strong\u003e, W566\u0026ndash;W570 (2015).\u003c/li\u003e\n\u003cli\u003eThe GIMP Development Team. GIMP (Version 2.10.18) [Computer software]. https://www.gimp.org (2020).\u003c/li\u003e\n\u003cli\u003eTuddenham, S. \u0026amp; Sears, C. L. The intestinal microbiome and health. \u003cem\u003eCurr. Opin. Infect. Dis.\u003c/em\u003e \u003cstrong\u003e28\u003c/strong\u003e, 464\u0026ndash;470 (2015).\u003c/li\u003e\n\u003cli\u003eZhan, G. et al. Abnormal gut microbiota composition contributes to cognitive dysfunction in SAMP8 mice. \u003cem\u003eAging\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, 1257\u0026ndash;1267. https://doi.org/10.18632/aging.101464 (2018).\u003c/li\u003e\n\u003cli\u003eMart\u0026iacute;nez-Haro, M., Triad\u0026oacute;-Margarit, X., Mateo, R., Vi\u0026ntilde;uela, J. A. \u0026amp; Casamayor, E. O. Taxonomic and functional fingerprints in the gut microbiota of Iberian hare populations (\u003cem\u003eLepus granatensis\u003c/em\u003e) inhabiting areas with organic and conventional farming. \u003cem\u003eEnviron. Pollut.\u003c/em\u003e \u003cstrong\u003e369\u003c/strong\u003e, 125862 (2025).\u003c/li\u003e\n\u003cli\u003eStothart, M. R. et al. Population dynamics and the microbiome in a wild boreal mammal: The snowshoe hare cycle and impacts of diet, season and predation risk. \u003cem\u003eMol. Ecol.\u003c/em\u003e \u003cstrong\u003e34\u003c/strong\u003e, e17629. https://doi.org/10.1111/mec.17629 (2025).\u003c/li\u003e\n\u003cli\u003eCao, Z. et al. Comparison and association of winter diets and gut microbiota using trnL and 16S rRNA gene sequencing for three herbivores in Taohongling, China. \u003cem\u003eGlob. Ecol. Conserv.\u003c/em\u003e \u003cstrong\u003e53\u003c/strong\u003e, e03041. https://doi.org/10.1016/j.gecco.2024.e03041 (2024).\u003c/li\u003e\n\u003cli\u003eKoskella, B. \u0026amp; Bergelson, J. The study of host\u0026ndash;microbiome (co)evolution across levels of selection. \u003cem\u003ePhilos. Trans. R. Soc. B\u003c/em\u003e \u003cstrong\u003e375\u003c/strong\u003e, 20190604 (2020).\u003c/li\u003e\n\u003cli\u003eBock, A. \u003cem\u003eLepus europaeus\u003c/em\u003e (Lagomorpha: Leporidae). \u003cem\u003eMamm. Species\u003c/em\u003e \u003cstrong\u003e52\u003c/strong\u003e, 125\u0026ndash;142. https://doi.org/10.1093/mspecies/seaa010 (2020).\u003c/li\u003e\n\u003cli\u003eAngerbj\u0026ouml;rn, A. \u0026amp; Flux, J. E. C. \u003cem\u003eLepus timidus\u003c/em\u003e. \u003cem\u003eMamm. Species\u003c/em\u003e \u003cstrong\u003e495\u003c/strong\u003e, 1\u0026ndash;11 (1995).\u003c/li\u003e\n\u003cli\u003eLey, R. E. et al. Evolution of mammals and their gut microbes. \u003cem\u003eScience\u003c/em\u003e \u003cstrong\u003e320\u003c/strong\u003e, 1647\u0026ndash;1651 (2008).\u003c/li\u003e\n\u003cli\u003eKnowles, S. C. L., Eccles, R. M. \u0026amp; Baltrūnaitė, L. Species identity dominates over environment in shaping the microbiota of small mammals. \u003cem\u003eEcol. Lett.\u003c/em\u003e \u003cstrong\u003e22\u003c/strong\u003e, 826\u0026ndash;837 (2019).\u003c/li\u003e\n\u003cli\u003eAnders, J. L. et al. Comparing the gut microbiome along the gastrointestinal tract of three sympatric species of wild rodents. \u003cem\u003eSci. Rep.\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, 19929. https://doi.org/10.1038/s41598-021-99379-6 (2021).\u003c/li\u003e\n\u003cli\u003eReichlin, T., Klansek, E. \u0026amp; Hackl\u0026auml;nder, K. Diet selection by hares (\u003cem\u003eLepus europaeus\u003c/em\u003e) in arable land and its implications for habitat management. \u003cem\u003eEur. J. Wildl. Res.\u003c/em\u003e \u003cstrong\u003e52\u003c/strong\u003e, 109\u0026ndash;118 (2006).\u003c/li\u003e\n\u003cli\u003eSchai-Braun, S. C. et al\u003cem\u003e.\u003c/em\u003e The European hare (\u003cem\u003eLepus europaeus\u003c/em\u003e): a picky herbivore searching for plant parts rich in fat. \u003cem\u003ePLoS ONE\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, e0134278. https://doi.org/10.1371/journal.pone.0134278 (2015).\u003c/li\u003e\n\u003cli\u003eSchai-Braun, S. C. et al. Diet niches of the Alpine \u003cem\u003eLepus timidus\u003c/em\u003e (\u003cem\u003eLepus timidus varronis\u003c/em\u003e) and European hare (\u003cem\u003eLepus europaeus\u003c/em\u003e) living in coexistence in the Alps. \u003cem\u003eMamm. Biol.\u003c/em\u003e \u003cstrong\u003e105\u003c/strong\u003e, 139\u0026ndash;145 (2025).\u003c/li\u003e\n\u003cli\u003eDavid, L. A. et al\u003cem\u003e.\u003c/em\u003e Diet rapidly and reproducibly alters the human gut microbiome. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e505\u003c/strong\u003e, 559\u0026ndash;563 (2014).\u003c/li\u003e\n\u003cli\u003eCarmody, R. N. et al. Diet dominates host genotype in shaping the murine gut microbiota. \u003cem\u003eCell Host Microbe\u003c/em\u003e \u003cstrong\u003e17\u003c/strong\u003e, 72\u0026ndash;84 (2015).\u003c/li\u003e\n\u003cli\u003eAlessandri, G., Rizzo, S. M., Ossiprandi, M. C., van Sinderen, D. \u0026amp; Ventura, M. Creating an atlas to visualize the biodiversity of the mammalian gut microbiota. \u003cem\u003eCurr. Opin. Biotechnol.\u003c/em\u003e \u003cstrong\u003e73\u003c/strong\u003e, 28\u0026ndash;33 (2022).\u003c/li\u003e\n\u003cli\u003eZoetendal, E. G. \u0026amp; de Vos, W. M. Effect of diet on the intestinal microbiota and its activity. \u003cem\u003eCurr. Opin. Gastroenterol.\u003c/em\u003e \u003cstrong\u003e30\u003c/strong\u003e, 189\u0026ndash;195 (2014).\u003c/li\u003e\n\u003cli\u003eSommer, F. et al. The gut microbiota modulates energy metabolism in the hibernating brown bear \u003cem\u003eUrsus arctos\u003c/em\u003e. \u003cem\u003eCell Rep.\u003c/em\u003e \u003cstrong\u003e14\u003c/strong\u003e, 1655\u0026ndash;1661 (2016).\u003c/li\u003e\n\u003cli\u003eXu, Z. et al. Cold-induced lipid dynamics and transcriptional programs in white adipose tissue. \u003cem\u003eBMC Biol.\u003c/em\u003e \u003cstrong\u003e17\u003c/strong\u003e, 74. https://doi.org/10.1186/s12915-019-0693-x (2019).\u003c/li\u003e\n\u003cli\u003eV\u0026aacute;zquez-Baeza, Y. et al. Dog and human inflammatory bowel disease rely on overlapping yet distinct dysbiosis networks. \u003cem\u003eNat. Microbiol.\u003c/em\u003e \u003cstrong\u003e1\u003c/strong\u003e, 16177 (2016).\u003c/li\u003e\n\u003cli\u003eLupp, C. et al\u003cem\u003e.\u003c/em\u003e Host-mediated inflammation disrupts the intestinal microbiota and promotes the overgrowth of \u003cem\u003eEnterobacteriaceae\u003c/em\u003e. \u003cem\u003eCell Host Microbe\u003c/em\u003e \u003cstrong\u003e2\u003c/strong\u003e, 119\u0026ndash;129 (2007).\u003c/li\u003e\n\u003cli\u003eShin, N. R., Whon, T. W. \u0026amp; Bae, J.-W. Proteobacteria: microbial signature of dysbiosis in gut microbiota. \u003cem\u003eTrends Biotechnol.\u003c/em\u003e \u003cstrong\u003e33\u003c/strong\u003e, 496\u0026ndash;503 (2015).\u003c/li\u003e\n\u003cli\u003ePraeg, N.\u003cem\u003e \u003c/em\u003eet al. The role of land management and elevation in shaping soil microbial communities: Insights from the Central European Alps. \u003cem\u003eSoil Biol. Biochem.\u003c/em\u003e \u003cstrong\u003e150\u003c/strong\u003e, 107951 (2020).\u003c/li\u003e\n\u003cli\u003eMouhamadou, B. et al. Differences in fungal communities associated to \u003cem\u003eFestuca paniculata\u003c/em\u003e roots in subalpine grasslands. \u003cem\u003eFungal Divers.\u003c/em\u003e \u003cstrong\u003e47\u003c/strong\u003e, 55\u0026ndash;63 (2011).\u003c/li\u003e\n\u003cli\u003eIrwin, N. A., Twynstra, C. S., Mathur, V. \u0026amp; Keeling, P. J. The molecular phylogeny of \u003cem\u003eChionaster nivalis\u003c/em\u003e reveals a novel order of psychrophilic and globally distributed Tremellomycetes (Fungi, Basidiomycota). \u003cem\u003ePLoS ONE\u003c/em\u003e \u003cstrong\u003e16\u003c/strong\u003e, e0247594. https://doi.org/10.1371/journal.pone.0247594 (2021).\u003c/li\u003e\n\u003cli\u003eStoppiello, G. A., Coleine, C., Moeller, R., Ripa, C., Billi, D., \u0026amp; Selbmann, L. Seasonality is the main determinant of microbial diversity associated to snow/ice around Concordia Station on the Antarctic Polar Plateau. \u003cem\u003eBiology\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, 1193. https://doi.org/10.3390/biology12091193 (2023).\u003c/li\u003e\n\u003cli\u003eKurtzman, C. P., Fell, J. W., \u0026amp; Boekhout, T. The Yeasts: A Taxonomic Study. Fifth Edition (2011).\u003c/li\u003e\n\u003cli\u003eOlsen, L., Thum, E., \u0026amp; Rohner, N. Lipid metabolism in adaptation to extreme nutritional challenges. \u003cem\u003eDevelop. Cell.\u003c/em\u003e \u003cstrong\u003e56 \u003c/strong\u003e(2021).\u003c/li\u003e\n\u003cli\u003eBauer, R., Begerow, D., \u0026amp; Oberwinkler, F. Ustilaginomycotina: The true smut fungi. \u003cem\u003eTree of Life Web Project\u003c/em\u003e available at: http://tolweb.org/Ustilaginomycotina/20530/2008.01.23 (2006).\u003c/li\u003e\n\u003cli\u003eMurphy, M. A., Waits, L. E., \u0026amp; Kendall, K. C. Quantitative evaluation of fecal drying methods for brown bear DNA analysis. \u003cem\u003eWildl. Soc. Bull.\u003c/em\u003e \u003cstrong\u003e28\u003c/strong\u003e, 951\u0026ndash;957 (2000).\u003c/li\u003e\n\u003cli\u003eElyasigorji, Z., Izadpanah, M., Hadi, F., \u0026amp; Zare, M. Mitochondrial genes as strong molecular markers for species identification. \u003cem\u003eNucleus\u003c/em\u003e \u003cstrong\u003e66\u003c/strong\u003e, 81\u0026ndash;93. https://doi.org/10.1007/s13237-022-00393-4 (2022).\u003c/li\u003e\n\u003cli\u003eThulin, C.-G., Stone, J., Tegelstr\u0026ouml;m, H., \u0026amp; Walker, C. W. Species assignment and hybrid identification among Scandinavian hares \u003cem\u003eLepus europaeus\u003c/em\u003e and \u003cem\u003eL. timidus\u003c/em\u003e. \u003cem\u003eWildl. Biol., \u003c/em\u003e\u003cstrong\u003e12\u003c/strong\u003e, 29\u0026ndash;38 (2006).\u003c/li\u003e\n\u003cli\u003eSchai-Braun, S. C., Schwienbacher, S., Smith, S., \u0026amp; Hackl\u0026auml;nder, K. Coexistence of European hares and Alpine mountain hares in the Alps: what drives the occurrence and frequency of their hybrids? \u003cem\u003eJ. Zool. \u003c/em\u003e\u003cstrong\u003e320\u003c/strong\u003e(1), 15\u0026ndash;25\u003cem\u003e \u003c/em\u003e(2023).\u003c/li\u003e\n\u003cli\u003eWang, J. et al. Analysis of intestinal microbiota in hybrid house mice reveals evolutionary divergence in a vertebrate hologenome. \u003cem\u003eNat. Commun.\u003c/em\u003e 6:6440. https//:doi:10.1038/ncomms7440 (2015)\u003c/li\u003e\n\u003cli\u003eHuggett, M. J. et al. Gut microbial communities of hybridising pygmy angelfishes reflect species boundaries. \u003cem\u003eCommun. Biol.\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, 542 (2023). \u003c/li\u003e\n\u003cli\u003eWei, L. et al. Hybridization alters red deer gut microbiome and metabolites. \u003cem\u003eFront. Microbiol.\u003c/em\u003e \u003cstrong\u003e15\u003c/strong\u003e, 1387957 (2024).\u003c/li\u003e\n\u003cli\u003eSuzuki, T. A., Martins, F. M., \u0026amp; Nachman, M. W. Altitudinal variation of the gut microbiota in wild house mice. \u003cem\u003eMol. Ecol.\u003c/em\u003e \u003cstrong\u003e28\u003c/strong\u003e (2019).\u003c/li\u003e\n\u003cli\u003eLi, H., Zhou, R., Zhu, J., Huang, X., \u0026amp; Qu, J. Environmental filtering increases with elevation for the assembly of gut microbiota in wild pikas. \u003cem\u003eMicrob. Biotechnol.\u003c/em\u003e \u003cstrong\u003e12 \u003c/strong\u003e(2019).\u003c/li\u003e\n\u003cli\u003eWang, X., et al. High-altitude drives the convergent evolution of alpha diversity and indicator microbiota in the gut microbiomes of ungulates. \u003cem\u003eFront. Microbiol. \u003c/em\u003e\u003cstrong\u003e13\u003c/strong\u003e https://doi.org/10.3389/fmicb.2022.953234 (2022).\u003c/li\u003e\n\u003cli\u003eZhao, J. et al. Diet and high altitude strongly drive convergent adaptation of gut microbiota in wild macaques, humans, and dogs to high altitude environments. \u003cem\u003eFront. Microbiol. \u003c/em\u003e\u003cstrong\u003e14 \u003c/strong\u003ehttps://doi.org/10.3389/fmicb.2023.1067240 (2023).\u003c/li\u003e\n\u003cli\u003ePerrotti, A. G., \u0026amp; van Asperen, E. Dung fungi as a proxy for megaherbivores: opportunities and limitations for archaeological applications. \u003cem\u003eVeg. Hist. Archaeobot.\u003c/em\u003e \u003cstrong\u003e28\u003c/strong\u003e, 93\u0026ndash;104 (2019).\u003c/li\u003e\n\u003cli\u003eTang, X., Zhang, L., Ren, S., Zhao, Y., Liu, K., \u0026amp; Zhang, Y. Stochastic processes derive gut fungi community assembly of plateau pikas (\u003cem\u003eOchotona curzoniae\u003c/em\u003e) along altitudinal gradients across warm and cold seasons. \u003cem\u003eJ. Fungi \u003c/em\u003e\u003cstrong\u003e9\u003c/strong\u003e\u003cstrong\u003e\u003cem\u003e \u003c/em\u003e\u003c/strong\u003ehttps://doi.org/10.3390/jof9101032 (2023).\u003c/li\u003e\n\u003cli\u003eAuchtung, T. A. et al. Investigating colonization of the healthy adult gastrointestinal tract by fungi. \u003cem\u003emSphere \u003c/em\u003e\u003cstrong\u003e3\u003c/strong\u003ehttps://doi.org/10.1128/mSphere.00092-18 (2018).\u003c/li\u003e\n\u003cli\u003eSuhr, M. J., \u0026amp; Hallen-Adams, H. E. The human gut mycobiome: pitfalls and potentials - a mycologist\u0026rsquo;s perspective. \u003cem\u003eMycologia\u003c/em\u003e \u003cstrong\u003e107\u003c/strong\u003e, 1057-73 (2015).\u003c/li\u003e\n\u003cli\u003eLi, J., Chen, L., Wang, H., Ouyang, S., Liu, X. \u0026amp; Lu, J. Pattern and drivers of soil fungal community along elevation gradient in the \u003cem\u003eAbies georgei\u003c/em\u003e forests of Segila mountains, Southeast Tibet. \u003cem\u003eGlob. Ecol. Conserv.\u003c/em\u003e \u003cstrong\u003e39\u003c/strong\u003e, e02446. https://doi.org/10.1016/j.gecco.2022.e02291 (2022).\u003c/li\u003e\n\u003cli\u003eV\u0026eacute;lez-Mart\u0026iacute;nez, G. A. et al. Soil bacteria and fungi communities are shaped by elevation influences in Colombian forest and p\u0026aacute;ramo natural ecosystems. \u003cem\u003eInt. Microbiol.\u003c/em\u003e \u003cstrong\u003e27\u003c/strong\u003e, 377-391 (2024).\u003c/li\u003e\n\u003cli\u003eDean, R. et al. The Top 10 fungal pathogens in molecular plant pathology. \u003cem\u003eMol. \u003c/em\u003e\u003cem\u003ePlant Pathol. \u003c/em\u003e\u003cstrong\u003e13\u003c/strong\u003e, 414\u0026ndash;430 (2012).\u003c/li\u003e\n\u003cli\u003eRodriguez, R. J., White, J. F. Jr., Arnold, A. E., \u0026amp; Redman, R. S. Fungal endophytes: diversity and functional roles. \u003cem\u003eNew Phytol. \u003c/em\u003e\u003cstrong\u003e182\u003c/strong\u003e, 314\u0026ndash;330 (2009).\u003c/li\u003e\n\u003cli\u003eSemchenko, M. et al. Fungal diversity regulates plant-soil feedbacks in temperate grassland. \u003cem\u003eSci. Adv. \u003c/em\u003e\u003cstrong\u003e4\u003c/strong\u003e https://doi.org/10.1126/sciadv.aau4578 (2018).\u003c/li\u003e\n\u003cli\u003eGrieneisen, L. E. et al. Genes, geology and germs: gut microbiota across a primate hybrid zone are explained by site soil properties, not host species. \u003cem\u003eProc. R. Soc. B\u003c/em\u003e\u003cem\u003e \u003c/em\u003e\u003cstrong\u003e286\u003c/strong\u003e, 20190431 (2019).\u003c/li\u003e\n\u003cli\u003eBanerjee, S., \u0026amp; van der Heijden, M. G. A. Soil microbiomes and one health. \u003cem\u003eNat. Rev. Microbiol.\u003c/em\u003e \u003cstrong\u003e21\u003c/strong\u003e, 6\u0026ndash;20 (2023).\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"scientific-reports","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"scirep","sideBox":"Learn more about [Scientific Reports](http://www.nature.com/srep/)","snPcode":"","submissionUrl":"","title":"Scientific Reports","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"stoa","reportingPortfolio":"Scientific Reports","inReviewEnabled":true,"inReviewRevisionsEnabled":true},"keywords":"Lepus europaeus, Lepus timidus, metataxonomy, mtDNA, 16S rRNA gene, ITS2","lastPublishedDoi":"10.21203/rs.3.rs-7375012/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-7375012/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"The mountain hare (Lepus timidus) is an arctic-alpine species with relictual populations in the Italian Alps, typically occurring at elevations above 2000 m a.s.l. This species is threatened by habitat loss and fragmentation, and declining snow cover due to climate warming. Moreover, as treelines shift upward, the European brown hare (L. europaeus) is expanding its distribution into areas previously dominated by the mountain hare, potentially leading to resource competition, and loss of local adaptation through hybridization and inter-specific gene flow. In particular, the consequences of sympatry on diversity and composition of prokaryote and fungal communities of the gut microbiota, which are critical to individual health, are currently unknown.\nHere, we compared the gut microbiota of these two hare species in an area of overlap in the central Alps by analysing fresh faecal pellets collected from Val Mazia/Matschertal, Italy along an elevational gradient (1000 to 2500 m a.s.l.). For the first time, we describe the prokaryote diversity and composition of L. timidus, and the fungal gut communities (mycobiota) of both Lepus species. Species identity was confirmed for 95 samples via mtDNA barcoding, while gut microbiota richness and composition were investigated using amplicon sequencing, targeting the V3-V4 region of the prokaryote 16S rRNA gene and fungal ITS2 regions.\nDistinct prokaryote and fungal communities were observed for each species, even in sympatry, indicating differences in their functional diversity. Interestingly, for both Lepus species, elevation influenced fungal but not prokaryote diversity. Therefore, sympatry appears to have had minimal impact on gut microbiota composition of either species thus far. Given the expected upward range shift of L. europaeus under climate warming and its continued restocking for hunting, our findings provide an important baseline for assessing the health and adaptability of L. timidus as well as the effectiveness of conservation efforts aimed at protecting L. timidus. However, expanding this research to other areas of sympatry will be essential to understand if gut microbial composition is indicative of L. timidus conservation status across its range.","manuscriptTitle":"Sympatric Lepus spp. in the central Italian Alps host significantly different gut microbiotas","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2026-03-02 14:08:33","doi":"10.21203/rs.3.rs-7375012/v1","editorialEvents":[{"type":"communityComments","content":0},{"type":"decision","content":"Accepted","date":"2026-03-12T07:29:11+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2026-03-09T17:29:56+00:00","index":"","fulltext":""},{"type":"editorInvited","content":"","date":"2026-03-04T11:26:57+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2026-02-07T05:21:18+00:00","index":"","fulltext":""},{"type":"submitted","content":"Scientific Reports","date":"2026-02-04T15:40:42+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"scientific-reports","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"scirep","sideBox":"Learn more about [Scientific Reports](http://www.nature.com/srep/)","snPcode":"","submissionUrl":"","title":"Scientific Reports","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"stoa","reportingPortfolio":"Scientific Reports","inReviewEnabled":true,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"58c860e7-1a7d-4646-846d-51328cfe016b","owner":[],"postedDate":"March 2nd, 2026","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"published-in-journal","subjectAreas":[],"tags":[],"updatedAt":"2026-04-07T16:01:12+00:00","versionOfRecord":{"articleIdentity":"rs-7375012","link":"https://doi.org/10.1038/s41598-026-44592-4","journal":{"identity":"scientific-reports","isVorOnly":false,"title":"Scientific Reports"},"publishedOn":"2026-04-04 15:58:16","publishedOnDateReadable":"April 4th, 2026"},"versionCreatedAt":"2026-03-02 14:08:33","video":"","vorDoi":"10.1038/s41598-026-44592-4","vorDoiUrl":"https://doi.org/10.1038/s41598-026-44592-4","workflowStages":[]},"version":"v1","identity":"rs-7375012","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-7375012","identity":"rs-7375012","version":["v1"]},"buildId":"XKTyCvWXoU3ODBz1xrDgd","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}