Seasonal Variation in the Organization of Dung Beetle Communities in the Moroccan Middle Atlas (Coleoptera: Scarabaeoidea)

: Dung beetles feed on and bury animal droppings, and their role is crucial in reducing the accumulation of manure, which diminishes the useful surface area of pastures. The aim of this research was to characterize the seasonal organization of dung beetle communities (Coleoptera: Scarabaeoidea) in the Middle Atlas region of Morocco in terms of core and satellite species. The beetles were collected using standard dung-baited traps. Four sites along a gradient of elevation were surveyed for one year every 7 to 10 days, depending on the season and local weather conditions. A total of 24,397 beetles were collected, belonging to 51 species. In most dung beetle communities, two to three species were found to be predominant, representing between 70 and 95% of all the individuals active at the same time but constituting only 10 to 30% of species diversity. The rapid succession of species at the same site limits the competition between species, allowing for the efﬁcient use of available trophic resources.


Introduction
Dung beetles (a group of Coleoptera that taxonomically includes the Scarabaeidae Scarabaeinae, most Scarabaeidae Aphodiinae, and most Geotrupidae) are among the insects primarily responsible for dung removal from the ground surface in many ecosystems.Factors affecting their abundance and distribution thus impact dung degradation [1][2][3][4].Their activity directly improves the fertility and physical characteristics of soils [5], helps to disperse seeds contained in dung [6], contributes to the reduction of dung breeding flies and gastrointestinal parasites of livestock [1,7], and in part reduces greenhouse gas emissions [8][9][10].Three main guilds have been identified to describe the feeding and/or nesting behavior of dung beetles: non-nesters (dwellers), paracoprids (tunnelers), and telecoprids (rollers) [11,12].These groups describe how species handle and consume dung, which determines where the dung is moved after it has been deposited.Non-nesters lay their eggs inside droppings, paracoprids dig nests below the droppings, and telecoprids cut off a parcel of dung and bury it at a distance.
Dung beetles are sensitive to habitat perturbations and are easily studied, making them an ideal focal taxon to evaluate the effects of anthropogenic disturbances [13][14][15].
Ecological functions are directly influenced by both community attributes and environmental variables and confirm the link between biodiversity, environment, and ecosystem functioning [16,17].Dung beetles are particularly vulnerable to the destruction of their habitat and the use of pesticides, especially veterinary medicines found in animal feces and resulting from the treatment of animals [18,19].About 644 dung beetle species are present in countries around the Mediterranean Basin [20], where the rate of endemism is very high [21,22].Studies on dung beetles in this biogeographical sub-region mostly focused on the spatio-temporal organization of their communities [23][24][25] but have, for the most part, not considered the combination of abundance and biomass of the species, which is a key factor in understanding the rapidity and level of degradation of animal dung.
The efficiency of dung beetles in using animal dung depends on both their abundance and size [26], the latter being correlated with their individual biomass [27].To reflect the importance of the species that compose a dung beetle community at a given time, some authors have proposed to divide them into the following three categories based on their abundance (number of individuals) and biomass: (1) core species (those representing at least 10% of the community, in both abundance and biomass [28,29]), (2) satellite species (those representing at least 10% of the community, in either abundance or biomass), and (3) accessory species (those that represent less than 10% of the community in both abundance and biomass) [30].Core and satellite species define the functionally most important species present at a given time within a community [31].Competition for trophic resources and reproduction sites in a dung beetle community can be mitigated if coexisting species of the same size belong to different guilds or if coexisting species of the same guild differ in size [28,32,33].In most situations, the core and satellite species include a few species whose cumulative frequency often ranges between 70 and 85% [32].
Janati-Idrissi [34] was the first author who investigated the spatio-temporal variation of core and satellite species in the dung beetle communities in Morocco.The present study extends this approach by examining the spatial and temporal structure of dung beetle communities in the Middle Atlas region of Morocco.The aim of this work is to go further by examining how diversity in dung beetle communities varied along an altitudinal gradient over an annual cycle, environmental variables contributing to the structural composition of communities.

Study Sites
The Middle Atlas is a mountain system that extends for some 350 km from the southwest to the northeast of Morocco, between the Rif and the High Atlas.The western side usually receives between 1000 and 1500 mm of rain a year, unlike the eastern valleys, which are far drier.More continental than the Rif and generally wetter than the High Atlas, the Middle Atlas experiences cold winters, with snowfalls above 2000 to 2500 m between November and April.Pastoralism is the main economic resource of this region.Dung beetle communities were studied in four stations located between the Sais Plain and the Middle Atlas: one close to Fez-Sais, one close to the town of Immouzzer, and two (Ifrane I and Ifrane II) near the town of Ifrane.These sites differed from each other in several factors (bioclimatic zones, altitude, soil, and vegetation cover) (Figure 1).
The Fez-Sais station (33 • 54 N-4 • 59 W; elevation 609 m) is located at the edge of the Sais plain and before the Middle Atlas plateau.It is a karstic formation of lacustrine origin composed mainly of massive limestone and travertine associated with conglomerates and marl, resulting in a hard and rocky soil.The site is situated in the semi-arid bioclimatic zone with temperate winters and consists of an open environment with a predominance of herbaceous vegetation, mainly grazed by sheep.The site has been degraded by the expansion of urbanization in the Fez conurbation, which is increasingly restricting the grazing area.

