Ninety (N=90) adult individuals of three species of the genus Onthophagus (Palaeonthophagus) were collected in the Western Italian Alps [Susa valley, Condove (TO), Piedmont (45.136° N, 7.296° E)] in May 2023. The three species, namely O. fracticornis (OFT), O. medius (OMD) and O. verticicornis (OVT), were identified considering their morphological traits (35) and through molecular characterization, using the mitochondrial cytochrome C oxidase subunit I (COI) gene sequence as a DNA barcode ( Supplementary File 1 , Supplementary Table S1 ).

We collected 30 individuals (15 males and 15 females) of each species: 10 from cow dung, 10 from sheep dung and 10 from donkey dung, with each dung type located in a distinct pasture. Two soil and two dung samples were collected from each of the three pastures, where cattle, sheep, and donkeys were being grazed, respectively: two dung pats from each pasture were chosen at random, and the soil samples were obtained from approximately 5 cm below the surface of the ground located next to these dung pats. The dung and soil samples were stored in 1.5 ml Eppendorf tubes. The two samples of the same dung or soil were subsequently homogenized in the laboratory prior to analysis. Absolute ethanol was added to the tubes containing dung samples to prevent dung fermentation. All samples were then stored at -20°C.

Dung beetles were housed in plastic terraria without food for at least 24 hours to make them excrete as much ingested dung as possible and thus clean out their guts. Individuals were euthanized by submersion in absolute ethanol and immediately dissected to extract the entire gut. The dissection tools were sterilized using a 30% sodium hypochlorite solution and then washed in distilled water. Once removed, the gut was preserved in absolute ethanol and stored at 4°C.

We extracted DNA from gut and dung samples using the CTAB modified method described in Natta et al. (19), while DNA from soil samples was extracted using the DNeasy^® PowerSoil^® Pro Kit (QIAGEN). A DNA metabarcoding approach was used to investigate microbiota: for the prokaryotic component, the 16S ribosomal RNA (rRNA) gene was amplified using the primer set 515fB (5′–GTGYCAGCMGCCGCGGTAA–3′) (36) and 806rB (5′–GGACTACNVGGGTWTCTAAT–3′) (37); for the fungal component, the nuclear ribosomal ITS2 region was amplified using the primer pair fITS7 (5′– GAACGCAGCRAAIIGYGA–3′) and ITS4 (5′–TCCTCCGCTTATTGATATGC–3′) (38). The following Illumina overhang adapter sequences were added to the primer pairs: forward overhang: 5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG-[locus specific target primer]; reverse overhang: 5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG-[locus specific target primer].

PCR reactions were run in a final volume of 25 μl using 1 U of XtraTaq Pol White DNA polymerase (GeneSpin Srl, Milano, Italy), 5x XtraWhite Buffer with MgCl₂, 0.2 μM of each dNTP, 0.5 μM of each primer and 20 ng of genomic DNA. For the prokaryotic community, the PCR cycling program consisted of an initial step at 94°C for 3 min, 35 cycles at 94°C for 45 s, 55°C for 60 s and 72°C for 90 s, followed by a final extension step of 72°C for 10 min. For the fungal community, the PCR cycling program consisted of an initial step at 94°C for 5 min, 35 cycles at 94°C for 30 s, 50°C for 30 s and 72°C for 30 s, and a final extension step of 72°C for 7 min.

Extracted DNA was amplified in triplicate and pooled before purification using Wizard^® SV Gel and the PCR Clean-Up System (Promega). PCR purified products were quantified using the Qubit dsDNA BR Assay kit and Qubit Fluorometer 2.0 following the manufacturer’s protocol and sent to IGA technologies (Udine, Italy) for Illumina NovaSeq sequencing (2 × 250 bp).

