In May, live adult dung beetles (N = 88) belonging to nine species across four tribes and characterised by three nesting strategies (Table 1) were collected from dung deposited by a herd of Piedmontese cattle. The collection pastures are located in a managed submontane agroecosystem characterised by semi-natural grasslands in the western Italian Alps [Mocchie (TO), Susa Valley, Piedmont (45°8.16'N, 7°17.76'E)], within a radius of approximately 200 m of the specified GPS point and at an elevation of 800–900 m a.s.l. At least eight individuals per species were collected to capture intraspecific variability in the gut microbiota (Natta et al. 2024). The species and sex of the individuals were determined based on morphological features (Baraud 1992; Dellacasa and Dellacasa 2006; Ballerio et al. 2014), and species identification was then confirmed through molecular characterisation using the mitochondrial cytochrome c oxidase subunit I (COI) gene sequence as a DNA barcode (https://boldsystems.org/) as described in Natta et al. (2025a).

Dung samples were collected from two randomly selected dung pats and soil samples at a depth of 5 cm from around the same pats. These environmental samples were stored in 1.5 mL Eppendorf tubes and homogenised upon arrival at the laboratory, and absolute ethanol was added to the tubes containing dung to prevent fermentation. All samples were stored at –20 °C until analysis.

To empty the insects’ guts as much as possible, the beetles were kept in a terrarium without food for at least 24 h (Zhang et al. 2021; Chen et al. 2024). They were then euthanised by submersion in absolute ethanol and immediately dissected to extract the entire gut. All dissection tools were sterilised using a 30% sodium hypochlorite solution and then washed in distilled water before use. Each removed gut was preserved in absolute ethanol and stored at 4 °C.

DNA was extracted from gut and dung samples using the modified cetyltrimethylammonium bromide (CTAB) method described by Natta et al. (2024). DNA from soil samples was extracted using the DNeasy® PowerSoil® Pro Kit (QIAGEN).

A DNA metabarcoding approach was used to investigate the microbiota. For the prokaryotic component, the 16S ribosomal RNA (rRNA) gene was amplified using the primer set 515fB (5'–GTGYCAGCMGCCGCGGTAA–3') (Parada et al. 2016) and 806rB (5'–GGACTACNVGGGTWTCTAAT–3') (Apprill et al. 2015). The polymerase chain reaction (PCR) cycling programme consisted of an initial denaturation step at 94 °C for 3 min, followed by 35 cycles at 94 °C for 45 s, 55 °C for 60 s, and 72 °C for 90 s, with a final extension step at 72 °C for 10 min.

To characterise the total fungal community, the nuclear ribosomal internal transcribed spacer 2 (ITS2) region was amplified using a semi-nested PCR approach. In the first PCR, the entire ITS region (ITS1–5.8S–ITS2) was amplified using the universal fungal primer pair ITS1F (5'–CTTGGTCATTTAGAGGAAGTAA–3') (Gardes and Bruns 1993) and ITS4 (5'–TCCTCCGCTTATTGATATGC–3') (Ihrmark et al. 2012). The amplification protocol consisted of an initial denaturation step at 94 °C for 5 min, followed by 35 cycles at 94 °C for 30 s, 54 °C for 45 s, and 72 °C for 60 s, with a final extension at 72 °C for 10 min. PCR products were visualised on an agarose gel, diluted 1:20, and used as templates for the semi-nested PCR targeting the ITS2 region. This second amplification employed the primer pair fITS7 (5'–GTGARTCATCGAATCTTTG–3') and ITS4 (5'–TCCTCCGCTTATTGATATGC–3') (Ihrmark et al. 2012). The semi-nested PCR cycling conditions consisted of an initial denaturation at 94 °C for 5 min, followed by 35 cycles at 94 °C for 30 s, 50 °C for 30 s, and 72 °C for 30 s, followed by a final extension step at 72 °C for 7 min.

Illumina overhang adapter sequences were added to the primers used for library preparation. Specifically, Illumina overhang adapter sequences (forward overhang: 5'–TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG–[locus-specific target primer]; reverse overhang: 5'–GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG–[locus-specific target primer]) were added to primers 515fB and 806rB, which were used for bacterial 16S rRNA gene amplification, and to primers fITS7 and ITS4, which were used in the semi-nested PCR for fungal ITS2 amplification.

PCR reactions were performed in a final volume of 25 μL using 5 μL of 5× FIREPol® Master Mix (12.5 mM MgCl₂; Solis BioDyne), 0.2 μM of each primer, and 20 ng of genomic DNA. Extracted DNA was amplified in triplicate and pooled before purification. Purification and quantification of the amplified DNA were conducted following Natta et al. (2025a). The negative controls for all PCR reactions were pooled, purified, and quantified in the same manner as the samples for both prokaryotic and fungal analyses. Amplicons were subsequently sent to IGA Technologies (Udine, Italy) for Illumina NovaSeq sequencing (2 × 250 bp).

Bioinformatics and statistical analyses were conducted separately for prokaryotes and fungi. The sequences were analysed using QIIME2 (Quantitative Insights Into Microbial Ecology 2, v. 2023.5) microbiome bioinformatics platform (Bolyen et al. 2019) after the removal of adapters and primers. Denoising and quality control, including chimaera removal, were performed using the DADA2 plugin (Callahan et al. 2016) through the “qiime dada2 denoise-paired” command. Chimaera detection was carried out using the “consensus” method. The taxonomic assignment of the prokaryotic community was achieved using the Greengenes2 database (McDonald et al. 2024), whereas the UNITE database v. 19.02.2025 (Abarenkov et al. 2025) was used for fungi. Phylogenetic trees were generated using 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 the corresponding metadata files, were imported into RStudio (Posit team 2025) using the R package qiime2R v. 0.99.6 (Bisanz 2022) to generate phyloseq objects. Two phyloseq objects were constructed, corresponding to the prokaryotic and fungal datasets, respectively.

The amplicon sequence variants (ASVs) detected in the negative control were removed from the phyloseq objects before analysis (additional information on the ASVs present in the negative control can be found in Suppl. material 1: fig. S1). ASVs shared by both the beetles and the environmental samples (i.e. dung and soil samples) were subsequently identified by analysing the results of the flower plots generated from all 92 samples using the plot_venn function in the R package microeco (Liu et al. 2021). Environmental ASVs were removed from the phyloseq objects, and all subsequent diversity analyses were conducted on the resulting filtered datasets. However, to assess the influence of these environmental ASVs on beta diversity and functional profiles, the results of the analyses carried out by also considering these ASVs in the phyloseq objects are included in the Supplementary materials.

To enable comparisons between samples with non-uniform coverage, the ASV tables were normalised using the rarefy_even_depth function in the R package phyloseq v. 1.36.0 (McMurdie and Holmes 2013). Rarefaction curves of the rarefied ASV tables were obtained using the rarecurve function in the R package vegan v. 2.6-4 (Oksanen et al. 2022). Absolute counts (i.e. the number of reads) and relative abundances (i.e. the ratio of the number of reads belonging to an ASV in a sample to the total number of reads in that sample) were then calculated for each ASV.

Alpha diversity was assessed using the “Observed ASVs”, “Shannon”, and “Faith’s Phylogenetic Diversity (Faith’s PD)” indices, calculated with the estimate_richness function in the R package phyloseq (McMurdie and Holmes 2013). Pairwise relationships among diversity indices were evaluated using Spearman’s rank correlation coefficients, computed with the cor function in the R package corrplot (Wei et al. 2024). As all three indices were highly correlated, only one, namely Faith’s PD, was used. The effects of species, tribe, nesting strategy, and sex on alpha diversity were then tested using analysis of variance (ANOVA), and variance homogeneity was checked using Levene’s test. In cases of unequal variances, a linear model was fitted using generalised least squares with the gls function in the R package nlme (Pinheiro et al. 2025).

Dissimilarity among individuals in terms of both prokaryotes and fungi was visualised by means of distance-based redundancy analysis (dbRDA) based on a Bray–Curtis distance matrix of ASV composition, using the R package vegan. Before running the dbRDA, the Hellinger transformation was applied to the ASV table. The dbRDA results were plotted using the R package collection tidyverse (Wickham et al. 2019). An envfit analysis was performed on the dbRDA ordination to determine which factors had the greatest ability to differentiate among individuals based on the R² values and p-values. An ANOVA was then performed on the full model plotted in the dbRDA.

Flower plots (generated using a custom script available upon request) were obtained to visualise the shared (i.e. core) and exclusive ASVs associated with each species, tribe, and nesting strategy. To visualise the taxonomic diversity of the different samples, percentage-stacked bar charts were created using the R package collection tidyverse (Wickham et al. 2019). To facilitate their representation in the bar charts, only the 15 most abundant genera (i.e. those with the highest number of reads out of the total) were considered, and the other genera were grouped under “Others”.

The predicted functions of the identified prokaryotic and fungal communities were assessed using FAPROTAX (Functional Annotation of Prokaryotic Taxa) (Louca et al. 2016) and FUNGuild (Nguyen et al. 2016), respectively, both of which were implemented in the R package microeco (Liu et al. 2021). As with taxonomic diversity, percentage-stacked bar charts were used to visualise differences in trophic behaviour among the samples, with only the 15 most abundant functions represented. Differences among species or tribes in taxonomic and functional composition were tested using t-tests or chi-squared tests.

Bioinformatic analysis identified 12,501 prokaryotic ASVs at a rarefaction depth of 6,167 reads per sample and 3,301 fungal ASVs at a rarefaction depth of 3,326 reads per sample. Rarefaction curves based on observed ASVs indicated that all 88 beetle samples reached an asymptote for both prokaryotes and fungi, indicating that the sequencing depth was sufficient to capture gut microbiota diversity (Suppl. material 1: fig. S2A, B).

Concerning prokaryotes, significant differences in the Faith’s PD index were detected among species (F = 8.85, P < 0.001), tribes (F = 14.17, P < 0.001), and nesting strategies (F = 11.27, P < 0.001), whereas no significant differences were observed between sexes (F = 0.948, P = 0.341). Regarding fungi, significant differences in the Faith’s PD index were detected among tribes (F = 3.39, P < 0.05) and nesting strategies (F = 4.92, P < 0.01), but no significant differences were found among species (F = 1.87, P = 0.077) or between sexes (F = 1.06, P = 0.307).

The dbRDA ordination of prokaryotes showed that the gut communities of the sampled individuals differed consistently according to host species and tribe. The nine dung beetle species were well separated across the two-dimensional plot, forming four clusters of species that corresponded to the four tribes (Fig. 1). The four clusters did not overlap, whereas the convex hulls of the species within each cluster showed partial overlap but still maintained a clear degree of interspecific differentiation. Specifically, A. pedellus and C. erraticus (tribe Aphodiini) were grouped together and clearly separated from the other three clusters formed by A. stercorosus, G. mutator, and T. pyrenaeus (tribe Geotrupini); O. fracticornis, O. medius, and O. verticicornis (tribe Onthophagini); and S. shaefferi (tribe Sisyphini). The envfit analysis confirmed this result, as species produced a significantly higher R² value (R² = 88.1%; P = 0.001), followed by tribe (R² = 79.7%; P = 0.001). Nesting strategy accounted for 44.7% (P = 0.001), and sex, the only non-significant factor, accounted for 2.74% (P = 0.087). Intraspecific dispersion varied among species, with some species showing compact clusters of individuals and others showing dispersed clusters of individuals. The ANOVA performed on the full model plotted in the dbRDA yielded a statistically significant result (F = 3.07, P < 0.001). The comparative analysis conducted while preserving the environmental ASVs largely mirrored these results, both in terms of spatial ordination and the envfit R² values (Suppl. material 1: fig. S3).

Conversely, the dbRDA ordination of fungal communities revealed no consistent differences in gut community composition among host species or tribes. The nine dung beetle species were not well separated across the two-dimensional plot: seven were plotted in the upper-right corner and showed extensive overlap, one was plotted in the right portion of the graph, and one was plotted in the left portion (Fig. 2). The plot revealed no evident differentiation among the four tribes. Specifically, A. pedellus and A. stercorosus formed well-defined and largely separate clusters, indicating marked differences in fungal community composition compared with the other species. The other seven species occupied overlapping regions of the ordination space, reflecting similar gut fungal assemblages. Nevertheless, when only these seven species were considered, reduced overlap among the polygons was observed (Suppl. material 1: fig. S4). Indeed, on the first constrained axis (dbRDA1), C. erraticus and T. pyrenaeus were clearly distinguished from other species, such as G. mutator, S. shaefferi, O. fracticornis, and O. medius, whereas the second axis (dbRDA2) further distinguished among the species, particularly within the genus Onthophagus, with O. fracticornis and O. medius grouped together and separated from O. verticicornis at higher dbRDA2 values. The envfit analysis, performed on the ordination of all nine species, identified species as the factor with the highest explanatory power (R² = 80.6%; P = 0.001). Tribe and nesting strategy were both significant factors, accounting for 26.4% (P = 0.001) and 23.2% (P = 0.001), respectively. Sex was the factor with the lowest R² value and was also non-significant (R² = 0.40%; P = 0.725). Intraspecific dispersion varied among species, with one species showing a very compact cluster of individuals (A. pedellus), whereas the individuals of another species were highly dispersed (A. stercorosus). The ANOVA performed on the full model returned a statistically significant result (F = 1.58, P < 0.001). The comparative analysis conducted while preserving the environmental ASVs showed similar results, except for the position of A. stercorosus in the main dbRDA plot (Suppl. material 1: fig. S5). Indeed, when the environmental ASVs were retained, A. pedellus showed a marked distance from the other species in the plot, which otherwise remained highly overlapping. This considerable similarity in the results between the two approaches (i.e. with and without the environmental ASVs) was also evidenced by the analogous envfit R² values (Suppl. material 1: fig. S5). Furthermore, when A. stercorosus and A. pedellus were excluded from the analyses, the similarity in the results between the two approaches was even more pronounced (Suppl. material 1: fig. S6).

The extent of exclusive and core ASVs was assessed by considering both (i) the number of ASVs and (ii) the number of reads per ASV, expressed as relative percentages (i.e. the relative abundance of exclusive and core ASVs). Exclusive ASVs were those unique to each of the nine species, whereas core ASVs were those shared by all nine species.

Overall, exclusive prokaryotic ASVs (Fig. 3, petals) accounted for 10,677 of the 12,501 detected ASVs (85.4%) but represented only 21.3% of the total reads. Similarly, exclusive fungal ASVs represented 2,761 of the 3,301 detected ASVs (83.6%) but represented only 34.1% of the total reads. The mean number of ASVs exclusive to each species was significantly higher in prokaryotes than in fungi (paired t-test: t₍₈₎ = 3.55, P < 0.01). However, the mean relative abundance of exclusive ASVs was significantly lower in prokaryotes than in fungi (paired t-test: t₍₈₎ = –3.33, P < 0.05). In addition, most species-exclusive ASVs were detected in a single individual of the species considered (on average, 96.6 ± 1.7% and 96.8 ± 1.8% for prokaryotes and fungi, respectively). The remaining fraction of these exclusive ASVs detected in two or more individuals of the same species was, therefore, very small. In particular, for prokaryotes, only one ASV was found in nine individuals of A. stercorosus, whereas for fungi, only one ASV was found in four individuals of G. mutator, and two ASVs were found in four individuals of S. shaefferi (Suppl. material 1: fig. S7). Core ASVs (i.e. those shared by all nine species) were rare, with only 47 detected in prokaryotes and nine in fungi (Fig. 3, centre), corresponding to 0.4% and 0.3% of all ASVs, respectively. Despite their low numbers, core ASVs accounted for a disproportionately higher fraction of the total reads, representing 19.5% of prokaryotic and 5.7% of fungal ASV abundances. The three most abundant prokaryotic genera among the core ASVs were Vagococcus_B, Apibacter, and Acinetobacter (Suppl. material 1: fig. S8A). With regard to fungi, over 50% of the core ASVs were assigned to the genus Candida, followed by an unclassified Ascomycota and Lithohypha (Suppl. material 1: fig. S8B).

Among the prokaryotes, Proteobacteria, Firmicutes_D, and Bacteroidota accounted for approximately 80% of all ASVs, whereas archaea were detected at very low abundance. At the genus level (Fig. 4), the three most abundant genera across all samples were Wolbachia, Apibacter, and Vagococcus_B. However, the relative abundances of the 15 most prevalent prokaryotic genera varied considerably among the nine host species, indicating pronounced species-specific structuring of the gut microbiota. In A. pedellus, C. erraticus, G. mutator, and S. shaefferi, low-abundance taxa grouped as “Others” comprised a substantial fraction of the community, representing approximately half of the total relative abundance. In contrast, O. fracticornis, O. medius, O. verticicornis, and T. pyrenaeus exhibited microbial profiles with strong dominance by one or two genera, which reached relative abundances of approximately 30%. Apibacter was the dominant component of the microbiota in two Onthophagus species, namely O. fracticornis and O. medius, where it accounted for approximately one third of the total (on average, 30 ± 8%). Wolbachia was a prominent component of the gut microbiota in O. verticicornis (24%) and was also present in O. fracticornis, O. medius, A. stercorosus, and S. shaefferi, accounting for 14 ± 4% of the community, while remaining negligible in the other species. Enterococcus_G was particularly abundant in T. pyrenaeus, contributing a quarter of the total microbiota (25%), whereas its relative abundance remained below 5% in all other species except A. stercorosus (9%). Additional genera, such as Vagococcus_B and Faecalibacter, displayed pronounced interspecific variability, reaching 15–20% in some hosts but remaining at low relative abundances in others. G. mutator and T. pyrenaeus exhibited a high abundance of Vagococcus_B (16% each), whereas Faecalibacter occurred at a high percentage only in C. erraticus (18%). Overall, despite the presence of shared dominant genera across species, their relative contributions differed markedly, resulting in distinct gut microbiota profiles for each species. Nevertheless, a pattern emerged when the host species were clustered into tribes according to their phylogenetic relationships (Suppl. material 1: fig. S9A). Specifically, the tribe Aphodiini displayed the highest abundance of Acinetobacter (χ² = 11.67, P < 0.01), the tribe Onthophagini had the highest abundance of Apibacter (χ² = 68.86, P < 0.001), the tribe Geotrupini had the highest levels of Enterococcus_G (χ² = 24.20, P < 0.01), and Sisyphini, represented only by S. shaefferi, presented a combination of genera that distinguished it from the remaining tribes. These results were confirmed by the analysis of the core ASVs among species of the same tribe, in which Acinetobacter accounted for 42% of the core ASVs of Aphodiini, Apibacter accounted for 38% of the core ASVs of Onthophagini, and Enterococcus_G and Vagococcus_B accounted for 24% and 22% of the core ASVs of Geotrupini, respectively (Suppl. material 1: fig. S10).

With regard to fungi, Ascomycota and Basidiomycota were the most abundant phyla, accounting for 95% of the ASVs. At the genus level (Fig. 4), the three most abundant genera across all samples were Apiotrichum, Epicoccum, and an unclassified genus of Agaricomycetes. The composition of the fungal gut microbiota varied considerably among the nine species analysed. Clear interspecific differences were observed in both the relative abundance and dominance of the 15 most prevalent fungal genera. In some species, the fungal community was strongly dominated by a limited number of taxa, with Apiotrichum representing 38% of the total relative abundance in G. mutator and Trichosporon accounting for 32% in S. shaefferi. In contrast, O. fracticornis exhibited a more diverse community structure, with no individual genus exceeding 10% relative abundance and low-abundance taxa, grouped collectively as “Others”, accounting for approximately 70% of the total community. Genera such as Papiliotrema, Malassezia, and Cutaneotrichosporon were detected across multiple species, generally contributing between 5% and 20%, depending on the host. Other genera, including Epicoccum, Alternaria, and Lithohypha, reached moderate relative abundances in specific hosts, in some cases accounting for approximately 15–20% of the fungal community. For instance, Epicoccum and Alternaria were mostly found in C. erraticus (18% and 20%, respectively), whereas Lithohypha was most dominant in O. medius (17%). The dbRDA ordination plot suggested that A. stercorosus and A. pedellus hosted distinct fungal communities. The former was the only species to show a substantial presence of Candida (35%), whereas the latter hosted genera such as Hannaella and Nigrospora, included in the “Others” group, which were abundant only in this species while representing less than 1% in a few other species. Overall, the relative contributions of fungal genera varied among species, indicating pronounced species-specific differences in gut microbiota composition in a few cases. When the host species were grouped into tribes, the tribe Aphodiini displayed the highest abundance of Alternaria (χ² = 17.20, P < 0.001), the tribe Geotrupini had the highest levels of Candida (χ² = 27.50, P < 0.001), and the tribe Sisyphini, represented only by S. shaefferi, showed a significantly higher abundance of Trichosporon than the other tribes (χ² = 82.60, P < 0.001) (Suppl. material 1: fig. S9B). Moreover, analysis of only the core ASVs among species of the same tribe revealed notable details. For example, Aphodiini were dominated by Hannaella (34%), Sisyphini were dominated by Trichosporon (32%), and Geotrupini were dominated by Candida and an unclassified genus of Agaricomycetes (26% and 21%, respectively). By contrast, the tribe Onthophagini showed a mixture of genera, with no clearly dominant taxon (Suppl. material 1: fig. S11).

Despite interspecific differences in taxonomic composition, the functional composition of prokaryotes and fungi remained relatively consistent across all samples, regardless of species (Fig. 5).

Fermentation was identified as the predominant metabolic process among the prokaryotes hosted in the guts of the nine species, with an average relative abundance of 28 ± 9%. Nitrogen-related metabolic pathways (i.e. nitrate reduction, nitrate respiration, nitrite denitrification, nitrite respiration, nitrogen fixation, nitrogen respiration, and ureolysis) represented nearly one quarter of the total functional profile in each species and were relatively constant among species, with an average contribution of 23 ± 4%. The comparative analysis conducted while preserving the environmental ASVs did not reveal any substantial differences from the results obtained through dataset cleaning (Suppl. material 1: fig. S12). The three species belonging to the tribe Onthophagini (O. fracticornis, O. medius, and O. verticicornis) and the three species belonging to the tribe Geotrupini (A. stercorosus, G. mutator, and T. pyrenaeus) exhibited very similar functional profiles, in which fermentation was predominant, accounting for at least 30% of the total abundance. Despite this relatively constant pattern, slight differences were also detected among species. For instance, A. pedellus, C. erraticus, and S. shaefferi displayed lower levels of fermentation than the other species (Welch’s t-test: t_(3.10) = –6.68, P < 0.01; mean = 16.33 vs. 33.83) but higher levels of aromatic compound degradation (Welch’s t-test: t_(2.85) = 4.93, P < 0.05; mean = 12.67 vs. 6.17).

Regarding fungi, from a broader functional perspective, the most common trophic mode (i.e. the way by which fungi obtain nutrients) was saprotrophy, accounting for approximately half of the total relative abundance (on average, 46 ± 4%), followed by pathotrophy (35 ± 4%) and symbiotrophy (19 ± 1%). Among the saprotrophs, the dominant fungal guild was wood saprotrophs (15 ± 2%), followed by soil saprotrophs (5 ± 3%), litter saprotrophs (3 ± 2%), and plant-, leaf-, and dung-associated saprotrophs, which were present at low relative abundances (1 ± 1%). These proportions were largely consistent across species, except for O. verticicornis, A. stercorosus, and G. mutator, which exhibited higher relative abundances of soil saprotrophs (Welch’s t-test: t_(2.51) = 4.71, P < 0.05; mean = 8.67 vs. 2.67). Plant and animal pathogens were abundant in all species, with comparable relative contributions (20 ± 1% and 17 ± 5%, respectively). The comparative analysis conducted while preserving the environmental ASVs did not reveal any substantial differences from the results obtained through dataset cleaning (Suppl. material 1: fig. S13).

Approximately 12,000 prokaryotic ASVs and 3,000 fungal ASVs were detected, indicating that prokaryotic alpha diversity was four times greater than fungal diversity. This finding aligns with previous studies indicating higher prokaryotic diversity in the dung beetle gut microbiota (Ziganshina et al. 2018; Natta et al. 2025a). Alpha-diversity analyses also revealed contrasting patterns between prokaryotic and fungal communities. Prokaryotic phylogenetic diversity varied significantly among species, tribes, and nesting strategies, whereas fungal diversity was significantly associated only with tribe and nesting strategy. The lack of sex-related effects on alpha diversity observed here further corroborates earlier findings in dung beetles, suggesting that sex plays a marginal role in shaping gut microbial diversity compared with host phylogeny and diet (Suárez-Moo et al. 2020; Chen et al. 2024; Natta et al. 2025a).

Beta-diversity analyses further emphasised the contrasting patterns between prokaryotes and fungi. Regarding prokaryotes, the dbRDA results showed that species from the four tribes considered clearly hosted different bacterial communities and that, even within the same tribe, the species tended to maintain distinct microbiota, demonstrating that host-species taxonomy was reflected in their microbiota. The presence of different bacterial microbiota among dung beetle tribes suggests a stable host–microbe association, potentially mediated through long-term cooperative association. Comparable patterns have been described in dung beetles and other insects, in which prokaryotic communities show evidence of host-specific filtering and partial evolutionary conservation (Yun et al. 2014; Brooks et al. 2016; Parker et al. 2019; Ebert et al. 2021). For example, Ebert et al. found that differences in the bacterial microbiota of dung beetles from the genera Onthophagus and Cephalodesmius were associated not only with phylogeny but also with gut morphophysiology (Ebert et al. 2021). Similarly, two species of Onthophagus, namely O. gazella and O. sagittarius, were found to harbour species-specific bacterial communities despite sharing the same diet (Parker et al. 2019). Similar results were also found for three other phylogenetically close species of Onthophagus that nonetheless exhibited distinct bacterial microbiota (Natta et al. 2025a).

In contrast, fungal communities appeared to be weakly structured. The dbRDA results showed that only A. pedellus and A. stercorosus formed well-defined and largely non-overlapping clusters, whereas the other seven species occupied overlapping regions. Although species identity explained a considerable proportion of the variation observed in the composition of the gut fungal community, it should be noted that A. pedellus and A. stercorosus determined the ordering of points in the dbRDA space and occupied a significant portion of the space, pushing the other species into the remaining space and thus influencing the R² values. In fact, when these two species were excluded from the analysis, the remaining species showed reduced overlap among polygons. The fungal microbiota appeared to be less strongly associated with host phylogeny and nesting strategy than the prokaryotic microbiota, consistent with the view that fungi are more frequently acquired from the environment (Franzini et al. 2016; Natta et al. 2025a). Given the intimate association of dung beetles with microbially rich substrates such as dung and soil, these environments likely act as major sources of fungal propagules, which the beetles ingest during feeding, resulting in more variable and less host-constrained fungal assemblages. Moreover, the results obtained when environmental ASVs were retained also confirmed these drivers of variability in prokaryotes and fungi: for prokaryotes, no substantial differences were observed between analyses with and without environmental ASVs – suggesting that these communities were more closely linked to phylogeny than to the environment and diet – whereas for fungi, some differences emerged, particularly in A. stercorosus. This species is known to have a highly diversified diet (Weithmann et al. 2020; Arias-Álvarez and Piquero-Rodriguez 2025), which probably explains why its fungal microbiota may be more strongly influenced by the presence of environmental ASVs commonly found in soil and dung.

For both prokaryotic and fungal gut communities, dung beetle nesting strategies showed moderate and significant R² values. Nevertheless, it is hypothesised that these results merely reflect the fact that individuals with different nesting strategies also belong to different species and tribes, which, as discussed above, exert a much greater influence on gut communities.

Similar to the results for alpha diversity, sex did not play a significant role in differentiating individuals in the dbRDA ordination for either prokaryotes or fungi, in agreement with previous results for dung beetles (Natta et al. 2024).

The observed variation in microbiota composition among individuals of the same species further suggests differences in the degree of host control over microbial communities. Species that form compact clusters in the ordination space may possess more specialised or tightly regulated gut microenvironments, whereas broader clusters may reflect greater dietary breadth, as previously hypothesised for beetles with diverse ecological niches (Franzini et al. 2016; Kolasa et al. 2019; de Beer et al. 2026).

The analysis of exclusive and core ASVs revealed a highly diverse and largely species-specific microbial assemblage. The majority of the detected ASVs were, in fact, exclusive to individual species, particularly among prokaryotes, yet these taxa accounted for only a minor fraction of the total reads. Similar patterns of high taxonomic turnover have previously been documented, suggesting that rare, host-specific microbial taxa contribute disproportionately to interspecific differentiation without dominating the overall community structure (Natta et al. 2025a). Moreover, most species-exclusive ASVs were detected in only a single individual of the species in question, casting doubt on the true biological exclusivity of these ASVs. In the case of ASVs represented by one or very few reads, their presence in only a single individual is inevitable, representing a form of “statistical exclusivity”. This high level of intraspecific variability, which is also supported by previous findings (Natta et al. 2024), serves as a reminder that it is advisable to consider at least 8–10 individuals per species when conducting analyses of this kind.

In contrast, the core microbiota shared across all nine species was numerically scarce but disproportionately abundant, especially among prokaryotes. This pattern closely resembles findings from previous work on dung beetles and other detritivorous insects, in which a limited set of widespread taxa dominated gut communities, in accordance with the minimal core model (Franzini et al. 2016; Natta et al. 2025a). The higher relative abundance of core ASVs in prokaryotes compared with fungi supports the idea that bacterial communities are more conserved across dung beetles, whereas fungal assemblages are more variable and potentially transient. These differences are consistent with the interpretation that bacteria form relatively stable mutualistic associations involved in digestion and nutrient cycling, whereas fungi are more often ingested from dung and soil and may act as opportunistic or short-lived gut associates (Estes et al. 2013; Schwab et al. 2016; Kucuk 2020; Šigut et al. 2022).

The prokaryotic gut microbiota was dominated by Proteobacteria, Firmicutes_D, and Bacteroidota, a pattern widely reported in dung beetles and other detritivorous beetles (Franzini et al. 2016; Kolasa et al. 2019; Suárez-Moo et al. 2020; Ebert et al. 2021; Chen et al. 2024). Previous studies have interpreted their dominance as evidence of adaptation to protein-rich and high-fibre diets such as dung (Chen et al. 2024; Jones et al. 2025). At the genus level, pronounced interspecific differences were observed despite the recurrence of a limited set of dominant taxa across hosts, as emphasised by the analysis of the core microbiota. Apibacter has been reported as a frequent and sometimes dominant component of the gut microbiota in Onthophagus beetles and other insects, and it has been hypothesised to play roles in carbohydrate metabolism and host nutrition (Kwong et al. 2018; Natta et al. 2025a). Acinetobacter, known to be a highly flexible decomposer of organic compounds, including dung, is also common in the microbiota of dung beetles (Kolasa et al. 2019; Jácome-Hernández et al. 2023; Chen et al. 2024; Natta et al. 2024, 2025a). Similarly, lactic acid bacteria such as Enterococcus and Vagococcus, which are also abundant in Thorectes lusitanicus (Geotrupini), are predominantly saccharolytic (i.e. they improve energy recovery by breaking down otherwise indigestible sugars), generate an acidic environment, and protect against putrefactive and sulphate-reducing bacteria that produce gas (Hernández et al. 2015).

The frequent detection of Wolbachia in several species also warrants consideration. Although Wolbachia is typically regarded as an intracellular endosymbiont rather than a gut-associated bacterium, its presence has been repeatedly reported in amplicon-based surveys of dung beetles and other insects (Hilgenboecker et al. 2008; Kolasa et al. 2019; Parker et al. 2020; Ebert et al. 2021). The effects of Wolbachia infections on host insects can be both beneficial, such as nutrient supplementation and protection against viruses, and deleterious, including feminisation and the killing of males (Parker et al. 2020).

In prokaryotes, a high degree of correspondence was observed between the dbRDA and taxonomic composition analyses, with both indicating that different dung beetle species harboured different microbial communities. By contrast, the results for fungi showed less concordance between the two approaches because, as mentioned above, the presence of two dung beetle species with very different microbiota reduced the dbRDA differences among the other species. This highlights the usefulness of the taxonomic composition results. In fact, despite the limitation that many groups of microorganisms were grouped into a generic “Others” category, the taxonomic composition analysis indicated that even species that completely overlapped in spatial ordination could still harbour partially distinct fungal communities.

Fungal communities were dominated by Ascomycota and Basidiomycota, in agreement with earlier work on dung beetles (Franzini et al. 2016; Chen et al. 2024; Natta et al. 2025a). Fungal assemblages showed weaker correspondence with host phylogeny than bacterial communities, as previously reported, but noteworthy differences were nevertheless highlighted. For example, the pronounced presence of Candida in A. stercorosus, coupled with its absence in the other species, likely explains the isolated position of this species in the two-dimensional dbRDA space. Candida yeasts are common inhabitants of beetle digestive tracts and are frequently reported in plant- or fungus-feeding beetles and dung beetles (Suh et al. 2008; Nwaefuna et al. 2023). They play important roles in nutrient breakdown and microbial community dynamics and, in some cases, form host-specific associations (Suh et al. 2008). Similarly, the isolated position of A. pedellus in the dbRDA space can be explained by the fact that this species harbours a heterogeneous assemblage of fungal taxa, including genera such as Hannaella and Nigrospora, which are almost exclusive to A. pedellus. Hannaella yeasts are widely distributed on the leaf surfaces of various plants (Li et al. 2021), and some species are believed to play an important role in promoting plant growth and biological control (Ji et al. 2024). Nigrospora, like many other fungi found in the guts of dung beetles, is a known plant pathogen that primarily infects leaves (Liu et al. 2024). Previous studies have reported Nigrospora to be a frequent component of the diets of various insects and to play a role in cellulose degradation (Zhang et al. 2018; Olivier et al. 2026). Hypothetically, Nigrospora may be acquired from plants consumed by cattle, subsequently excreted in dung, and finally ingested by individuals of A. pedellus. Following this hypothesis, it can be speculated that the dung beetle food chain may contribute to controlling these pathogenic fungi, promoting their removal from dung and potentially limiting their spread to other plants. The case of A. pedellus highlights the fact that two dung beetle species may have similar microbiota in terms of the 15 most abundant genera while differing greatly in the composition of the genera grouped into the “Others” category. In fact, the presence of high percentages of “Others” could limit the usefulness of these analyses unless a species-by-species study of the most abundant taxa in the microbiota is conducted, as in this case.

The absence of a clear tribal signal suggests that fungi are less strongly influenced by host evolutionary history and are instead shaped mostly by environmental acquisition, i.e. from dung, soil, or plant material. Nevertheless, the recurrent presence of genera such as Apiotrichum, Trichosporon, and Epicoccum across multiple species indicates that some fungi may be frequently ingested and are able to persist, at least transiently, in the gut environment. If gut fungi are acquired through the diet, differences in trophic preferences among dung beetle species (Urrutia et al. 2022; Rolando et al. in preparation) would explain how a wide range of beneficial fungal taxa is introduced into their gut microbiota. For example, fungal genera such as Apiotrichum and Trichosporon are known to be capable of assimilating or degrading lignocellulose (Akamine and Mishima 2023).

Despite the taxonomic differences among host species, the predicted functional profiles of both prokaryotic and fungal communities remained remarkably stable. This finding, although based on predicted functional profiles rather than measurements obtained through metatranscriptomics or metabolomics, corroborates previous studies reporting a decoupling between taxonomic composition and functional potential (i.e. functional redundancy) in dung beetles and other insects. These findings led to the hypothesis that microbial consortia could be assembled from different taxonomic pools yet still maintain functional outcomes (Moya and Ferrer 2016; Natta et al. 2025a). In prokaryotes, fermentation emerged as the dominant metabolic process across all species, underscoring the central role of gut bacteria in mediating the breakdown of complex organic substrates present in dung. The relative stability of nitrogen-related pathways among species further supports a conserved contribution of gut microbiota to nitrogen transformation and recycling. Similarly, fungal functional guilds were largely conserved among species, with saprotrophy being the most dominant trophic mode in all hosts. This pattern aligns with previous observations that dung beetle-associated fungi largely reflect saprotrophic lifestyles linked to dung and soil environments (Akamine and Mishima 2023; Natta et al. 2024, 2025a).