Sampling
From September 2018 to August 2019, the study sites were surveyed every 7 to 10 days, depending on the season and local weather conditions.Three CSR-type baited traps [35] were set up simultaneously at each of the study sites.The traps, 15 cm deep and 30 cm in diameter, were buried at ground level about 50 m apart and covered with a largemesh (1 × 1 cm) metal grid supporting about 250 g of fresh cattle dung.The dung beetles attracted fell into the container, which had been three-quarters filled with a preservative liquid (solution of 80% water and 20% of ethyl alcohol (95%), and a few drops of liquid detergent).The traps were collected one week after they had been installed.

Identification of Dung Beetles
The collected beetles were preserved in a 70% ethyl alcohol solution and identified to the species level using a Zeiss Stemi 508 stereomicroscope and counted.The species identification was mainly based on the reference work on the Scarabaeoidea of North

Sampling
From September 2018 to August 2019, the study sites were surveyed every 7 to 10 days, depending on the season and local weather conditions.Three CSR-type baited traps [35] were set up simultaneously at each of the study sites.The traps, 15 cm deep and 30 cm in diameter, were buried at ground level about 50 m apart and covered with a large-mesh (1 × 1 cm) metal grid supporting about 250 g of fresh cattle dung.The dung beetles attracted fell into the container, which had been three-quarters filled with a preservative liquid (solution of 80% water and 20% of ethyl alcohol (95%), and a few drops of liquid detergent).The traps were collected one week after they had been installed.

Identification of Dung Beetles
The collected beetles were preserved in a 70% ethyl alcohol solution and identified to the species level using a Zeiss Stemi 508 stereomicroscope and counted.The species identification was mainly based on the reference work on the Scarabaeoidea of North Africa [36].For each species, 10 to 200 individuals (depending on their size) were airdried for 24 to 48 h, then placed in a fan-assisted oven heated to 70 • C for 12 h and weighed using a high-precision balance (APX-200, Denver Instrument GmbH, Göttingen, Germany; precision 0.1 mg).The mean individual biomass (dry weight) was calculated for each species.The total biomass of each species was calculated using the total number of individuals collected at a sampling event multiplied by the mean individual biomass for that species.A matrix of 51 species and 48 samples was drawn up from these surveys (51 species/48 samples [= 4 sites × 12 months]).

Statistical Analysis
The analyses of the data were carried out using the software R in version 4.1.3,under the Ade4, vegan, leaflet, labdsv, ggplot2, mvpart, MASS, and agricolae packages.We also used some additional scripts from the book Numerical Ecology with R [37,38].
Species were classified as core species (>10% total number and biomass), satellite species (>10% total number or biomass), or accessory (<10% total number and biomass).Core and satellite species were used in the analyses as these species would likely play a significant role in the recycling of dung.Cumulative frequencies were calculated for the total number and biomass for each month.A distribution map was built by processing all the gathered data in ArcGIS software.
The distribution of diversity requires an additive measurement of the variance within the sample, and linear diversity indices such as Hill numbers [39] facilitate a clearer interpretation of the obtained results.Many community ecologists currently prefer to use Hill numbers for taxonomic diversity and Hill ratios for regularity rather than Shannon entropy and Pielou regularity, respectively [40].This is because these Hill numbers, also known as "equivalent numbers" are simpler to understand, representing "the number of equally probable elements (individuals, species, etc.) needed to produce the observed value of the diversity index" [40].Hill's numbers [41] differ only in the importance assigned to relatively rarer species [39].
In most cases, N 1 (where N 1 = exp(H) [H: Shannon diversity (base e)]) and N 2 (where N 2 = (1/λ) [λ: Simpson diversity], the inverse of Simpson's index), are sufficient for answering any questions that a heterogeneity index can answer [42].N 1 offers greater importance to the evenness of the rarest species compared to N 2 .N 0 , which represents the number of species (i.e., species richness), is a diversity index that completely ignores the distribution of individuals among species [39].Note that the N 2 index, calculated as 1/λ, reduces sensitivity to variation in the abundance of highly abundant (generally few in numbers) species and hence indicates the number of dominant species present in a sample [40].
Hill's numbers (N), all measured in equivalent units of species, were recorded for each of the four collection sites over the course of the twelve-month sampling period.
Hill ratios (E) are calculated from Hill numbers and are denoted by E 1 , E 2 [40], E 1.0 [43], and E 2.0 [38], respectively.E 1 = N 1 /N 0 represents Shannon's regularity (Hill ratio), while E 2 = N 2 /N 0 represents Simpson's regularity (Hill ratio).They measure the regularity of species within a sample.When the species in the sample have the same abundance or dominance, or when there is only one species present, these two ratios are equal to 1 [40].It is important to note that regularity cannot be measured in a singular way, and its imprecision requires caution when applied [43].
To compare diversity values between sites, we used ANOVA followed by Bonferroni tests if conditions of normality and/or homoscedasticity were satisfied; otherwise, we used the Kruskal Wallis test followed by Wilcoxon tests.
The analysis of community composition was made using a redundancy analysis (RA), which combines regression and principal component analysis (PCA).It is a direct extension of regression analysis to model multivariate response data.Redundancy analysis is a two-step process [38].The first step is a multiple regression, where each object in Y is regressed on the explanatory variables in X, resulting in a matrix of fitted values Yfit.This step is calculated using the following linear equation: Yfit = X[X X] − 1X Y.In the second step, we applied a principal component analysis (PCA) to the fitted Yfit matrix to reduce dimensionality, i.e., to obtain the eigenvalues and eigenvectors.We then obtained a Z matrix containing the canonical axes, which correspond to linear combinations of the explanatory variables in the space of X.In the analysis of community composition, these canonical axes are interpreted as complex environmental gradients [38].
Before performing the redundancy analysis, the environmental or explanatory data were standardized (normalized).Moreover, to resolve the issue of double zeros in the matrix, we used the Hellinger transformation.To identify the significant explanatory variables, a progressive selection was conducted [38], with 1000 permutations to assess the significance of our ultimate model.
Relationships between samples were explored using a non-metric multidimensional scaling (NMDS) with the Bray-Curtis dissimilarity coefficient, which is a variation of Sørensen's index that accounts for species abundance and is unaffected by double absences.NMDS uses an iterative algorithm to locate items in the specified number of dimensions with the aim of minimizing a stress function, which ranges from 0 to 1, and evaluates the adjustment of the distance between items in the ordination space.The greater the precision of object representation in the ordination space, the lower the stress value.Object ordinations can be generated from any dissimilarity through the use of NMDS [38,40].
Finally, a correspondence analysis (CA) was used to create a temporal biotypology.This method has proven particularly effective as an ordination method [38,40].Its use in the field of biological zonation is justified in that its biocenoses are generally organized in a continuum.To homogenize the variances and to minimize the influence of the most abundant taxa, all our data underwent a logarithmic transformation: y = log(x + 1).

Environmental Data
Monthly meteorological data, including the mean, maximum, and minimum temperature, average amount of rainfall, relative humidity, and wind speed, were provided by the Sebou water basin agency (www.abhsebou.ma;accessed on 4 June 2023)).
Additionally, we used a GPS device (Garmin eTrex 10) to determine elevation and geographic coordinates.To identify Emberger's bioclimatic levels, we consulted the map produced by Brignon and Sauvage [44] at the Rabat Scientific Institute (formerly known as the Institut Scientifique Chérifien) (Table S6).