The following bioinformatics and statistical analyses were conducted for prokaryotes (i.e. archaea and bacteria) and fungi separately. Sequencing adapters and primers were removed, and then the sequences were analyzed using the microbiome bioinformatics platform QIIME2 (Quantitative Insights Into Microbial Ecology 2, v. 2021.2 (39). Denoising and quality control, including chimaera removal, were performed using the DADA2 plugin (40) through the qiime dada2 denoise-paired command, with chimaera detection carried out using the “consensus” method. The taxonomic assignment of the prokaryotic community was achieved using the Silva 138 99% OTUs full-length sequences database (41, 42), whereas for fungi we used the UNITE Community (2019): UNITE QIIME release for fungi v.04.04.2024 (43). Phylogenetic trees were generated by the QIIME2 plugin “qiime phylogeny align-to-tree-mafft-fasttree”. The outputs of the QIIME2 pipeline “taxonomy.qza”, “otu_table.qza” and “rooted-tree.qza”, together with their metadata files, were then imported into Rstudio (44) to create phyloseq objects using the R package qiime2R v.0.99.6 (45). Shared amplicon sequence variants (ASVs) between the beetles and the environmental samples (i.e., dung and soil samples) were then sought by analyzing the results of the flower plots generated from all 96 samples, using the function plot_venn in the R package microeco (46). The shared ASVs were removed from the phyloseq object.

Following the removal of shared ASVs, we used the phyloseq objects for the following diversity analyses. To allow for the comparisons of samples with non-uniform coverage, we normalized the ASV tables using the rarefy_even_depth function of the R package phyloseq v.1.36.0 (47). We used the rarecurve function of the R package vegan v2.6-2 (48) to obtain rarefaction curves of the rarefied ASV table. From the ASV table, we calculated absolute counts and relative abundances (i.e. the ratio between the number of reads belonging to the ASV in a specific sample and the total number of reads in the sample) for each ASV.

We evaluated alpha diversity using “Observed ASVs”, “Shannon” and “Faith’s Phylogenetic Distance (Faith PD)” indices using the estimate_richness function of the R package phyloseq. We then tested the effect of species, dung type, sex and their interaction on alpha diversity by analysis of variance (ANOVA). We checked the normal distribution of the model dataset by performing the Shapiro-Wilk test on the residuals. This analysis was not performed on the Observed ASVs because of the highly positive correlation with Shannon and Faith PD.

For both prokaryotes and fungi, we visualized dissimilarity between individuals by means of non-metric multidimensional scaling (nMDS) based on a Bray-Curtis distance matrix of ASV composition using R package vegan. Stress was used as a measure of goodness of fit. We plotted the results of the nMDS using the tidyverse collection of R packages (49), first in a general plot for all 90 individuals and then by dividing the plot into three distinct sub-plots to reveal any differences better, considering the 30 individuals belonging to each of the three species or collected in the three dung types. We performed an envfit analysis on the nMDS ordination to determine which of the factors had the greatest influence in differentiating individuals based on the R² values and p-values.Differentiating microbial variability due to dung preference from that due to species is very difficult because these two factors were inseparably associated (we collected the individuals of a species in a certain type of dung). Therefore, the best way to study them was to investigate the interindividual variability in all possible associations between each species and each type of dung. To do so, we considered the Cartesian product of two sets, A and B, which consists of all the ordered pairs that can be constructed with the first element coming from the first set, A, and the second element coming from the second set, B. In our case, one set comprised the three Onthophagus species, and the other set comprised the three dung types. The 9 pairs (3 species x 3 dung types) were therefore: OFT-cow, OFT-sheep, OFT-donkey, OMD-cow, OMD-sheep, OMD-donkey, OVT-cow, OVT-sheep and OVT-donkey. We considered these 9 pairs in the analyses of shared, exclusive and indicator ASVs, and in taxonomic and functional analyses of the gut microbiota of dung beetles.

To visualize the shared (i.e. core) and exclusive ASVs of the nine pairs, we used flower plots generated using the R package microeco (46). In addition, we investigated the indicator ASVs of each of the nine pairs using the multipatt function in the R package indicspecies (50). To visualize the taxonomic diversity of the different samples, we used percent stacked bar charts generated using the tidyverse collection of R packages (49).

The trophic behaviors of the identified communities of prokaryotes and fungi were assessed using FAPROTAX (Functional Annotation of Prokaryotic Taxa) v. 1.2.6 (51) and FUNGuild (52), respectively, both of which were implemented in the R package microeco (46). As done for the taxonomic diversity, we used percent stacked bar charts to visualize the differences in trophic behavior between the samples. The effect of species and dung type on taxonomic and functional composition was tested using ANOVA and the Tukey honestly significant difference (HSD) test.

COI sequences confirmed that the three species (OFT, OMD and OVT) were correctly identified using external morphological traits. All three species considered in this work fell into the same cluster, namely the subgenus Palaeonthophagus ( Supplementary File 1 , Supplementary Figures S1 , S2 ).

The bioinformatics analysis gave rise to 4659 prokaryotic ASVs and 4549 fungal ASVs. The “Observed ASVs” rarefaction curves showed that all 90 beetle samples reached the asymptote for both prokaryotes and fungi. Thus, the depth of sequencing was adequate to represent the diversity of the gut microbiota ( Supplementary Figures S3A, B ). The rarefied number of reads per sample was 1193 for prokaryotes and 1901 for fungi.

Regarding the prokaryotes, ASVs ranged from a minimum of 24 in a male OMD individual collected from donkey dung to a maximum of 395 in a male OMD collected from cow dung ( Supplementary Table S2 ). We found significant differences in the Shannon index based on sex (F = 4.37, p < 0.05) and the interaction term sex-dung type (F = 4.12, p < 0.05). Significant differences were also highlighted in the Faith PD index based on the interaction term species-sex (F = 5.02, p < 0.01) and dung type (F = 3.40, p < 0.05).

Regarding the fungi, the number of ASVs detected ranged from 20 in a male OFT individual collected from donkey dung to a maximum of 399 in a male OMD, once again collected from donkey dung ( Supplementary Table S3 ). The Shannon index varied significantly according to dung type (F = 11.75, p < 0.001), the interaction term species-sex-dung type (F = 5.12, p < 0.01), and on the interaction term species-sex (F = 4.26, p < 0.05). We found significant differences in the Faith PD index based on dung type (F = 8.64, p < 0.001) and the interaction term species-sex-dung type (F = 4.36, p < 0.01).

The ordination plot of the 90 individuals according to the prokaryotic component of their microbiota showed a clearer separation of points in the nMDS space when we grouped individuals according to species ( Figure 1A ) instead of dung type ( Supplementary Figure S4A ). The envfit analysis confirmed this result as species showed the highest R² value (R² = 35.7%; p = 0.001). Dung type and sex accounted for 15% (p = 0.001) and 10.1% (p = 0.001), respectively. The R² value indicated that species was the strongest factor in differentiating individuals, but at the same time, the significant p-value also found for dung type and sex showed that these two other factors also contribute to differentiating individuals even if to a smaller extent.

Conversely, for fungi, the ordination plot of the 90 individuals showed a clearer distribution of the points in the nMDS space when we grouped the individuals according to dung type ( Figure 1B ) instead of species ( Supplementary Figure S4B ). R² values confirmed these results as we found that the dung type showed the highest R² value (R² = 30.4%; p = 0.001). Species and sex accounted for 5.2% (p = 0.051) and 1.3% (p = 0.282), respectively. Furthermore, in this case, species was only nearly significant, and sex was not significant, so these two factors do not seem to contribute to the differentiation of individuals. However, we observed a clear separation of individuals collected from cow dung and instead a general overlap between individuals collected from sheep and donkey dung.

The ordinations in Figure 1 showed good diversification between samples with regard to prokaryotes when we grouped individuals by species, and with regard to fungi when we grouped individuals by dung type. Thus, using the same nMDS results, the distinction between individuals was even greater when we considered the effect of the most important factor only in each of the two cases ( Figure 2 ). When we separated the prokaryotes according to dung type, we noticed the effect of species in differentiating individuals. For example, OFT was well separated from the other species in all dung types ( Figure 2A ). In contrast, the distinction between OMD and OVT was less marked albeit still evident, especially in sheep dung, where the distinction between the three species was particularly good despite the presence of an OVT outlier. The effect of sex was mainly visible between individuals collected in donkey dung for all three species. We also observed an effect of sex in OFT in cow dung. However, the effect of sex in differentiating individuals was not appreciable in the other species or dung types. Considering, once more, the R² values, sex was the factor with the lowest value.

With regard to fungi, when we separated the samples according to species, we noticed the effect of dung type in differentiating individuals. Mainly, we observed a clear separation between cow and sheep dung in all species ( Figure 2B ). Donkey dung samples, on the other hand, were always slightly mixed with the other dung types, especially in OMD where there was a large overlap between donkey dung and cow and sheep dung. For fungi, the effect of sex was even less evident than what was observed for prokaryotes. The only group with good separation between males and females was the OMD collected in donkey dung. R² and p-values identified sex as a weak and non-significant factor.

We evaluated the quantity of exclusive and core ASVs by considering both i) the number of ASVs and ii) the number of reads per ASV, as relative percentages (i.e. the relative abundance of exclusive and core ASVs).

