Ancient human colonization explains dung beetle species richness in the Mediterranean and Macaronesian islands

Different hypothesis have been proposed to explain differences in species richness among islands. However, few studies have attempted to compare the explanatory power of multiple hypotheses using a large data set. Here, we analyse how different types of predictors (energetic/climatic, environmental heterogeneity, island biogeography and anthropogenic) affect variation in dung beetle species richness on Mediterranean and Macaronesian islands.


| INTRODUC TI ON
Islands could be considered a model system for ecological studies (Vitousek, 2002).Due to their clear boundaries, great number and variability in shapes, size, geology and environments, islands are natural laboratories for biogeographers and ecologists who can study biodiversity patterns in apparently simpler conditions avoiding the difficulties imposed by the continuity in environmental and spatial conditions that often appear in continental areas (Losos & Ricklefs, 2009;Whittaker & Fernández-Palacios, 2007).Spatial variation in species richness has been one of the most studied biodiversity variables, and a plethora of literature has been published discussing what factors determine such patterns both within continental areas (Keil & Chase, 2019) and between islands (Whittaker & Fernández-Palacios, 2007).Hypotheses that attempt to explain the causes of variation in species richness can be summarized into four main ones depending on the explanatory variables used: energetic/climatic, environmental heterogeneity, island biogeography and historical/anthropogenic.The energetic/climatic variables are fundamental to account for variation in species richness (Wright, 1983).This hypothesis considers that areas or islands with more available energy and more primary productivity can support more individuals, allowing species to reach higher population sizes, thereby reducing extinction rates, increasing diversification rates and promoting species richness (Evans et al., 2005;Gaston, 2000;Hawkins et al., 2003;Kalmar & Currie, 2006;Kreft et al., 2008).Other hypotheses that attempt to explain species richness variations are based on niche theory, which emphasizes the role of environmental heterogeneity and habitat diversity as key factors structuring ecological communities (Hortal et al., 2009;Kadmon & Allouche, 2007) because a more heterogeneous environment could support more species through niche partitioning (Klopfer & MacArthur, 1960).Island biogeography variables include some of those traditionally used in the Equilibrium Theory of Island Biogeography (ETIB) (MacArthur & Wilson, 1963, 1967), plus others that try to measure the colonization capacity of the islands by the species.This theory states that the number of species on islands is determined by an equilibrium between immigration and extinction rates, both of which depend on the area of the islands and their geographical isolation.Despite its great explanatory capacity using simple derived variables, the ETIB was also criticized from an empirical and theoretical point of view because of its oversimplification and because the omission of factors such as the dynamic nature of islands over time or habitat diversity and status (Bush & Whittaker, 1991;Gilbert, 1980).This led the way to consider other factors to explain species richness variation, such as time and islands ontogeny, thus developing the so-called general dynamic model of oceanic island biogeography (Fine, 2015;Fine & Ree, 2006;Whittaker et al., 2008).Finally, it is not only contemporary factors that can affect the variation in species richness.The contingent history of different regions and taxa may also be influencing current species richness values.
Among the historical contingent factors, it is particularly important to highlight the role of human influences on ecological assemblages, both in a positive or negative sense (Castro et al., 2010;Chown et al., 1998;Fattorini et al., 2016;Rizali et al., 2010;Walentowitz et al., 2022).Human activities can influence species richness patterns by increasing the probability of propagule transfer, the likelihood of introduction of invasive species, or by modifying habitat and resource availability (Nakamura et al., 2015).Some authors suggest that anthropogenic factors may come to play a pre-eminent role relative to other biogeographical processes (Helmus et al., 2014;Nakamura et al., 2015;Walentowitz et al., 2022).
Empirical evidence has shown that climatic/energetic variables are important in determining dung beetle species richness.These variables directly affect the physiology of dung beetles (water balance, thermoregulation and thermal tolerance), also influencing trophic resource use or nesting in the soil (wet or dry) (Bogoni et al., 2019;Chefaoui et al., 2005;Davis et al., 2008;Gebert et al., 2020;Krell et al., 2003;Lobo et al., 2002;Roslin et al., 2009;Verdú et al., 2007;Vessby, 2001).In the case of environmental heterogeneity, several studies have confirmed the general positive effect of increasing available niche space on dung beetle diversity (Negro et al., 2011;Numa et al., 2009;Tonelli et al., 2019).Island biogeography variables have also been shown to be relevant in shaping dung beetle assemblages (Fattorini, 2010).Regarding contingent and historical factors, here we are specifically interested in determining the explanatory capacity of some variables representing human actions on the species richness of island dung beetles, as previously suggested (Tonelli et al., 2016).Anthropogenic activities are recognized as one of the greatest threats to species diversity, and the intensity of this threat is often associated with human history and population growth (Davies et al., 2006;McKee et al., 2003;McKinney, 2001;McNeely et al., 1995).However, human presence can exert a beneficial effect on some species (Gaston, 2005;Luck, 2007), especially when anthropogenic influence is moderate.
In the case of dung beetles, a decline in their populations has recently been observed (Lobo et al., 2006;Lumaret et al., 2022;Tonelli et al., 2018) as a consequence of the abandonment of extensive livestock farming and the use of anti-parasitic treatments in livestock, which are excreted unaltered and poison the beetles (Verdú et al., 2018).However, the close link between these species and the livestock practices in the Mediterranean since the Neolithic (Ethier et al., 2017) must have promoted the expansion of the populations of some species (Kang et al., 2018), but not the disappearance of native dung beetle fauna.Pleistocene fossil evidence indicates that the dung beetle species present before the development of animal anthropic impact, Aphodiidae, Geotrupinae, Scarabaeinae, species richness hypotheses, survey effort husbandry are the same as those found today (Coope, 2001(Coope, , 2010;;Sandom et al., 2014).Dung beetles have been and are a key taxon in Mediterranean grassland ecosystems (Numa et al., 2020), and their feeding and nesting activities are linked to the use of mammal dung (Halffter & Matthews, 1966).Although some species may have a degree of preference for specific types of dung (Giménez Gómez et al., 2023;Tonelli et al., 2021;Verdú & Galante, 2004), they have a highly generalist and opportunistic feeding behaviour, using both domestic and wild mammal excrements interchangeably (Frank et al., 2018;Martín-Piera & Lobo, 1996).
Different hypotheses able to explain species richness difference on islands were tested using different focal taxa.However, to our knowledge, few studies have evaluated their comparative importance (but see Aranda et al., 2014;Cardoso et al., 2010;Fattorini, 2011;Kalmar & Currie, 2006).Comparative evaluation of multiple competing hypotheses is crucial to promote objectivity, minimize bias, stimulate discovery, strengthen theories and advance our understanding of the natural world (Betini et al., 2017).Moreover, the majority of this research was based on studying a single archipelago, generating idiosyncratic results that are difficult to generalize to other contexts (Chiarucci et al., 2011).Additionally and importantly, almost all the biogeographical studies based on data derived from the literature assume the exhaustiveness of the data and their completeness, considering that the used inventories are reliable and that the available data do not need to be controlled for survey effort disparities (e.g.Borges & Hortal, 2009;Chown et al., 1998).
However, the addition of the survey effort could be a key factor to manage the unevenness in data quality and geographical coverage of the different areas (Aranda & Lobo, 2011;Gray & Cavers, 2014;Hortal et al., 2007;Santos et al., 2010).In this study, we analyse for the first time the effects of different types of variables representative of the aforementioned hypotheses (energetic/climatic, environmental heterogeneity, island biogeography and anthropogenic) on the variation in dung beetle species richness across Mediterranean and Macaronesian islands.Our main aim is to compare the explanatory capacity of these variables in these two bioregions in order to assess the relevance of the different hypothesis and, specifically, to examine whether anthropogenic activity can be considered a key factor in accounting for variation in dung beetle richness.

| Geographical extent and taxonomic groups
The study focused on the islands of the Mediterranean and Macaronesian regions, totalling 147 islands (Figure 1; Appendix S1 and S2).
Islands were included in this study when they had faunistic studies or published records of dung beetle species.The great variability in the environmental features of all these islands maximize the capacity of finding geographical species richness patterns able to be explained by the used explanatory variables and factors.
The dung beetle species considered were those belonging to the superfamily Scarabaeoidea, which mainly feed on mammal faeces at any stage of their development (Krell & Moon, 2019).In other words, dung beetles strictly belong to the subfamilies Geotrupinae, Scarabaeinae and Aphodiinae (Krell & Moon, 2019).However, in our analysis, we used data of dung beetles sensu lato because we added data of all the family Aphodiidae, which also includes the subfamilies Psammodiinae, Aegialiinae and Eremazinae, which have F I G U R E 1 Geographical location of the 147 Mediterranean and Macaronesian studied islands.The photos of beetles represent (from left to right) the typical aspect of two species of Scarabaeinae, one of Aphodiidae and another of Geotrupinae.a saprophagous feeding behaviour (Dellacasa & Dellacasa, 2006).Thus, we studied here the geographical variation in the species richness of dung beetle species belonging to Geotrupinae, Scarabaeinae and Aphodiidae.
Geotrupinae is phylogenetically the oldest taxonomic group (Ahrens et al., 2014;Cunha et al., 2011;Gunter et al., 2016) and includes around 150 species of medium/big size (≈20 mm) with a relatively homogeneous morphology.During nesting, they essentially show a paracoprid behaviour (Tonelli, 2021; but see Zunino & Palestrini, 1986), excavating tunnels and constructing their nest in the soil beneath the food source.Geotrupinae could be defined as trophic generalist due to their ability to exploit a variety of dung types and some plant detritus, fungus and even acorns (Sánchez-Piñero et al., 2019).Overall, we found 22 Geotrupinae species in the considered islands.Aphodiidae include more than 3000 species distributed worldwide, which are characterized by their small body size (<20 mm) (Dellacasa & Dellacasa, 2006).We can consider Aphodiidae, with a few exceptions, as non-nesters because they directly use the trophic source without any manipulation, relocation or nest building, and the larvae live freely inside dung pats (Tonelli, 2021).Aphodiidae can consume a wide range of trophic resources, and it is possible to find both coprophagous and saprophagous species.Overall, we found 135 Aphodiidae species in the studied islands.Finally, the Scarabaeinae subfamily includes more than 6000 species, showing a great diversity of morphological, ecological and behavioural adaptations to the consumption of mammal faeces (Halffter & Matthews, 1966;Tonelli, 2021).Overall, in our study, we included 63 Scarabaeinae species in the considered islands.
All available published manuscripts concerning dung beetles in the studied areas were compiled, such as atlases, books, peerreviewed papers and grey literature.Most of these references come from the authors' exhaustive library.A total of 362 bibliographic reference sources published from 1838 to 2019 were used to extract dung beetle occurrence data (see Appendix S1 for more details).All these data were checked, and the nomenclature was standardized following an up-to-date taxonomical criterion (Löbl & Löbl, 2016) avoiding subspecies, recent species splitting and doubtful citations.
All dung beetles species finally listed (n = 220) are considered native to the biogeographic region in which the islands studied are located.

| Explanatory predictors
We tested four types of causal hypotheses able to explain the variation in species richness: energetic/climatic, environmental heterogeneity, island biogeography and anthropogenic.
Energetic/climatic variables (EC) come from WorldClim 2.1.(Fick & Hijmans, 2017; https://www.worldclim.org) and ENVIREM (Title & Bemmels, 2018; https://envir em.github.io/)totalling the following 11 variables: annual mean temperature, annual precipitation, annual potential evapotranspiration, Thornthwaite aridity index, growing degree days with mean temperature >5°C, maximum temperature of the warmest month, minimum temperature of the coldest month, precipitation of the wettest month, precipitation of the driest month, total solar radiation and total wind speed.All these variables were rescaled to a 2.5-minute resolution (a cell of approximately 17 km 2 at the 39° latitude) and averaged values calculated for each island using the capabilities to handle digital mapping of the software ModestR version 6.3 (García-Roselló et al., 2013).The obtained EC values for each island were submitted to a principal component analysis (PCA) to reduce the dimensionality of all these variables by computing orthogonal uncorrelated variables able to represent these frequently collinear variables, while preserving as much of the variability in the original data as possible.Three components with eigenvalues higher than one, accounting for 87.9% of total variability, were obtained (58.8%, 18.1% and 11.0%, respectively) and their factor scores were selected as the variables representing the EC influence.Factors loadings indicated that the first component (temperature component) is positively correlated with temperature variables (e.g.0.907 for annual mean temperature).
The second component (evapotranspiration component) is positively correlated with annual potential evapotranspiration (0.760), while the third component (aridity component) is negatively correlated with precipitation variables (e.g.−0.917 for annual precipitation).
Two variables were used to account for the environmental heterogeneity (EH) hypothesis: island area and island maximum elevation.Island area has been also included within island biogeography variables since it can affect species richness by their direct effect on populations (larger islands intercept more immigrants and support more individuals).However, island area can also be relevant in explaining species richness variations because of their positive relationship with the available environmental niche space and their capacity to facilitate coexistence by diminishing interspecific competition (Hortal et al., 2009;Stein et al., 2014).The term environmental heterogeneity is used to refer to the environmental variability within an area, including topography, habitat type and climate (Hortal et al., 2009).Island maximum elevation is also one of the most used proxy variables for environmental heterogeneity (e.g.see Aranda et al., 2014;Barajas-Barbosa et al., 2020;Weigelt et al., 2013).Elevational gradients can cause marked differences in temperature, orographic precipitation regimes and rain shadow effects over short geographical distances (Barajas-Barbosa et al., 2020).
Five variables were used to test island biogeography (IB) hypothesis: island connection to continent during last glacial maximum (LGM), island geological origin, island area, distance from continent or from the nearest big island and nearest continent.We consider that all these variables may affect island species richness representing the classical ETIB variables and others trying to measure the possibilities of colonization of these islands by species.
Distances to the nearest mainland were calculated using the Google Maps Distance Measurement Tool as the distance in kilometres from the coastline of each island to the coastline of the nearest continental mainland or nearest big island.The nearest continent (Europe or Africa) was determined using the Google Maps Distance Measurement Tool and was also included as a dummy variable because the environmental conditions of the source regions (Wiens & Donoghue, 2004) and the species pool (Ricklefs, 2004) inhabiting these source regions may have conditioned the current distribution of species.Island connection to mainland during last glacial maximum is widely recognized as an important variable in explaining species richness because increased connectivity with the mainland may result in higher establishment rate of continental species (Fernández-Palacios et al., 2016;Hammoud et al., 2021).The island geological origin (volcanic or continental) is a variable that takes into account the ontogenetic evolution of the island, as proposed in the general dynamic theory of oceanic islands (Whittaker et al., 2008).
Islands originated as unsubmerged parts of the continental shelf, and islands formed as the result of volcanic actions may have completely different evolutionary dynamics affecting species richness (Ali, 2017;Whittaker et al., 2017).Island area and the distance from the nearest mainland are included because they are classical key variables able to determine extinction and immigration rates, which regulate the species richness reached at a dynamic equilibrium (MacArthur & Wilson, 1967;Whittaker et al., 2017).Lastly, the dummy variable representing the nearest continent can be a decisive factor that aims to represent the possible role played by the different biogeographical history of African and European landmasses and its effects on species richness (Sanmartín, 2003).
To test the possible role of anthropogenic (A) influences on species richness, two variables have been used: human density and years since first human colonization.Human density (current inhabitants per km 2 ) was extracted from Wikipedia (https://es.wikipedia.org).
Human population density and species richness usually show a significant positive correlation up to a certain threshold beyond which species richness decreases (Paradis, 2018).Thus, this variable can have both a positive and a negative effect of species richness.Years since first human colonization was obtained from archaeological literature and the kind assistance of expert archaeologists (Dr.Helen Dawson and Prof. Stašo Forenbaher; cfr Appendix S1 for more details).When possible, the time of the first stable colonization by human population has been included.However, in the case of some small islands, the available archaeological data do not allow a clear determination of the probable date of the first settlement.In these cases, the listed year for first presence of humans is often based on a single archaeological characteristic (such as a decorated potsherd or stone axe) that can be unambiguously attributed to a specific time period.In total, we have been able to extract data of early human colonization for 132 islands.
In addition to the above-mentioned variables, we included as a covariate the number of published papers studying dung beetles in each of the islands.This variable may act as a surrogate of the survey/study effort carried out in each island (the more published studies, the more species are observed).However, the number of published studies can be also related with other characteristics of the islands such as its area or its own species richness (e.g. the greater the species richness, the greater the number of studies that are carried out).Thus, we compared the results provided by regressions using the selected predictors against those using as response variable the residuals of a regression between species richness and the number of published papers.That is, models considering that the possible difference in survey effort between the islands hardly influences species richness variations, and models in which these differences in survey effort are fundamental to understand species richness patterns.We assumed that a hypothesis maintaining its comparative relevance in the two approaches would be more likely to causally affect species richness variations.

| Statistical analyses
Separate generalized linear models (GLM; McCullagh & Nelder, 1989) were fitted to summarize the relationship between the species richness of the three main groups of dung beetles (Scarabaeinae, Aphodiidae and Geotrupinae) and the above-mentioned explanatory variables.The purpose was to examine how much of the variation in the species richness of islands is explained by the different sets of explanatory variables, which aim to represent alternative/complementary causal hypotheses.A Poisson error distribution for the number of dung beetle species was assumed, which is linked to the predictor variables via a logarithmic link function.All the explanatory variables are previously standardized to mean zero and one standard deviation to eliminate the effect of the differences in the measurement scale.Subsequently, the number of species in each island is firstly related with each predictor separately to examine the individual predictive capacity of each variable.As the relationship between species richness and environmental variables is often curvilinear, either a linear or a quadratic function of each environmental variable was selected.Quadratic functions are selected when both the linear and the quadratic terms are statistically significant at the 0.01 level.The goodness of fit of the models was measured by the statistic deviance, comparing the reduction in deviance of each predictor with a full model in which the number of parameters was equal to the total number of observations (Dobson, 1990).The Wald statistic was used for testing the statistical significance of the regression coefficients based on maximum likelihood estimates.
Next, the explanatory capacity of the four groups of variables representing causal hypotheses (EC, A, IB and EH) was calculated also using GLMs.As above, we are especially interested in the percentage of deviance accounted for each group of predictors.These regressions have been repeated using as response variable the residuals of a GLM model between species richness and the number of published papers in each island.In this case, we assumed an identity link function and a Normal error distribution for the dependent variable.A dummy variable representing the Macaronesian islands was also added to these analyses to examine whether the explanatory capacity of the different groups of variables increased substantially.
The statistical significance (p ≤ 0.01) of the interactions between each of the explanatory variables and this dummy variable was also estimated to find out whether the effect of each of these predictors depends on the Mediterranean or Macaronesian distribution of the islands.Finally, a forward selection procedure was used to iteratively derive a complete model in which all the predictors are included if their statistical significance is equal or lower than 0.01.The residuals of this complete or saturated model (the variability unexplained by the used predictors) were regressed on a third-order polynomial of latitude and longitude (trend surface analysis; see Legendre & Legendre, 2012) to assess whether one or more important spatially structured explanatory variables have been omitted (Legendre & Legendre, 2012) or, conversely, whether the variables that have been used correctly represent the pattern of spatial variation.Geographical coordinates were previously transformed in a Euclidean coordinate system (UTM) and then standardized to zero mean and one standard deviation.All these statistics were performed using StatSoft's STATISTICA v10.0 (StatSoft Inc., Tulsa, Oklahoma, USA).

| RE SULTS
Individual GLMs suggest that the years since first human colonization (explaining 47%, 38% and 41% of total deviance), the number of published papers (explaining 24%, 48% and 18% of total deviance) and island area (23%, 41% and 17%) were the predictors with the highest explanatory capacity for Scarabaeinae, Aphodiidae and Geotrupinae species richness respectively (Table 1; see Appendix S3).
The two first variables were always positively related to species richness, so that the number of species increases when the time elapsed since the first human colonization is greater and the number of studies carried out is higher (Figure 2).Island area was curvilinearly related to Scarabaeinae and Aphodiidae richness, so that the number of species increases with the size of the island up to a certain threshold (i.e.+−; see Table 1 and Appendix S4).The volcanic character of the islands and the distance to mainland had some relevance in the case of Scarabaeinae and Geotrupinae (less species in volcanic and distant islands in the two taxonomic groups), and maximum elevation appear as relevant in the case of Scarabaeinae and Aphodiidae species richness (more mountainous islands tend to have more species).Of the remaining predictors, only the first PCA component (temperature component) exceeded 10% of explanatory capacity in Scarabaeinae and Geotrupinae, so that increases in temperature would imply a decrease in species richness (Table 1; Appendix S3).
The geographical variation of species richness in the considered islands (Figure 3a) and the geographical variation of the residuals between species richness and the number of published papers ( Note: Provided results are each explanatory variable separately considering both their linear (L) and quadratic (Q) functions against species richness as well as the sign of the relationship (negative or positive).% ExpDev = percentage of explained deviance on the full model in which no predictors are included.
p < 0.0001).Anthropogenic and island biogeography hypotheses about the variation in species richness were the ones that would have the strongest support due to their high predictive capacity in the three dung beetle groups, with percentages equal to or >40% each one (Figure 4).When the surrogate of survey effort is firstly included as a covariate (the number of published papers), these two sets of predictors continue to be comparatively highly relevant, although their joint predictive capacity diminishes reaching values oscillating from 43% to 53%.Energetic/climatic variables were relatively important in the case of Scarabaeinae and Geotrupinae (≈20%), but their relevance in the three dung beetle groups is similar and hardly decreases when the effect of survey effort is removed (≈18%).Finally, the set of variables representing the environmental heterogeneity hypothesis would also be relatively influential, but their relevance decreases abruptly when survey effort is accounted for (Figure 4).F I G U R E 3 (a) Geographical variation in species richness of dung beetles (Aphodiidae, Scarabaeinae and Geotrupinae) for the islands considered.The size of the circles is related to the number of species in a range of richness values from minimum (0 species) to maximum (114 species).(b) Geographical variation of residuals in the relationship between species richness and the number of published papers (used as surrogate of the survey/ study effort carried out in each island).Smaller circles correspond to negative residuals (minimum value = −26) where the number of species observed is lower than the number of species that should be present according to the effort, and larger circles correspond to positive residuals (maximum value = 35) in which the number of species observed is higher than the number of species that should be present according to the effort.
Complete GLMs using all the predictors included seven variables in the case of Scarabaeinae and Aphodiidae, explaining 62.1% and 63.6% of total deviance respectively (not showed).In the case of Geotrupinae, these complete models accounted for 53.3% of total variability in species richness, being composed by only three explanatory variables.A forward stepwise procedure using as predictors the nine terms of a trend surface analysis do not allow us to select any statistically significant parameter able to explain the geographical variation in the residuals of these complete models (F 9,122 = 1.08, p = 0.38 for Scarabaeinae; F 9,122 = 1.02, p = 0.38 for Aphodiidae; and F 9,122 = 0.47, p = 0.89 for Geotrupinae) (see Appendix S6).

| Determinant factors
To the best of our knowledge, our study is the first attempt to infer the main factors responsible for the variation in dung beetle species richness in Mediterranean and Macaronesian islands.In total, we found 220 dung beetle species, 200 species in Mediterranean islands and 38 species in Macaronesian islands (18 species common to both bioregions).This figure represents almost one-third of the 644 dung beetle species estimated to inhabit the whole Mediterranean region (Numa et al., 2020).All species occurring only on the Macaronesian islands are Aphodiidae, only three species of Scarabaeinae and no species of Geotrupinae have been reported from these islands.The available biogeographic knowledge about Palaearctic dung beetles (Löbl & Löbl, 2016) indicates that only 12 species can be considered endemic to some of the larger islands such as Sicily, Sardinia, Corsica, Menorca and Cyprus, which represents <2% of the total number of Mediterranean dung beetle species.This low rate of endemicity suggests that the populations of the dung beetle species inhabiting Mediterranean and Macaronesian islands probably have a recent origin (Veron et al., 2019).
The explanatory variables appearing as relevant in our study also this group present on these islands (0.6 species per island).The result is interesting because in general species richness on islands seems to suffer an impoverishment and homogeneity due to human influence (Cardoso et al., 2010;Longman et al., 2018;Russell & Kueffer, 2019;Sayol et al., 2021).The negative effect of human impact can sometimes be even greater than that due to 'natural' impacts such as volcanic eruptions or climatic changes (Connor et al., 2012).To survive, early settlers modify island ecosystems for agricultural and pastoral purposes, introduce animals and plants and exploit native species, often resulting in the extinction of native species (de Nascimento et al., 2020;Steadman, 1995).However, negative human impacts on biodiversity are more evident on recently colonized islands than on previously colonized islands, probably because the effects of introduced species, land use practices and technology deployed by early settlers are less impactful than those of later settlers (Nogué et al., 2021).
Although dung beetle richness on Mediterranean and Macaronesian islands is strongly and positively correlated with the years since first human colonization, it is poorly correlated with current human density.Therefore, we suspect that the explanation of the between-island variation in dung beetle species richness must be related to the ancient dispersal of human populations and the se- ( Blondel, 2008).These human movements could have promoted the differential extinction of some taxa and improved dung beetle colonization by involuntary human transport, as well as facilitating populations settling as a consequence of the presence of trophic resources derived from livestock and human faeces (Tonelli et al., 2016(Tonelli et al., , 2017)).Indeed, it is well documented that human populations dispersed domestic livestock during their voyages (Zeder, 2017), and this is also true for islands (Dawson, 2014).A region in which the landscape and the composition and abundance of its mammalian fauna have been modified by human actions since Neolithic times must necessarily and profoundly affect the biogeography of the main group of species associated with the consumption of mammal faeces (Martín-Piera & Lobo, 1996;Verdú & Galante, 2002).It is well known that the quantity of trophic resources in the Mediterranean region is a much more relevant variable to explain dung beetle species richness variations than the variety of trophic resources (Lumaret et al., 1992).This implies that moderate growth in agricultural practices and increased livestock density would have had a beneficial effect on dung beetle populations and species richness.Furthermore, both fossil data (Coope, 2001(Coope, , 2010;;Sandom et al., 2014) and the general lack of high trophic specialization (Frank et al., 2018) suggest that human actions have not led to the replacement of a pre-existing native dung beetle fauna, but to the restructuring of dung beetle assemblages by promoting dispersal and favouring certain types of species and functional morphotypes.
Island biogeography variables also showed a relative high explanatory capacity, especially island area, island geological origin (volcanic or continental) and distance from mainland.Larger sized islands might have a higher availability and heterogeneity of resources, host larger populations of dung beetles and show a lower extinction risk than smaller islands, where stochastic events may play important roles (Qie et al., 2011).Furthermore, we cannot exclude the possibility that larger islands may have also facilitated ancient human settlement and allowed for a higher density of domestic and wild mammals (Connor et al., 2000), thus indirectly promoting dung beetle richness.The fact that volcanic islands tend to host less dung beetle species could be related to the disturbance of periodic eruptive events, the comparatively lower age (Whittaker et al., 2008) and the difficulty of being colonized and sustaining large and prolonged human activity.The negative role of the volcanic character of the islands on dung beetle species richness is particularly evident in the Macaronesian islands.
The other studied variables also show a relatively important predictive capacity, but their relevance oscillates depending of the taxonomic group (see Table 1 and Figure 4).Maximum elevation is relatively relevant in Scarabaeinae, but accounts for a significant proportion of the variability in species richness in the case of Aphodiidae, especially in Mediterranean islands.Although Aphodiidae is basically a Holarctic group well adapted to cold-temperate conditions, there are also many warm-adapted lineages with a probable Afrotropical origin (Cabrero-Sañudo & Lobo, 2009;Dellacasa et al., 2001).Due to this, heterogeneous mountain areas should be especially suitable to harbour a high diversity of Aphodiidae species with different environmental preferences, as well as acting in the Mediterranean region as refuges during Pleistocene glaciations (Cabrero-Sañudo & Lobo, 2006).Similar reasoning could be adduced to explain the greater influence of island area on Aphodiidae relative to Scarabaeinae and Geotrupinae, thus showing the positive influence of environmental heterogeneity and mountain areas to maintain a high number of species.In the case of Geotrupinae, environmental temperature appears negatively associated with the number of species of this eminently Holarctic group, which only in its more recent evolutionary history showed adaptations to arid and semi-arid environments (Cunha et al., 2011).

| The importance of survey effort
Species richness is the most used measure of biodiversity, representing the milestone in conservation biology (Fleishman et al., 2006), community ecology (Magurran & McGill, 2011), ecological research (Hooper et al., 2005) and a mechanistic understanding of species richness patterns remains the holy grail of modern biogeography and macroecology (Gotelli et al., 2009).Unfortunately, species richness is very sensitive to sampling intensity and species abundance distribution (Chao et al., 2014), so it could be strongly biased across a large geographical scale.Hence, in macroecological studies that use distributional data from different sources, with different aims and without a standardized sampling protocol, we cannot know whether inventories are complete and comparable (Hortal et al., 2007), unless covariates such as the number of database records have been used to try to control for this problem (Hortal et al., 2006;Lobo, 2008).Our complete model results show a relatively important explanatory capacity (≈60% of total variability) without leaving any apparent spatial pattern unexplained.This is partially due because we included the number of published papers as a surrogate for survey effort, and its strong correlation with dung beetle species richness demonstrates the importance of controlling this factor in large-scale biogeographical studies.Indeed, all the hypotheses tested lose explanatory power after controlling for survey effort, but not in an equal manner.Environmental heterogeneity variables loss about 80% of their explanatory power among the three dung beetle groups, and becomes the least important variable in explaining species richness when the effect of survey effort is ruled out.This result is probably due to the increasing level of survey carried out in those islands, thus showing a higher heterogeneity.We consider that using some surrogate of the survey effort carried out should always be an option in macroecological and biogeographical studies.

| CON CLUS IONS
The main conclusion of our study is that the most relevant factors grazing practices, as has been widely and recently recognized (Numa et al., 2020).

ACK N O WLE D G E M ENTS
We are very grateful to Prof. Helen Dawson and Prof. Stašo Forenbaher for their great help in the management and interpretation of archaeological data.We also want to thank Prof. Mario Zunino for his great help in finding and sharing bibliographic references.This study is based on data previously collected and published by several authors, so no collection permissions were required.We acknowledge support of the publication fee by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).

CO N FLI C T O F I NTE R E S T S TATE M E NT
None.
Figure 3b) are positively and significantly correlated (Pearson r = 0.657;TA B L E 1 Results of generalized linear models explaining the variation in island species richness for the three main dung beetle groups using as predictors the 12 considered variables. The explanatory capacity of the different hypotheses hardly increased (between 1% and 5%) when the Macaronesian or Mediterranean distribution of the islands is included in the models (see Appendix S5).The only exception is the environmental heterogeneity hypothesis in the case of Aphodiidae species, which increases from 5.6% to 17.9% due to the differential effect of maximum elevation.In Macaronesia, the number of Aphodiidae species hardly increases with maximum elevation (1.6 species per km), whereas in the Mediterranean islands, species richness increases considerably with maximum island elevation (13.5 species per km).In Aphodiidae, three other predictors showed opposite relevance between islands of the two bioregions: years since first human colonization, PCA evapotranspiration component and nearest continent (see Appendix S5).Thus, the increase in the number of Aphodiidae species with the age of human colonization is greater in Macaronesia (4.8 species per 1000 years) than in the Mediterranean islands (2.8 species per 1000 years), although human presence is much more recent.The increase in species richness with increasing evapotranspiration is almost four times higher in the Mediterranean islands, and the proximity to Europe has a slightly smaller negative effect on Aphodiidae in the Macaronesian islands.In the case of Scarabaeinae, the years since the first human colonization have little influence on the species richness in the Macaronesian islands (−0.6 species per 1000 years) while they positively influence species richness in the Mediterranean islands (1.6 species per 1000 years).Scarabaeinae species richness decreases about three times faster with increasing temperature in the Mediterranean islands, while increasing evapotranspiration decreases species richness in Macaronesia but increases it in the Mediterranean islands (see Appendix S5).

F
Between island variation in the species richness of Aphodiidae (red squares), Scarabaeinae (blue circles) and Geotrupinae (green triangles) with the years since first human colonization (in calibrated years before present).Broken lines represent the curvilinear relationship of the data (quadratic polynomial).
support the deep and frequent influence of continental faunas on the extant dung beetle diversity of Mediterranean and Macaronesian islands.The most important variables accounting for dung beetle species richness are those related to anthropogenic and island biogeography/colonization, particularly the years since first human colonization and to a lesser extent island area.The years since first human colonization seeks to reflect the age of the agricultural and livestock practices with which the activity of these beetles is associated (Jay-Robert et al., 2008; Tonelli et al., 2019), being the most powerful predictor of the explained deviance for the three groups of dung beetles.The comparative importance of this variable is maintained when the effect of survey effort differences is ruled out, always showing a positive effect on dung beetle richness.The positive effect of the age of human colonization is only not observed in the case of the Macaronesian Scarabaeinae as a consequence of the low number of species of lection regimens and habitat alterations promoted by their actions in the Mediterranean basin during prehistoric and historic periods F I G U R E 4 Percentage of explained deviance for the four considered hypotheses on the variation in species richness of the three main dung beetle groups on studied islands.Percentage values of anthropogenic (grey), energetic/climatic (red), island biogeography (blue) and environmental heterogeneity (yellow) predictors are included within each rectangular bar.Upper bars represent the explained deviance of these variables while bottom bars represent their explanatory capacity on the residuals of a GLM regression between species richness and the number of published papers on dung beetle fauna in each island.
explaining variation in dung beetle species richness in Mediterranean and Macaronesian islands are those related to anthropogenic actions and, to a lesser extent, those variables that try to represent the species' capacity to colonize islands, such as those classically considered in island biogeography.The years since first human colonization and, to a lesser extent, island area or island heterogeneity seem to be the most relevant variables.Thus, our results suggest a colonization process favoured by the movements of humans and their associated domestic fauna.This process of colonization from the mainland must have occurred repeatedly in view of the low endemicity of these islands.Trophic resources derived from livestock and human faeces, together with frequent human transport of materials and animals, may have facilitated the unintentional immigration of dung beetles on the Mediterranean and Macaronesian islands.Classical biogeographical factors related to area, environmental heterogeneity and distance from the mainland must also have played a role.However, the high dispersal capacity of these species adapted to an ephemeral habitat and the close relationship of this fauna with livestock activities overshadow the role played by these classical factors.The importance of anthropogenic factors in determining the species richness of dung beetles on Mediterranean and Macaronesian islands suggests that these same factors must be promoted if the biodiversity of this fauna and the ecosystem functions it performs are to be preserved.Conservation efforts should thus focus on favouring traditional human actions related to pastoralism and extensive