Faunal Composition of Assemblages and Organization of the Dung Beetle Communities
During the study, 24,397 individuals and 51 species were collected in total, including 27 Scarabaeinae species, 21 Aphodiinae species, and 3 Geotrupidae species (Table 1).In all the monitored sites, the core and satellite species (usually two to five) represented more than 70% and often more than 80% of the total number and/or biomass of the beetles in a community at any time of the year.Each of these was rarely dominant for more than a few weeks or months, being followed by another combination of core and satellite species, enabling almost 50% of the total number of species recorded during the year to acquire a core or satellite species status: 19 species (9 core and 10 satellite species) out of a total of 39 at Fez-Sais (Figure 2); 20 (13 core and 7 satellite species) out of 35 at Immouzzer (Figure 3); 17 (10 core and 7 satellite) out of 42 at Ifrane I (Figure 4); and 18 (12 core and 6 satellite) out of 43 at Ifrane II (Figure 5).
In autumn and winter, some non-nester species became dominant in the dung beetle communities at all sites.At Fez-Sais in September and October, Anomius castaneus dominated the community in terms of numbers (97.9% of individuals in September).In September, two species dominated, Cheironitis ungaricus and A. castaneus, which together accounted for almost 87% of the cumulative biomass of all the individuals in the community (Figure 2).At Immouzzer, A. castaneus and Aphodius fimetarius were seasonal species that appeared in large numbers between September and November.The winter non-nester species Amidorus cribricollis was also well represented between December and February (Figure 3).At Ifrane I, Nimbus affinis dominated the community throughout winter, both in numbers and biomass, in association with Thorectes trituberculatus in November and December.Nimbus affinis was the only active species in January (100% both in numbers and biomass) (Figure 4).At Ifrane II, during the winter season and from November to January, N. affinis was the dominant species both in terms of numbers and biomass.To a lesser extent, this species was accompanied by several satellite species, based either on the biomass they represented in the community (Scarabaeus sacer, Scarabaeus laticollis, and Onthophagus marginalis ssp.andalusicus) or the number of their individuals (A.cribricollis and A. fimetarius).
In all the monitored sites, the core and satellite species (usually two to fiv represented more than 70% and often more than 80% of the total number and/or bioma of the beetles in a community at any time of the year.Each of these was rarely domina for more than a few weeks or months, being followed by another combination of core an satellite species, enabling almost 50% of the total number of species recorded during th year to acquire a core or satellite species status: 19 species (9 core and 10 satellite specie out of a total of 39 at Fez-Sais (Figure 2); 20 (13 core and 7 satellite species) out of 35 Immouzzer (Figure 3); 17 (10 core and 7 satellite) out of 42 at Ifrane I (Figure 4); and 18 (1 core and 6 satellite) out of 43 at Ifrane II (Figure 5).satellite species are underlined.Squares and diamonds correspond to the relative abundanc biomass, respectively, of species that were part of the monthly groups of core species.The s the squares and diamonds was adjusted to the relative part of each species.In autumn and winter, some non-nester species became dominant in the dung b communities at all sites.At Fez-Sais in September and October, Anomius cast dominated the community in terms of numbers (97.9% of individuals in Septembe September, two species dominated, Cheironitis ungaricus and A. castaneus, which tog accounted for almost 87% of the cumulative biomass of all the individuals in community (Figure 2).At Immouzzer, A. castaneus and Aphodius fimetarius were sea species that appeared in large numbers between September and November.The w non-nester species Amidorus cribricollis was also well represented between Decembe February (Figure 3).At Ifrane I, Nimbus affinis dominated the community throug winter, both in numbers and biomass, in association with Thorectes trituberculat November and December.Nimbus affinis was the only active species in January (100% in numbers and biomass) (Figure 4).At Ifrane II, during the winter season and November to January, N. affinis was the dominant species both in terms of number biomass.To a lesser extent, this species was accompanied by several satellite sp based either on the biomass they represented in the community (Scarabaeus Scarabaeus laticollis, and Onthophagus marginalis ssp.andalusicus) or the number of individuals (A.cribricollis and A. fimetarius).
In spring and summer (March to August), the non-nester species disappeared were replaced by paracoprids and telecoprids.At Fez-Sais, three species were dom in May in terms of numbers: Euonthophagus crocatus (paracoprid), Gymnopleurus flage satellite species are underlined.Squares and diamonds correspond to the relative abundance and biomass, respectively, of species that were part of the monthly groups of core species.The size of the squares and diamonds was adjusted to the relative part of each species.

Core Species and Satellite Species
At all the sites, non-nester species became dominant in autumn, and this trend was confirmed in winter due to their numbers and biomass.At Fez-Sais, several small paracoprids replaced the non-nester species in early spring, especially Onthophagus melitaeus and E. crocatus, which were in turn supplanted from May until August by several telecoprids, especially G. flagellatus and more notably G. sturmi.At Immouzzer, the winter non-nester species were replaced by a cascading succession of small paracoprids from March to August, followed by large paracoprids in July and August (Cheironitis furcifer and C. irroratus).At Ifrane I and Ifrane II, the winter non-nester species were replaced by a combination of several small paracoprids (O.vacca and O. m. andalusicus) and telecoprids (S. sacer and S. schaefferi).

Diversity in the Communities and Environmental Variables
Calculations of Hill's diversity numbers N 0 , N 1, and N 2 and Hill's ratios E 1 and E 2 showed that all six indices change over time (Figure 6).The N 0 index expresses taxonomic richness, whereas the E 1 and E 2 indices (Hill ratios) measure the diversity and equitability of samples, respectively.

Diversity in the Communities and Environmental Variables
Calculations of Hill's diversity numbers N0, N1, and N2 and Hill's ratios E1 and E2 showed that all six indices change over time (Figure 6).The N0 index expresses taxonomic richness, whereas the E1 and E2 indices (Hill ratios) measure the diversity and equitability of samples, respectively.Figure 6 provides a clear visual representation of how abundance, diversity indices (N 0 , N 1 , and N 2 ), and regularity (E 1 and E 2 ) changed over different sampling sites and months, both spatially and temporally.For a more detailed analysis of the statistical significance of these spatial changes, see Figures 7 and 8.
There was a high abundance of individuals at higher elevations, particularly at Ifrane II in September, as well as in spring and early summer at Ifrane I and Ifrane II.For Ifrane I, the N 1 and N 2 values in March were particularly high (equaling 12 and 9, respectively).The difference in these diversity indices between the four sites was not statistically significant, except for N 0 (richness), which showed a clear difference between Ifrane II and (Fez-Sais) and (Immouzzer) (Figure 8).
Figure 9 is a summary of species frequency for each site over the twelve-month sampling period.It illustrates fluctuations in dung beetle communities between sites, as well as species that are sensitive or insensitive to seasonal changes.Species activity was highly variable from one month to the next.This is expressed by the presence or absence of some species in the surveys.Some species can be considered indifferent to seasonal fluctuations (e.g., Onthophagus opacicollis, present 11 times over 12 months at Ifrane II), while others were rare and absent from most sites, such as Plagiogonus esymoides (only present at Fez-Sais), which can be regarded as quite occasional (Figure 9).
The results of the full redundancy analysis (RA) suggested that the ten environmental explanatory variables (Table S6) retained in this study explained 41.05% (38.44% after adjustment) of the variance in the temporal and spatial structure of the dung beetle communities studied.Out of all the variables scrutinized, only four proved statistically relevant, namely temperature, rainfall, relative wind speed, and elevation.The final model of the canonical redundancy analysis (RDA) indicates that the chosen four environmental factors account for 38.86% (35.17% after adjustment) of the variance in the spatio-temporal structure of the dung beetle community examined.
significance of these spatial changes, see Figures 7 and 8.
There was a high abundance of individuals at higher elevations, particularly at Ifrane II in September, as well as in spring and early summer at Ifrane I and Ifrane II.For Ifrane I, the N1 and N2 values in March were particularly high (equaling 12 and 9, respectively).The difference in these diversity indices between the four sites was not statistically significant, except for N0 (richness), which showed a clear difference between Ifrane II and (Fez-Sais) and (Immouzzer) (Figure 8). Figure 9 is a summary of species frequency for each site over the twelve-month sampling period.It illustrates fluctuations in dung beetle communities between sites, as well as species that are sensitive or insensitive to seasonal changes.Species activity was highly variable from one month to the next.This is expressed by the presence or absence of some species in the surveys.Some species can be considered indifferent to seasonal fluctuations (e.g., Onthophagus opacicollis, present 11 times over 12 months at Ifrane II), while others were rare and absent from most sites, such as Plagiogonus esymoides (only present at Fez-Sais), which can be regarded as quite occasional (Figure 9).The results of the full redundancy analysis (RA) suggested that the ten environmental explanatory variables (Table S6) retained in this study explained 41.05% (38.44% after adjustment) of the variance in the temporal and spatial structure of the dung beetle communities studied.Out of all the variables scrutinized, only four proved statistically relevant, namely temperature, rainfall, relative wind speed, and elevation.The final model of the canonical redundancy analysis (RDA) indicates that the chosen four environmental factors account for 38.86% (35.17% after adjustment) of the variance in the spatio-temporal structure of the dung beetle community examined.The results of the full redundancy analysis (RA) suggested that the ten environmental explanatory variables (Table S6) retained in this study explained 41.05% (38.44% after adjustment) of the variance in the temporal and spatial structure of the dung beetle communities studied.Out of all the variables scrutinized, only four proved statistically relevant, namely temperature, rainfall, relative wind speed, and elevation.The final model of the canonical redundancy analysis (RDA) indicates that the chosen four environmental factors account for 38.86% (35.17% after adjustment) of the variance in the spatio-temporal structure of the dung beetle community examined.The findings indicate that the ultimate model is statistically significant (p = 0.0009), and each variable incorporated in this model is also significant (p = 0.0001).In addition, the first two canonical axes generated via the RDA are also statistically significant (p = 0.001) (Table S7).The adjusted RDA model explained 27.63% of the variation in these four explanatory variables distributed over the first two axes of the RDA.These play an important role in the dispersion of the species sampled at the four sites throughout the year.The RDA triplot (scaling 1,2) opposes along the first axis, from left to right, the coldest months with high rainfall (Group 1) to the warmest and driest months (Group 2).The second axis separates sites according to elevation (Figure 10).
An analysis using the non-metric multidimensional scaling (NMDS), with the Bray-Curtis dissimilarity matrix of the dung beetle species matrix (51 species/48 samples [= 4 sites × 12 months]), was used to graphically represent seasonal changes in the species communities.The correlation between the observed distance and the ordination distance (R 2 > 0.91) and the low-stress value (0.103) suggests that the NMDS fit is satisfactory.Figure 11 shows three seasonal clusters: one on the left, comprising the hottest months with no or reduced rainfall (from May to September); in contrast, another cluster comprising the coldest months with high rainfall (from November to February) was set up on the right.Between these two clusters, a third groupx shows the transition between them; on the one hand, this group shows the transition from the cold, rainy season to the dry, hot season (March and April); on the other, the opposite transition (October).0.001) (Table S7).The adjusted RDA model explained 27.63% of the variation in these four explanatory variables distributed over the first two axes of the RDA.These play an important role in the dispersion of the species sampled at the four sites throughout the year.The RDA triplot (scaling 1,2) opposes along the first axis, from left to right, the coldest months with high rainfall (Group 1) to the warmest and driest months (Group 2).The second axis separates sites according to elevation (Figure 10).An analysis using the non-metric multidimensional scaling (NMDS), with the Bray-Curtis dissimilarity matrix of the dung beetle species matrix (51 species/48 samples [= 4 sites x 12 months]), was used to graphically represent seasonal changes in the species communities.The correlation between the observed distance and the ordination distance (R 2 > 0.91) and the low-stress value (0.103) suggests that the NMDS fit is satisfactory.Figure 11 shows three seasonal clusters: one on the left, comprising the hottest months with no or reduced rainfall (from May to September); in contrast, another cluster comprising the coldest months with high rainfall (from November to February) was set up on the right.Between these two clusters, a third groupx shows the transition between them; on the one hand, this group shows the transition from the cold, rainy season to the dry, hot season (March and April); on the other, the opposite transition (October).To exhibit groups or factor levels on ordination diagrams, we used the ordihu functionality from the vegan package.The groups were derived from the Bray-Curt metric.Species can also be categorized into three groups: (1) species that prefer cold an wet months; (2) species that prefer hot and dry months, each group replacing the othe during the transition months (months 3, 4, and 10); and (3) species that were indifferen to the seasonal changes.The ordination of the first axis of the correspondence analys (CA) (Figure 12) confirms the results of this analysis.To exhibit groups or factor levels on ordination diagrams, we used the ordihull functionality from the vegan package.The groups were derived from the Bray-Curtis metric.Species can also be categorized into three groups: (1) species that prefer cold and wet months; (2) species that prefer hot and dry months, each group replacing the other during the transition months (months 3, 4, and 10); and (3) species that were indifferent to the seasonal changes.The ordination of the first axis of the correspondence analysis (CA) (Figure 12) confirms the results of this analysis.

Discussion
The organization of dung beetle communities depends on the combination of species traits (individual size, life cycles, and food preferences) and environmental variables (soil type, habitat characteristics, and elevation) [32,[45][46][47][48][49].Locally, the abundance of a species depends on the quantity and quality of the trophic resources available [28], the quantity of dung required for reproduction being directly linked to the size of the individuals [45,50].This relationship makes it possible to identify the most important species involved in the use of dung at any given time.Smaller species compensate for their size by many individuals to achieve the same level of functional efficiency as larger species, hence the combination of the criteria size (or weight) and number of individuals to identify the core and satellite groups.The organization of dung beetle communities in the Middle Atlas, characterized by very few dominant species at the same time, is like those observed in the mountains of alpine areas [23,30,[46][47][48], as well as in Mediterranean environments [29,46,47,[51][52][53][54][55], and in North America [32,49,56].
However, despite the elevational gradient of the sites under study (elevational variation of 1000 m), few differences were observed in the organization of their dung beetle communities year-round.In general, along an elevational gradient, their communities tend to become less diverse due to less favorable environmental conditions [48].In the present study, the gradient was too weak to observe a comparable situation.Non-nester species were active during the coldest and wettest periods of the year (late autumn and winter) and were later replaced by paracoprids and/or telecoprids during the rest of the year.This sharing of the activity periods by the guilds can be related to climatic

Discussion
The organization of dung beetle communities depends on the combination of species traits (individual size, life cycles, and food preferences) and environmental variables (soil type, habitat characteristics, and elevation) [32,[45][46][47][48][49].Locally, the abundance of a species depends on the quantity and quality of the trophic resources available [28], the quantity of dung required for reproduction being directly linked to the size of the individuals [45,50].This relationship makes it possible to identify the most important species involved in the use of dung at any given time.Smaller species compensate for their size by many individuals to achieve the same level of functional efficiency as larger species, hence the combination of the criteria size (or weight) and number of individuals to identify the core and satellite groups.The organization of dung beetle communities in the Middle Atlas, characterized by very few dominant species at the same time, is like those observed in the mountains of alpine areas [23,30,[46][47][48], as well as in Mediterranean environments [29,46,47,[51][52][53][54][55], and in North America [32,49,56].
However, despite the elevational gradient of the sites under study (elevational variation of 1000 m), few differences were observed in the organization of their dung beetle communities year-round.In general, along an elevational gradient, their communities tend to become less diverse due to less favorable environmental conditions [48].In the present study, the gradient was too weak to observe a comparable situation.Non-nester species were active during the coldest and wettest periods of the year (late autumn and winter) and were later replaced by paracoprids and/or telecoprids during the rest of the year.This sharing of the activity periods by the guilds can be related to climatic conditions (temperature and precipitations), which influence the moisture content of the dung and its quality.Nonnester species quickly colonize dung once it has been deposited.Their eggs are laid directly in the fecal matter, and the larvae eat and develop within the dung mass.This reproductive strategy is observed in most Aphodiinae, where females compensate for the non-protection of their eggs with a high fecundity [57].Larval development often takes several weeks, and it is crucial that the dung mass retains sufficient moisture content throughout this period.In this sense, autumn and winter were particularly favorable periods for non-nester species.They also avoid competition with the paracoprids and telecoprids, which require higher temperatures to reproduce.Paracoprids and telecoprids use a substantial proportion of dung, which is unfavorable for the non-nester species whose larvae develop unprotected inside the dung pats.On the other hand, paracoprids reduce the direct constraints due to dry conditions by burrowing under the dung and accumulating the material in the soil profile, where it retains its initial moisture content.The accumulated reserves enable their larvae to develop in a protected environment without high competition.This parental effort is compensated by reduced fecundity in the females, as most of their offspring reach maturity.There are different ways of storing the food mass for larvae.Dung may be packed at the blind end of a gallery like a sausage (most Geotrupidae; Bubas and Cheironitis species) or like pellets (e.g., Onthophagus species).In the Coprini tribe, the brood mass is remodeled into several brood balls containing an egg (Copris species).The telecoprid species quickly compact fresh dung, which is shaped into a ball and rolled several meters away from the initial mass.It is then buried in an underground chamber where conditions ensure that the brood ball will retain sufficient humidity to support the development of the larvae [58].
Competition between paracoprids and telecoprids was restricted mainly to sharing the trophic resources available in the pastures, with no real competition for egg-laying sites: paracoprids exploited the volume of space below the droppings to dig their nests, while telecoprids dispersed in search of egg-laying sites.The rapid succession of species making up the core and satellite groups probably restricted interspecific competition.Only around 20% of species can be considered to have played a decisive role in the functioning of each site studied (Figures 2-5).Despite structural similarities between the core and satellite groups, their composition was often different, even though some species were generalists and participated simultaneously in the formation of the core and satellite groups at several sites, such as S. sacer (two sites), A. fimetarius (three sites), and E. crocatus (four sites).Conversely, some species were dominant and characteristic of a single site, such as Onitis ion and O. melitaeus at Fez-Sais.Interestingly, even sites 1 km apart showed differences in their core and satellite groups, as seen in Ifrane I and Ifrane II, indicating that local conditions can play an important role in structuring communities.
The strategy of certain species to continue their activity over a long period could possibly impede other species, which have insufficient trophic resources to develop large populations.This is potentially the case for G. sturmi at Fez-Sais, which remained active for four consecutive months (May to August) during the warm period, while A. castaneus (September-October) and A. cribricollis (November-February) relayed each other to use most of the trophic resources in autumn and winter, these three species offering limited opportunities to the other species in the community (Figure 2).Similarly, E. crocatus was extremely dominant at Immouzzer, where there were no other species present to compete for four months (March to June) (Figure 3).At Ifrane I, two core species were dominant: S. schaefferi (small roller), which remained active for half the year, and N. affinis (dweller), taking over between November and January (Figure 4).At Ifrane II, N. affinis and S. laticollis (roller) successively dominated the dung beetle community and were the more distinctive species in the functioning of the site (Figure 5).
As expected, differences in the life traits of species can contribute to their spatiotemporal segregation, allowing each of them to exploit the available trophic resources in turn.Although the abundance and/or biomass of accessory species (more than fifty percent of total species) were very low compared with the other active species at a given time to have a substantial ecosystem impact, their presence remained important, as some of them could constitute a functional reserve in the future in the event of disturbance to the system.The populations of some dominant species could collapse under unfavorable conditions, and these accessory species, which in the meantime have found only narrow trophic niches in which to maintain themselves, could eventually take over and provide the same functions as those provided by the current dominant species [28].
Biodiversity fluctuates over time and space, resulting from intricate interactions among abiotic factors [59].For insects, seasonality primarily hinges on three factors, namely resource availability, temperature, and precipitation [60].This is true even for dung beetle communities [61,62].The beetle populations changed over time due to dominant factors such as temperature, precipitation, and elevation, as illustrated by the redundancy analysis.This analysis divided the surveys into two main groups, winter and summer, respectively, with an additional group created to highlight the transition between the other two.The transition group encompassed months 3 and 4 (March-April), delineating the shift from the cold, wet season to the dry, hot period, whilst month 10 (October) showed the opposite change.
This change is illustrated by the composition of the two main taxonomic groups.The first one consisted of species present during the winter months, from November to February, when average monthly temperatures were below 14.5 • C at all sites.The second group comprised thermophilic species that tended to be present during summer, from May to September, when average monthly temperatures at all sites were above 20 • C. The months of March, April, and November, which were positioned between these two groups, were transitional periods.A third group consisted of indifferent and eurythermic species that may occur at any time of the year, with varying frequency.This group was quite diverse, with some species displaying a distinct inclination toward the summer season while others were more typical of the transitional months.
Our findings contribute to the understanding of the regional dynamics of dung beetle communities in the Moroccan Middle Atlas and enhance our insight into their diversity within an environment facing multiple climatic constraints.

Supplementary Materials:
The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/d15111138/s1.Table S1.Monthly participation from September 2018 (month M9) to August 2019 (month M8) of core and satellite species in the Fez-Sais site.For each month, the left-hand column corresponds to the cumulative numbers (in %) of individuals, and the right-hand column corresponds to their biomass (%).Table S2.Monthly participation from September 2018 (month M9) to August 2019 (month M8) of core and satellite species in Immouzzer.For each month, the left-hand column corresponds to the cumulative numbers (in %) of individuals, and the right-hand column corresponds to their biomass (%).Table S3.Monthly participation from September 2018 (month M9) to August 2019 (month M8) of core and satellite species in the Ifrane I site.For each month, the left-hand column corresponds to the cumulative numbers (in %) of individuals, and the right-hand column corresponds to their biomass (%).Table S4.Monthly participation from September 2018 (month M9) to August 2019 (month M8) of core and satellite species in the Ifrane II site.For each month, the left-hand column corresponds to the cumulative numbers (in %) of individuals, and the right-hand column corresponds to their biomass (%).Table S5.Monthly Headcount of Core and Satellite Species Recorded from September 2018 to August 2019.Table S6.Data Matrix of Mesological Observations for the Period September 2018-August 2019.Table S7.Numerical Results for Figures 8 and 10.
program administered through the Department of Agriculture, Fisheries and Forestry.Contract Service Agreement between CSIRO and Faculty of Sciences Dhar el Mahraz n • 2019030623.

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki and approved by the General Director of Office National de Sécurité Sanitaire des Produits Alimentaires (ONSSA), Rabat (Morocco) for the delivery of successive collect permits.

Figure 1 .
Figure 1.Location of the study sites in Morocco.(A) Fez-Sais; (B) Immouzzer; (C) Ifrane I; and (D) Ifrane II.Red stars correspond to the position of the study sites.The Immouzzer station (33°47′ N-4°59′ W; elevation 898 m) is located on the Middle Atlas plateau, 3 km north of the town of Immouzzer Kandar.The substratum consisted of superficial red clay soil with outcrops of limestone blocks.The scrubland alternated between very open areas with a predominance of herbaceous vegetation and a few patches of Chamaerops humilis L. (locally called doum) and more closed areas (holm oak (Quercus ilex L.) coppice).Two sites, one kilometer apart, respectively, Ifrane I (33°32′ N-5°09′ W; elevation 1631 m) and Ifrane II (33°33′ N-5°10′ W; elevation 1613 m) were surveyed on the Middle Atlas plateau, four kilometers from the town of Ifrane.The sites were in the humid bioclimatic zone with cold winters.The landscape was very open and rocky, with a superficial red soil of silty clay resting on a bedrock of Lias limestone and dolomite.The vegetation was a mixture of herbaceous steppe formations with numerous clumps of Chamaerops humilis, indicating intensive grazing by sheep flocks.These open areas were bordered by a holm oak forest.

Figure 1 .
Figure 1.Location of the study sites in Morocco.(A) Fez-Sais; (B) Immouzzer; (C) Ifrane I; and (D) Ifrane II.Red stars correspond to the position of the study sites.The Immouzzer station (33 • 47 N-4 • 59 W; elevation 898 m) is located on the Middle Atlas plateau, 3 km north of the town of Immouzzer Kandar.The substratum consisted of superficial red clay soil with outcrops of limestone blocks.The scrubland alternated between very open areas with a predominance of herbaceous vegetation and a few patches of Chamaerops humilis L. (locally called doum) and more closed areas (holm oak (Quercus ilex L.) coppice).Two sites, one kilometer apart, respectively, Ifrane I (33 • 32 N-5 • 09 W; elevation 1631 m) and Ifrane II (33 • 33 N-5 • 10 W; elevation 1613 m) were surveyed on the Middle Atlas plateau, four kilometers from the town of Ifrane.The sites were in the humid bioclimatic zone with cold winters.The landscape was very open and rocky, with a superficial red soil of silty clay resting on a bedrock of Lias limestone and dolomite.The vegetation was a mixture of herbaceous steppe formations with numerous clumps of Chamaerops humilis, indicating intensive grazing by sheep flocks.These open areas were bordered by a holm oak forest.

Figure 2 .
Figure 2. Monthly organization of the dung beetle community in Fez-Sais from September 20 (month M9) to August 2019 (month M8).Core species names are mentioned in bold characte satellite species are underlined.Squares and diamonds correspond to the relative abundance an biomass, respectively, of species that were part of the monthly groups of core species.The size the squares and diamonds was adjusted to the relative part of each species.

Figure 2 .
Figure 2. Monthly organization of the dung beetle community in Fez-Sais from September 2018 (month M9) to August 2019 (month M8).Core species names are mentioned in bold characters; satellite species are underlined.Squares and diamonds correspond to the relative abundance and biomass, respectively, of species that were part of the monthly groups of core species.The size of the squares and diamonds was adjusted to the relative part of each species.In spring and summer (March to August), the non-nester species disappeared and were replaced by paracoprids and telecoprids.At Fez-Sais, three species were dominant in May in terms of numbers: Euonthophagus crocatus (paracoprid), Gymnopleurus flagellatus (telecoprid), and Gymnopleurus sturmi (telecoprid), which together accounted for 96.4% of the total number of individuals in the community, while G. flagellatus and G. sturmi together accounted for 93.4% in biomass.At Immouzzer, the small paracoprids were mainly represented by E. crocatus, which was active from March to June, with 62% of the total biomass of active species in April.At Ifrane I, in August, two telecoprids (G.sturmi and Sisyphus schaefferi) and two paracoprids (Bubas bison and Onthophagus maki) together accounted for 100% of the active dung beetles.At Ifrane II, from March to August, telecoprids (S. sacer, S. laticollis, G. flagellatus, G. sturmi, and S. schaefferi) and four small paracoprids (O.maki, Onthophagus vacca L., O. m. andalusicus, and E. crocatus) were the main species that composed the dung beetle community.

Figure 3 .
Figure 3. Monthly organization of the dung beetle community in Immouzzer from September 2018 (month M9) to August 2019 (month M8).Core species names are mentioned in bold characters; satellite species are underlined.Squares and diamonds correspond to the relative abundance and biomass, respectively, of species that were part of the monthly groups of core species.The size of the squares and diamonds was adjusted to the relative part of each species.

Figure 4 .Figure 3 . 20 Figure 3 .
Figure 4. Monthly organization of the dung beetle community in Ifrane I from September 2018 (month M9) to August 2019 (month M8).Core species names are mentioned in bold characters;

Figure 4 .
Figure 4. Monthly organization of the dung beetle community in Ifrane I from September 2018 (month M9) to August 2019 (month M8).Core species names are mentioned in bold characters;

Figure 4 .
Figure 4. Monthly organization of the dung beetle community in Ifrane I from September 2018 (month M9) to August 2019 (month M8).Core species names are mentioned in bold characters; satellite species are underlined.Squares and diamonds correspond to the relative abundance and biomass, respectively, of species that were part of the monthly groups of core species.The size of the squares and diamonds was adjusted to the relative part of each species.

Figure 5 .
Figure 5. Monthly organization of the dung beetle community in Ifrane II from Septembe (month M9) to August 2019 (month M8).Core species names are mentioned in bold chara satellite species are underlined.Squares and diamonds correspond to the relative abundanc biomass, respectively, of species that were part of the monthly groups of core species.The s the squares and diamonds was adjusted to the relative part of each species.

Figure 5 .
Figure 5. Monthly organization of the dung beetle community in Ifrane II from September 2018 (month M9) to August 2019 (month M8).Core species names are mentioned in bold characters;satellite species are underlined.Squares and diamonds correspond to the relative abundance and biomass, respectively, of species that were part of the monthly groups of core species.The size of the squares and diamonds was adjusted to the relative part of each species.

Figure 7 .
Figure 7. Scatter plots and pairwise correlation matrix of alpha diversity indices (Hill's diversity numbers N and Hill's ratio E) obtained from beetle's abundance data.Significant correlations are in bold, with stars at the bottom.p-values: ** ≤ 0.01, *** ≤ 0.001 Red lines represent the lowess smoothing.

Figure 7 . 20 Figure 8 .
Figure 7. Scatter plots and pairwise correlation matrix of alpha diversity indices (Hill's diversity numbers N and Hill's ratio E) obtained from beetle's abundance data.Significant correlations are in bold, with stars at the bottom.p-values: ** ≤ 0.01, *** ≤ 0.001 Red lines represent the lowess smoothing.5, x FOR PEER REVIEW 12 of 20

Figure 8 .
Figure 8. Boxplots of six diversity indices grouped according to the sites: Fez-Sais (A); Immouzzer (B); Ifrane I (C); and Ifrane II (D).The significant overall difference between sites after analysis of variance.The letters above indicate whether there is a significant difference between sites in pairwise post hoc comparisons (the same letter indicates no difference); ns: not significant.p-values: * ≤ 0.05; ** ≤ 0.01.

Figure 8 .
Figure 8. Boxplots of six diversity indices grouped according to the sites: Fez-Sais (A); Immouzzer (B); Ifrane I (C); and Ifrane II (D).The significant overall difference between sites after analysis of variance.The letters above indicate whether there is a significant difference between sites in pairwise post hoc comparisons (the same letter indicates no difference); ns: not significant.p-values: * ≤ 0.05; ** ≤ 0.01.

Figure 8 .
Figure 8. Boxplots of six diversity indices grouped according to the sites: Fez-Sais (A); Immouzzer (B); Ifrane I (C); and Ifrane II (D).The significant overall difference between sites after analysis of variance.The letters above indicate whether there is a significant difference between sites in pairwise post hoc comparisons (the same letter indicates no difference); ns: not significant.p-values: * ≤ 0.05; ** ≤ 0.01.

Figure 11 .
Figure 11.Non-metric multidimensional scaling (NMDS) of samples categorized into three groups defined via the Bray-Curtis hierarchical clustering and Ward's D2 method.The letters identify the four sites (Fez-Sais (A), Immouzzer (B), Ifrane I (C), and Ifrane II (D)) and the numbers of the sampling months (e.g., A9: site A, sampling month 9).

Figure 12 .
Figure 12.Ordination of samples along the F1 axis of a correspondence analysis (CA).The matrix is composed of 48 samples/51 species-samples (months x sites).The circles reflect the species' abundance transformed into log(x + 1), and their sizes are proportional to the abundance.

Figure 12 .
Figure 12.Ordination of samples along the F1 axis of a correspondence analysis (CA).The matrix is composed of 48 samples/51 species-samples (months x sites).The circles reflect the species' abundance transformed into log(x + 1), and their sizes are proportional to the abundance.