Overall, exclusive prokaryotic ASVs (i.e. those exclusive to each of the nine species-dung type pairs, Figure 3 , see petals) comprised 78.8% of the 4659 detected ASVs and 12.5% of the 1193 detected reads. Exclusive fungal ASVs comprised 82.2% of the 4549 detected ASVs and 30.3% of the 1901 detected reads. The number of prokaryotic and fungal ASVs which were exclusive in each species-dung type pair were similar (on average 408 for prokaryotes, and 415 for fungi, 8.8% and 9.1%, respectively), but the ASV relative abundances were significantly lower in prokaryotes (on average 1.4%) than in fungi (3.4%). Core ASVs (i.e. those shared between all nine species-dung type pairs) were very few: 30 in prokaryotes and only 5 in fungi, corresponding to 0.6 and 0.1% of the detected ASVs, respectively; however, the percentage of core ASV abundances were higher with regard to their number, corresponding to 14.7% and 2.1%, respectively.

We detected very few indicator ASVs (i.e. those significantly associated with species, dung type or pairs) ( Table 1 ), and ASVs associated with only one of the three dung beetle species were more numerous in prokaryotes than in fungi (on average 73.7 and 5.3, respectively). Vice-versa, ASVs associated with just one of the three dung types were more numerous in fungi than in prokaryotes (on average 45.7 and 24, respectively). We detected very few indicator ASVs across the nine species-dung type pairs in both prokaryotes and fungi, on average 7.4 and 5.7, respectively.

Among the prokaryotes, Proteobacteria, Bacteroidota and Firmicutes accounted for more than 80% of the ASVs ( Supplementary Figure S5A ), while Archaea were barely detectable. At the genus level ( Figure 4A ), considering all samples, the top three genera in terms of abundance were Apibacter, Enterococcus and Pseudomonas. However, the taxonomic composition of each pair proved to be relatively constant regardless of the dung type, whilst varying considerably from one host species to another. In particular, the only significant difference between dung types was found for Stenotrophomonas (F_2,6 = 20.6, p < 0.01), which was significantly more abundant in cow dung, whereas several genera showed significantly different abundances between the dung beetle species. For example, the relative abundances of Pasteurella varied between the three host-species (F_2,6 = 64.478, p < 0.001), with higher values in OFT than in OMD (17% ± 1% vs 4% ± 2%; p < 0.001) or OVT (17% ± 1% vs 4% ± 2%; p < 0.001). At the same time, OFT hosted significantly fewer Dysgonomonas than OMD (2% ± 0% vs 11% ± 3%; p < 0.05) or OVT (2% ± 0% vs 13% ± 6%; p < 0.05). Apibacter was less abundant in OVT than in OFT (4% ± 4% vs 12% ± 4%; p < 0.05) or OMD (4% ± 4% vs 19% ± 2%; p < 0.01). Acinetobacter was observed in all samples, with no significant differences detected between species (F_2,6 = 1.768, p = 0.249). Nevertheless, the highest percentage of abundance for Acinetobacter was observed in OFT collected from sheep dung (21%). The genus Wolbachia was found in all samples, with no significant differences between species (F_2,6 = 0.497, p = 0.632); however, Wolbachia did account for a considerable percentage in OVT collected from cow dung (19%).

With regards to fungi, Ascomycota and Basidiomycota were the most abundant phyla, accounting for more than 95% of the ASVs ( Supplementary Figure S5B ). At the genus level ( Figure 4B ), considering all samples, the top three genera in terms of abundance were Cladosporium, Alternaria, and Trichosporon. Nevertheless, the taxonomic composition of each pair was relatively constant regardless of beetle species, whilst varying from one dung type to another. In particular, we found significant differences in the abundance of Malassezia in relation to dung type (F_2,6 = 25.731, p < 0.01), with higher values in sheep dung than in cow dung (10% ± 2% vs 2% ± 2%; p < 0.01) and donkey dung (10% ± 2% vs 0% ± 1%; p < 0.01). Overall, while donkey dung and sheep dung samples showed a more constant internal set-up, we observed large internal differences in cow dung, even at the species level; for example, OFT was dominated by Cladosporium (11%), OMD showed a high abundance of Lithophyla and Bullera (21% and 19%, respectively), and OVT was dominated by Peniophora (18%). On the other hand, beetles collected from donkey dung exhibited high relative abundances of Apiotrichum and Fusarium, while Cladosporium and Alternaria were, on average, more abundant in beetles collected from sheep dung, although none of these differences were statistically significant (ANOVA).

For the sake of simplicity, all genera outside the ten most abundant have been pooled and labelled as ‘others’ in the bar charts ( Figure 4 ). For prokaryotes, ‘other’ genera accounted for about 25% of the total microbial genera in each pair, whereas the percentages were much higher for fungi (on average more than 50%). For example, in the OFT-cow pair, ‘other’ genera accounted for more than 75%.

The functional composition of prokaryotes and fungi was relatively constant in all samples regardless of the species and type of dung consumed.

More than 95% of the detected prokaryotes had metabolisms based on aerobic or anaerobic chemoheterotrophy ( Supplementary Figure S6A ). The three host species showed approximately the same functional set-up ( Figure 5A ). The only prokaryotic function with a significantly different distribution in abundance between host species was fumarate respiration (F_2,6 = 67.0, p < 0.001), which was more abundant in OFT than in OMD (5% ± 0% vs 2% ± 0%; p < 0.001) or OVT (5% ± 0% vs 3% ± 1%; p < 0.001). Fermentation, a typical carbon metabolic pathway, was the most common metabolic process in the prokaryotic microbiota, with an average of 37% ± 7% of reads, irrespective of the species considered. However, some differences in the abundances of fermentation function across dung types should be noted (F_2,6 = 5.378, p < 0.05). These weak significant differences mainly seemed to be related to cow dung for which fermentation was relatively low, whereas the nitrogen-related metabolic pathways were more prevalent (i.e. nitrate respiration: F_2,6 = 5.15, p < 0.05; nitrate reduction: F_2,6 = 7.389, p < 0.05; nitrogen respiration: F_2,6 = 5.15, p < 0.05). However, when examining the Tukey HDS results, all comparisons between dung types resulted non-significant, except for nitrate reduction, which was more abundant in beetles collected from cow dung than from sheep dung (9% ± 1% vs 4% ± 2%; p < 0.05). A relevant percentage in each pair consisted of prokaryotes identified as parasites or animal symbionts (on average 14% ± 2%, regardless of the species and dung type considered).

Regarding fungi, from a broader perspective, the most common trophic modes were saprotroph, followed by pathotroph and symbiotroph ( Supplementary Figure S6B ). Plant pathogens, animal pathogens and wood saprotrophs were the most common guild in the fungal microbiota, averaging 22% ± 1%, 22% ± 3% and 18% ± 2%, respectively, regardless of the host species and dung type considered ( Figure 5B ). No significant differences were found in relation to dung type or host species, with the exception of wood saprotrophic fungi (F_2,6 = 6.067, p < 0.05), which were slightly more abundant in OMD than in OVT (20% ± 2% vs 17% ± 1%, respectively; p < 0.05). The relative abundance of the dung saprotroph guild was equal in all pairs (1% ± 1% on average; no significant differences in relation to dung type or host species).