Potential effects of climate change on the distribution of Scarabaeidae dung beetles in Western Europe

Page 1

ORIGINAL PAPER

Potential effects of climate change on the distribution

of Scarabaeidae dung beetles in Western Europe

E. Dortel • W. Thuiller • J. M. Lobo •

H. Bohbot • J. P. Lumaret • P. Jay-Robert

Received: 18 March 2013 / Accepted: 29 August 2013

Ó Springer Science+Business Media Dordrecht 2013

Abstract Dung beetles are indispensable in pasturelands,

especially when poor efficiency of earthworms and irreg-

ular rainfall (e.g. under a Mediterranean climate) limit pad

decomposition. Although observed and projected species

range shifts and extinctions due to climate change have

been documented for plants and animals, little effort has

focused on the response of keystone species such as the

scarab beetles of dung beetle decomposers. Our study aims

to forecast the distribution of 37 common Scarabaeidae

dung beetle species in France, Portugal and Spain (i.e.

more than half of the western European Scarabaeidae

fauna) in relation to two climate change scenarios (A2 and

B1) for the period leading to 2080. On average, 21 % of the

species should change in each 50-km UTM grid cell. The

highest faunistic turnover rate and a significant increase in

species richness are expected in the north of the study area

while a marked impoverishment is expected in the south,

with little difference between scenarios. The potential

enrichment of northern regions depends on the achieve-

ment of the northward shift of thermophilous species, and

climate change is generally likely to reduce the current

distribution of the majority of species. Under these con-

ditions, the distribution of resource—i.e. the extent and

distribution of pastures—will be a key factor limiting

species’ responses to climate change. The dramatic aban-

donment of extensive grazing across many low mountains

of southern Europe may thus represent a serious threat to

dung beetle distribution changes.

Keywords Dung beetles Á Scarabaeidae Á Climate

change Á Land use change Á Europe

Introduction

Although sudden disruptions to biogeochemical cycles and

land use changes have severe impacts on biodiversity, an

increasing number of studies show that climatic change has

already affected many species and ecosystems during the

last century (McCarty 2001; Parmesan and Yohe 2003;

Root et al. 2003; Parmesan 2006). The latitudinal and

altitudinal distributions of terrestrial biota prove that tem-

perature may play an important role in the response of

many organisms. Changes in species distributions modify

interspecific interactions and could lead to the extinction of

more specialized species (Hughes 2000; McCarty 2001;

Parmesan 2006), whereas generalist, thermophilous or

invasive species could expand in new habitats (Thuiller

2007). Marked species turnover could strongly affect eco-

system functioning. In this way, the disruption of regula-

tion processes provided by species (e.g. carbon storage)

will exacerbate the anthropogenic global warming by

Electronic supplementary material The online version of this

article (doi:10.1007/s10841-013-9590-8) contains supplementary

material, which is available to authorized users.

E. Dortel Á H. Bohbot Á J. P. Lumaret Á P. Jay-Robert

UMR CEFE 5175, 1919 route de Mende,

34293 Montpellier Cedex 5, France

W. Thuiller

Laboratoire d’Ecologie Alpine, UMR-CNRS 5553, Université

Joseph Fourier, BP 53, 38041 Grenoble Cedex 9, France

J. M. Lobo

Museo Nacional de Ciencias Naturales CSIC, c/José Gutierrez

Abascal 2, 28006 Madrid, Spain

J. P. Lumaret Á P. Jay-Robert (&)

Université Montpellier 3, Route de Mende,

34199 Montpellier Cedex 5, France

e-mail: pierre.jay-robert@univ-montp3.fr

123

J Insect Conserv

DOI 10.1007/s10841-013-9590-8

Page 2

positive retroaction (Friedlingstein 2008; Heimann and

Reichstein 2008).

In terrestrial ecosystems, where C and N are stored in

perennial vegetation and in the soil, mineralisation by

decomposers (microorganisms and fungi) is the last step of a

complex decomposition process carried out by detritivores

(invertebrates). Among detritivores, some keystone species

have a strong impact on decomposition (Hättenschwiler

et al. 2005). In grazed ecosystems, for example, dung beetles

(Scarabaeoidea: Scarabaeidae, Aphodiidae, Geotrupidae)

play an essential role in tearing up faeces, burying, and

sowing them with microorganisms (Hanski and Cambefort

1991). Their absence can cause major problems and, in

Australia, the introduction about twenty Scarabaeidae spe-

cies from the Mediterranean basin and South Africa was

required to bury the dung of the domestic stock, and thereby

benefit pasture production and reduce the numbers of pest

flies (Doube et al. 1991).

In Europe, both species richness and endemism of

Scarabaeidae are concentrated around the Mediterranean

basin (Lumaret and Lobo 1996; Lobo and Martın-Piera

2002; Lobo et al. 2002) and one may expect that climate

warming could induce a global northward shift of the

distribution of most of these species. Nevertheless, during

the twentieth century, a significant rarefaction of Scara-

baeidae was observed in Europe (Lumaret 1990; Biström

et al. 1991; Barbero et al. 1999; Roslin 1999; Lobo 2001;

Lobo et al. 2001; Carpaneto et al. 2007). This worrying

trend should be attributed to the drastic changes observed

in land use and farm practices.

Consequently, the magnitude of potential changes in

species distribution driven by climate change needs to be

evaluated both at medium latitudes where the burying

activity of Scarabaeidae could become crucial if aridity

prevents the activity of mesophilous dung beetles (Aph-

odiidae, Geotrupidae) and earthworms, and in the

southernmost extremities of the continent (Iberian, Italian

and Balkan peninsula) where the rapid replacement of

European species by a more thermophilous fauna is

questionable. Our ultimate goal is actually not to pre-

cisely map the future distribution of species (an

unreachable goal considering the complexity of biotic

and abiotic interactions; Duncan et al. 2009) but to

examine the probable effect of climatic change on the

current species distribution in order to evaluate the pos-

sible variations in species richness and regional faunas.

By comparing the consequences of extreme A2 and B1

scenarios developed by the IPCC for the twenty-first

century (Nakicenovic and Swart 2000), we would also

determine if the differences in climatic conditions related

to these alternative economic scenarios could have very

different consequences for this family of beetles.

Methods

Distribution data

The study area extends from 36°N to 51°N across the

Iberian Peninsula and France, i.e. the main centre of spe-

cies richness for European dung beetles and contiguous

temperate extension. Data were extracted from the French

Scarabaeoidea Laparosticti database (Lumaret 1990; Lobo

et al. 1997a; http://inpn.mnhn.fr) and the Iberian Scara-

baeidae database (Lobo and Martın-Piera 1991; BANDA-

SCA free on www.biogeografia.org). The French database

includes more than 42,000 records (compiled from 1762 to

2006) for nearly 190 species (49 Scarabaeidae) observed in

mainland France and on Corsica. The Iberian database

includes more than 15,900 records (from 1872 to 2001) for

55 Scarabaeidae species observed in the Iberian Peninsula

and Balearic islands. The 61 Scarabaeidae species with

records later than 1950 (in the French or Iberian database)

were selected. For each species, we extracted and mapped

the records from 1950 with ArcGis 9.1 (ESRI Corp.,

Redlands, California).

Estimation of adequately surveyed areas

Adequately surveyed cells were discriminated to increase

the reliability of the used absence data. The study area

(France and the Iberian Peninsula), was divided into

10 9 10 km UTM grid cells, to give a total of 11,995

cells. Because the number of species per cell depends on

environmental conditions, the main biogeographic regions

were defined according to the classification proposed by

the European Topic Centre on Biological Diversity

(http://biodiversity.eionet.europa.eu/) for France and by

Lobo and Martın-Piera (2002) for the Iberian Peninsula.

Each 100 km2 UTM grid cell was attributed to a bio-

climatic subregion (Figs. 1, 2). For each bioclimatic

subregion, well surveyed cells were identified. The ade-

quacy of sampling in each cell was determined by a

negative exponential function relating the number of

species (Sr) to the number of records from 1950 (r)

(Soberón and Llorente 1993):

Sr ¼ Smax 1 À exp Àbr

¼ a 1 À exp Àbr

Þ/b with Smax ¼ a/b

where Smax, the asymptote, was the estimated total number

of species per cell, a corresponded to the increase rate at

the beginning of the inventory (r = 0) and b characterized

the shape of the curve. Because the value of the asymptote

depended on the subregion, the function was fitted for each

subregion by the Quasi-Newton method (Jiménez-Valverde

and Hortal 2003). Then, we calculated the number of

J Insect Conserv

123

Page 3

Fig. 1 Distribution of the five

biogeographic climatic areas of

France (A Alpine, B Atlantic,

C Continental, D Corsica,

E Mediterranean)

Fig. 2 Map of the six principal

physioclimatic subregions of the

Iberian Peninsula determined by

Lob and Martın-Piera (2002)

(A Eurosiberian, B Montane,

C North Plateau, D South

Plateau, E East Mediterranean,

F West Mediterranean,

G Balearic Islands)

J Insect Conserv

123

Page 4

records required for a rate of species increment B0.01 (less

than one added species for each 100 database records):

r0.01 ¼ 1/b ln 1 þ b/0.01

The cells for which the number of records r was higher

than r0.01 were considered well surveyed and kept for the

species distribution modelling part. This analysis showed

that 188,100 km2 cells out of 11,995 can be considered

well surveyed (1.57 %; Table 1; Fig. 3). These cells were

widely distributed along the latitudinal gradient and most

of them were located in hilly subregions (Mediterranean,

South Plateau, Continental and Alpine) characterized by a

high climatic and edaphic heterogeneity. France was

slightly more prospected than Spain. For Corsica, well

inventoried cells were deduced from the results (r0.01)

obtained for the Mediterranean subregion. All calculations

were performed with STATISTICA 6.0 (Stat Soft 2001).

Explanatory variables

Nine climatic variables were used as predictors: annual

spring rainfall (mm; RAINSP), annual summer rainfall

(mm; RAINSU), annual winter rainfall (mm; RAINW),

annual autumnal rainfall (mm; RAINA), annual mean tem-

perature (°C; T), minimum and maximum monthly mean

temperature (°C; TMIN and TMAX respectively), annual

mean net radiation (MJ m-2; RAD) and mean annual real

evapotranspiration/potential

evapotranspiration

ratio

(EVA). Considering that Scarabaeidae species are soil-dig-

ger beetles with a complete endogenous larval development,

we also used 7 edaphic variables: 4 variables reflecting the

texture of soil (fine, medium fine, medium, coarse) and 3

variables reflecting the soil water regime (soil very dry, dry,

or wet). The inclusion of these non-climatic variables in the

modelling process may help to determine the true indepen-

dent contribution of climatic variables, thus enhancing the

reliability of future projected distributions (Luoto and Hei-

kkinen 2008; Real et al. 2010; Aragón et al. 2010). The

edaphic data came from the Soil Geographical Database of

Eurasia provided by the European Commission (http://

eusoils.jrc.ec.europa.eu/). We considered that soil charac-

teristics remain constant for the studied period.

The climatic data were obtained from the Climate

Research Unit (CRU data center). The current values cor-

responded to 1961–1990 averages. The forecasted values

corresponded to 2070–2099 averages (Hadley Climate

Model) established according to both extreme scenarios

(A2 and B1) developed by the IPCC Special Report

(Nakicenovic and Swart 2000). Scenario A2 considers a

very heterogeneous world with a strong population growth,

slow economic and technological developments and a

reinforcement of current regional inequalities. Scenario B1

assumes a converging world with global solutions oriented

at durability, the decline in inequality between regions, the

decrease in population from 2050. Thus, 3.4 and 1.8 °C

worldwide increase in temperature are forecasted for A2

and B1 scenarios, respectively. For each variable, the mean

current and forecasted values were calculated for each

100 km2 UTM cell (performed with ArcGis 9.1).

Forecasting models

To summarize the relationship between the presence/

absence of species in UTM-cells (response variable) and all

Table 1 Asymptotic relationship between number of species and sampling effort for the biogeographic areas and physioclimatic subregions

Number of records

Total UTM-cells

Smax

r0.01

Well surveyed cells

France

Alpine

4,159

464

0.909

35.47

1.206

0.034

43.577

22 (4.74 %)

Atlantic

3,571

3,031

0.968

55.45

1.109

0.020

54.931

4 (0.13 %)

Continental

8,121

1,775

0.936

40.20

1.206

0.03

46.21

33 (1.86 %)

Mediterranean

11,892

634

0.901

53.60

1.340

0.025

50.11

63 (9.94 %)

Corsica

1,477

126

0.943

37.60

1.504

0.04

40.236

5 (3.97 %)

Total

29,220

6,030

127 (2.11 %)

Iberian peninsula

Eurosiberian

337

487

0.958

23.092

1.247

0.054

34.376

Montane

1,832

1,005

0.953

34.125

1.365

0.04

40.236

12 (1.19 %)

North Plateau

2,376

1,200

0.954

29.104

1.397

0.048

36.622

10 (0.83 %)

South Plateau

6,039

1,946

0.930

36.459

1.349

0.037

41.826

34 (1.75 %)

East Mediterranean

199

841

0.984

37.933

1.138

0.03

46.21

West Mediterranean

994

703

0.947

30.909

1.360

0.044

38.327

5 (0.71 %)

Balearic Islands

120

0.974

22.222

1.200

0.054

34.376

Total

11,897

6,277

61 (0.9 %)

J Insect Conserv

123

Page 5

the formerly mentioned explanatory variables, we used

three kind of modelling techniques implemented into the

BIOMOD library (Thuiller 2003; Thuiller et al. 2009)

under the R software: Generalised Linear Model (GLM),

Generalised Additive Model (GAM) and Boosted Regres-

sion Trees (BRT). GLMs were fitted using linear and

quadratic functions. GAMs were fitted using 4 degree of

smoothing. We used a forward stepwise selection for var-

iable selection into the GLM and GAM using the AIC

criteria. We used 100-fold cross-validation to select the

optimal number of trees for the BRT with an interaction-

depth of 2.

The models were run for each species, using all avail-

able presence data but only absences coming from cells

previously considered as well surveyed. By using these

absences in both calibration and evaluation processes we

try to minimize the effects of false absences. We assume

that the estimated parameters of the so obtained predictive

functions are capable of reflecting the effect of climate on

current distributions, and that these effects remain similar

under the future climatic conditions. The obtained contin-

uous probabilities were converted into binary presence/

absence data, selecting the recommended threshold that

maximized the percentage of presences and absences cor-

rectly predicted (see Jiménez-Valverde and Lobo 2007).

The models were finally performed for 40 species (35

species common to France and Iberian Peninsula, 1 species

strictly French and 4 species strictly Iberian) observed in at

least 40 presence cells and 40 absence cells considered as

well surveyed.To be reliable, predictions must be validated

using independent data but such data did not exist in our case.

We thus randomly divided the original dataset into two sub-

sets: the calibration data (70 % of total presence/absence

data) used to calibrate the models and the evaluation data

(30 %) used to examine the reliability of the predictions.

The best model among the outputs of the three model-

ling techniques was determined, for each species, by using

the area under the Receiver Operating Characteristic curve

(AUC Hanley and McNeil 1982) on the evaluation data.

This curve represents the fraction of presences correctly

predicted (sensitivity) as a function of commission errors

(1-specificity) for a range of thresholds, being used as a

standard discrimination measure. AUC method cannot be

used when the comparisons differ in the relative occurrence

area (the ratio between the species extent and the whole

extent of the region of study), as occur in the case of dis-

tributional models carried out at the same extent for species

differing in range size (Lobo et al. 2008). That was not the

case here as AUC was used to compare models built with

the same species data. Furthermore, we assume that the use

Fig. 3 Map of 1,199,510 km

UTM-cells (adequately sampled

cells in black)

J Insect Conserv

123

Page 6

of absence data selected from well surveyed cells should

increase the capacity of minimizing both commission and

omission errors in the evaluation process. We overlaid both

observed individual species maps and predicted estima-

tions in order to obtain two species richness maps at a

10 9 10 km cell resolution, calculating the correlation

between observed and predicted species richness in the

well surveyed cells as a measure of the general accuracy of

predictions.

Changes in species’ distribution

The number of species by 50 9 50 km UTM grid cells

(pools of 25,100 km2 UTM grid cells) was calculated for the

twentieth century and for the twenty-first century according

to A2 and B1 scenarios. To do that we overlap each one of

the individual predictions being the species considered

present in a 50 9 50 km cell if it is predicted as present in

any one of their constituent 10 9 10 km cells. The use of

50 9 50 km UTM grid cells should allow to reduce the bias

induced by the relative lack of precision of forecasted spe-

cies distributions. b-diversity was used to evaluate the spe-

cies turnover between the end of the twentieth century and

the end of the twenty-first century. The temporal turnover

was estimated with the Wilson and Shmida’s b-diversity

index (1984) rescaled at a 0–100 range:

b ¼ 100 Ã Sg þ Sl

(

)

/ S20 þ S21

where Sg is the number of species gained between the

dates, Sl the number of species lost, S20 the predicted

number of species for the twentieth century and S21 the

predicted number of species for the twenty-first century.

b = 0 if fauna does not change between the dates, b = 50

if there is as many gained and lost species as shared spe-

cies, b = 100 if fauna completely changes.

Results

Choice of models

The mean of the obtained AUC values is 0.898 ± 0.017

(±95 % confidence interval; minimum = 0.783, maxi-

mum = 0.997; see Table 2) which may be considered a high

value when compared to those obtained in similar studies

(Elith et al. 2006). The Pearson correlation value between

observed and predicted species richness in the well surveyed

10 9 10 km cells is positive and statistically significant

(r = 0.66, n = 188, p\ 0.001). In spite of high AUC val-

ues, three species (O. stylocerus, S. semipunctatus and S.

typhon) were excluded from the analysis because of too

peculiar forecasted distributions. The present work finally

concerned 37 Scarabaeidae species (maps in Annex).

Future forecasted distributions

Important changes in distribution appeared for both sce-

narios: the mean change (sum of gain and loss) reached

70 % for A2 and 63 % for B1 (Table 3) and only four

species did show a change in distribution lower than 30 %

Table 2 AUC values obtained in the evaluation data (30 % of total

data) for all the considered species

Species

GAM

BRT

GLM

Bubas bison

0.962

0.972

0.964

Bubas bubalus

0.959

0.935

0.942

Caccobius schreberi

0.836

0.866

0.798

Cheironitis hungaricus

0.851

0.777

0.630

Copris hispanus

0.947

0.958

0.954

Copris lunaris

0.840

0.878

0.841

Euoniticellus fulvus

0.832

0.846

0.837

Euoniticellus pallipes

0.966

0.911

0.941

Euonthophagus amyntas

0.899

0.912

0.903

Euonthophagus gibbosus

0.831

0.828

0.755

Gymnopleurus flagellatus

0.871

0.891

0.869

Gymnopleurus sturmi

0.687

0.783

0.606

Onitis belial

0.947

0.945

0.927

Onitis ion

0.902

0.879

Onthophagus baraudi

0.976

0.985

0.824

Onthophagus coenobita

0.835

0.871

0.840

Onthophagus emarginatus

0.842

0.891

0.837

Onthophagus fracticornis

0.917

0.943

0.898

Onthophagus furcatus

0.821

0.845

0.821

Onthophagus grossepunctatus

0.791

0.816

0.782

Onthophagus illyricus

0.741

0.818

0.799

Onthophagus joannae

0.883

0.873

0.838

Onthophagus latigena

0.944

0.973

0.956

Onthophagus lemur

0.863

0.892

0.852

Onthophagus maki

0.829

0.888

0.792

Onthophagus nuchicornis

0.864

0.938

0.912

Onthophagus opacicollis

0.898

0.889

0.901

Onthophagus ovatus

0.915

0.862

0.911

Onthophagus punctatus

0.927

0.897

0.882

Onthophagus ruficapillus

0.767

0.819

0.792

Onthophagus similis

0.850

0.838

0.861

Onthophagus stylocerus

0.919

0.941

0.861

Onthophagus taurus

0.875

0.870

0.808

Onthophagus vacca

0.814

0.956

0.914

Onthophagus verticicornis

0.853

0.890

0.872

Scarabaeus laticollis

0.846

0.850

0.842

Scarabaeus sacer

0.958

0.963

0.946

Scarabaeus semipunctatus

0.995

0.997

0.869

Scarabaeus typhon

0.850

0.851

0.712

Sisyphus schaefferi

0.808

0.871

0.785

J Insect Conserv

123

Page 7

of their current distribution (O. vacca, O. illyricus, O.

punctatus and O. ovatus). This change in distribution

mainly corresponded to a northward shift: on average 1.16°

for A2 and 1.01° for B1 (max. *5° for E. amyntas).

Results obtained for the two scenarios were not very dif-

ferent (slightly higher changes for A2). The geographical

boundaries of our study might lead to an underestimation

of the latitudinal shift. Moreover, an altitudinal shift (dif-

ficult to detect with 10-km UTM grid cells) was also

observed (e.g. O. ovatus, O. verticicornis).

Two parameters were important to distinguish: (1) the

loss of present favourable cells and (2) the net balance

Table 3 Estimated changes in the distribution of species according to the A2 and B1 IPCC scenarii

Species

Model Estimated current

range size (nb cells)

Scenario A2

Scenario B1

Loss (%) Gain (%) DLon

DLat

Loss (%) Gain (%) DLon

DLat

Bubas bison

BRT

3,658

110.69

1.4

2.74

0.14

81.25

0.99

2.06

Bubas bubalus

GAM

3,789

9.77

83.5

0.62

1.03 12.83

78.65

0.35

0.89

Caccobius schreberi

BRT

8,556

32.82

20.47

-0.03

0.51 26.05

23.19

0.39

0.72

Cheironitis hungaricus

GAM

3,178

76.49

12.24

-0.41

1.73 51.1

4.37

0.57

1.05

Copris hispanus

BRT

3,844

0.05

103.46

1.12

2.3

0.03

71.96

0.79

1.66

Copris lunaris

BRT

6,676

84.17

14.92

1.94

3.1

71.39

18.05

1.79

2.5

Euoniticellus fulvus

BRT

8,560

6.99

25.34

0.58

0.39 15.02

26.13

0.61

0.36

Euoniticellus pallipes

GAM

3,115

137.98

1.51

2.64

0.51

95.76

1.13

1.98

Euonthophagus amyntas

GAM

5,086

75.01

56.04

3.15

5.44 73.48

60.48

2.87

4.8

Euonthophagus gibbosus

GAM

2,943

23.82

96.47

-3.08 -0.87

6.56

140.47

-1.65 -0.06

Gymnopleurus flagellatus

BRT

4,369

29.96

28.66

1.37 22.98

28.36

0.92

1.21

Gymnopleurus sturmi

BRT

3,592

0.22

149.33

1.95

2.94

0.33

118.62

1.95

2.58

Onitis belial

BRT

3,447

3.16

42.3

0.66

0.39 12.01

40.59

0.75

0.2

Onitis ion

BRT

2,818

0.71

44.18

0.75

0.51

2.66

40.92

0.9

0.44

Onthophagus baraudi

BRT

310

74.52

1.96

0.54 62.58

1.48

0.28

Onthophagus coenobita

BRT

6,074

64.82

11.46

0.88

0.79 43.48

10.9

0.87

0.86

Onthophagus emarginatus

BRT

6,147

19.44

32.44

0.66

0.79 18.32

31.02

0.55

0.87

Onthophagus fracticornis

BRT

7,833

75.32

0.79

0.94 -0.2

56.49

0.84

1.08 -0.76

Onthophagus furcatus

BRT

5,975

1.82

31.62

0.89

1.19

1.77

29.84

0.84

1.1

Onthophagus grossepunctatus BRT

5,294

27.92

31.6

1.01

1.05 26.14

30.9

1.01

1.19

Onthophagus illyricus

BRT

6,344

14.52

7.52

0.49

0.49 14.56

7.22

0.43

Onthophagus joannae

BRT

6,204

56.88

5.14

0.93 -0.2

5.08

0.55 -0.41

Onthophagus latigena

BRT

1,675

18.57

64.96

1.12

1.06 16.96

1.1

1.09

Onthophagus lemur

BRT

6,298

55.83

12.38

1.54

0.77 48.89

12.67

1.1

0.55

Onthophagus maki

BRT

4,519

19.36

48.59

0.51

1.11 21.04

53.44

0.87

1.59

Onthophagus nuchicornis

BRT

5,694

16.05

15.1

-0.13

0.17 18

13.12

-0.08

0.25

Onthophagus opacicollis

GLM

5,575

0.04

69.29

1.25

1.12

1.42

49.69

0.89

0.68

Onthophagus ovatus

GAM

6,688

15.68

12.87

0.42

0.43

9.46

18.53

0.02

0.18

Onthophagus punctatus

BRT

4,912

4.05

19.89

0.98

0.92

3.14

25.47

1.07

1.13

Onthophagus ruficapillus

GLM

6,503

3.48

42.98

1.43

1.12

4.81

43.4

1.29

1.13

Onthophagus similis

GLM

7,859

43.87

10.89

1.25

1.96 41.9

13.3

1.94

1.96

Onthophagus taurus

GAM

8,521

34.26

1.38

1.22

0.04

32.75

1.35

1.12

Onthophagus vacca

BRT

10,136

0.01

14.7

0.67

0.59

0.21

13.21

0.62

0.56

Onthophagus verticicornis

BRT

5,708

61.63

13.58

1.26

0.42 53.01

13.89

0.97

0.49

Scarabaeus laticollis

BRT

4,322

23.3

43.48

0.32

1.18 22.61

38.89

0.25

0.94

Scarabaeus sacer

BRT

3,792

13.16

72.49

1.25

1.43 14.06

58.18

1.05

1.12

Sisyphus schaefferi

BRT

2,936

85.49

13.28

0.65

0.75 74.9

18.15

0.43

0.76

Loss and Gain are expressed in % of estimated current range size

DLon and DLat (decimal degrees) = difference in central location of each species’ longitude and latitude between the current and the forecast

estimated distributions (positive = shift toward East/North; negative = shift toward West/South)

J Insect Conserv

123

Page 8

between gain and loss of cells. The loss of current

favourable areas could lead to the disappearance of popu-

lations and, consequently, could weaken the capacity of

species to face global warming by reaching new favourable

territories. The net gain of favourable territories should be

considered with caution because the actual presence of

species in such new territories depends on many factors

(landscape structure, population dynamics, etc.).

Three categories of species can be distinguished (values

for A2 scenario; Table 3; Fig. 4):

(1) Species with no clear trend: twelve species could lose

or gain no more than 20% of their predicted

distribution (C. schreberi, E. fulvus, E. amyntas, G.

flagellatus, O. emarginatus, O. grossepunctatus, O.

illyricus, O. nuchicornis, O. ovatus, O. punctatus, O.

vacca, S. laticollis). This relatively low change was

associated with a reduced loss of present favourable

surfaces, except for E. amyntas (75 %).

(2) Looser species: ten species could lose at least one

third of their distribution in the studied area. This

global decrease corresponded to a loss of current

favourable cells comprised between 44 % (O. similis)

and 86 % (S. schaefferi). For seven of these species

(C. lunaris, O. coenobita, O. fracticornis, O. joannae,

O. lemur, O. similis, O. verticicornis), an expansion

north of the studied area is possible while for C.

hungaricus, S. schaefferi (Mediterranean) and O.

baraudi (alpine endemic) the disappearance of pres-

ent favourable sectors could not be compensated and

Fig. 4 Change in distribution between present and future (A2

scenario). Each circle represents a species (positive values for gain,

negative values for loss)

Fig. 5 The distribution of

estimated dung beetle species

richness at the present time (a),

and at the end of the twenty-first

century according to the two

climatic scenarii A2 (b) and B1

(c)

J Insect Conserv

123

Page 9

the species could lose more than 2/3 of their current

distribution.

(3) Winner species: fifteen species could increase their

current predicted distribution in the area by more than

29 % (B. bison, B. bubalus, C. hispanus, E. gibbosus,

E. pallipes, G. sturmi, O. belial, O. ion, O. furcatus,

O. latigena, O. maki, O. opacicollis, O. ruficapillus,

O. taurus, S. sacer). For E. gibbosus, O. latigena and

O. maki this increase was accompanied by a loss of

more than 15 % of present favourable cells, whereas,

at the opposite, the net gain exceeded 100 % for B.

bison, C. hispanus, E. pallipes and G. sturmi.

Changes in biodiversity

The map of current predicted species richness deriving from

the overlay of individual distribution maps (Fig. 5a) showed

a similar spatial pattern that the estimations of species

richness distribution previously obtained by another way

(Lobo and Martın-Piera 2002; Lobo et al. 2002). By com-

parison with the present, the species richness maps corre-

sponding to the scenarios A2 and B1 showed a general

increase in diversity in France, and a significant decrease in

the centre of the Iberian Peninsula for the end of the twenty-

first century (Fig. 5b, c). A very similar regional species

turnover was expected under both scenarios (Fig. 6a, b): on

average 21% of the fauna could change in each

50 9 50 km UTM grid cell (20.8 % with A2; 20.9 % with

B1; t = -0.22, df = 1,116, p = 0.82) with highest values

along the eastern French boundary and in the north-western

half of this country (up to 68 % for A2 and 55 % for B1).

Discussion

In the Western Palaearctic region, the Scarabaeidae fauna

contains 162 species split into 12 genera (Cabrero-San˜udo

and Lobo 2003). The present study focused on 37 of the 61

species from the Iberian Peninsula and France. These spe-

cies represent eleven genera, are all widely distributed in the

studied area (Lumaret and Lobo 1996) and one can assume

that the 15° latitudinal gradient/4,000 m altitudinal gradient

covered by our study captured a well part of their environ-

mental niche. The choice of these species, primarily

imposed by the necessity to have a large number of obser-

vations, allowed us to avoid the problems related to micro-

endemism and historical contingencies (Guisan 2003; Ara-

újo et al. 2008). The widespread distribution of the 37

studied species proved that they have been able to disperse

and, consequently, one may expect that they should be able

to respond to future climate change by modifications in their

distribution. Because local abundance and distribution range

are generally correlated in dung beetles (Lobo 1993), one

may consider that our work dealt mainly with core species

that constitute the bulk of local assemblages (Hanski 1991).

Our results forecast a general northward shift of Scara-

baeidae with a related increase in species richness at

intermediate and northernmost latitudes. In parallel, an

altitudinal increase—hard to depict with 100 km2 UTM-

cells—may operate for some species. The static nature of

our modelling, as well as the scale and resolution consid-

ered, does not allow to derive more detailed conclusions,

especially about the effect of climate change on the

modalities of species coexistence (Guisan and Thuiller

2005; Duncan et al. 2009).

Fig. 6 The spatial distribution

of estimated dung beetle species

turnover for the end of the

twenty-first century according

to the A2 (a) and B1

(b) climatic scenarii

J Insect Conserv

123

Page 10

This shift was expected because Scarabaeidae constitute

the main thermophilous group of dung beetles (Lumaret

and Kirk 1991; Lobo and Martın-Piera 2002; Lobo et al.

2002), and some recent empirical evidences support it.

Onitis belial, a species formerly restricted to some places

of the Mediterranean seashore in France (Paulian and Ba-

raud 1982), was recently observed at an altitude of 900 m

at the eastern end of the Pyrenees (Jay-Robert unpubl.).

The mean northward shift in the predicted distributions

reached ca. 13 km per decade under the A2 scenario and,

for some species, the northern range boundary could move

around 75 km per decade. An extensive analysis of recent

changes in European butterfly communities suggests that

these predictions are conservative (Devictor et al. 2012).

Nevertheless, differences between A2 and B1 scenarios

were slight and, for both scenarios, the mean change in

regional compositions (50 9 50 km UTM grid cells) could

exceed 20 % of species. Several core species—that can

constitute more than 75 % of local dung beetle assem-

blages (Lumaret and Kirk 1991; Errouissi et al. 2004; Jay-

Robert et al. 2008a)—would lose a large part of their

current distribution (e.g. O. fracticornis, O. joannae, O.

lemur or O. similis). Even if little is known about the

respective functional efficiency of species in dung removal

(Rosenlew and Roslin 2008), one may fear that the local

disappearance of such very common species may induce a

depletion of ecosystem functioning (Cardinale et al. 2006).

At intermediate latitude the northward shift of

mesophilous species might be compensated by the arrival

of more thermophilous ones. This replacement requires a

very rapid response of populations (Parmesan and Yohe

2003), one that is questionable for low-fecundity taxa like

Bubas sp., Copris sp. or Gymnopleurus sp. For most spe-

cies, the potential gain of new territories would be con-

comitant to a significant loss of present habitats, and

population dynamics could be weakened by this rapid

turnover of favourable areas (Keith et al. 2008). Addi-

tionally, the migration towards these new favourable areas

would require a good connectivity among pastured habitats

(especially along the French Atlantic coast and the Rhone

valley). The significance of the connectivity will depend on

species characteristics. In Aphodiidae—the dominant

group of dung beetles in northern Europe (Cabrero-San˜udo

and Lobo 2003)—pasture specialist species had lower

migration abilities than generalist species (Roslin 2000;

Roslin and Koivunen 2001). Unfortunately, available sce-

narios forecast the continuation of the polarisation of cattle

breeding during the twenty-first century, with intensifica-

tion in favourable areas versus abandonment in harsh

regions like southern mountain ranges (Schröter et al.

2005). This polarisation and several concomitant changes

in practices (the stalling of cattle, the increasing using of

parasiticides toxic for beetles…), that begun in the middle

of the twentieth century, were invoked to explain the

decrease in dung beetle diversity which was already

observed from northern latitudes to the Mediterranean

region (Lumaret 1990; Biström et al. 1991; Barbero et al.

1999; Roslin 1999; Lobo 2001; Lobo et al. 2001; Carpa-

neto et al. 2007). Although Scarabaeidae appear to be less

sensitive to habitat heterogeneity than Aphodiidae or

Geotrupidae in southern Europe (Lobo et al. 1997b; Lobo

and Martın-Piera 1999), the expansion of ungrazed wooded

habitats could also severely limit the movements of beetles

(Kadiri et al. 1997; Jay-Robert et al. 2008c) and exacerbate

the role of barrier naturally played by longitudinal Euro-

pean mountain ranges (Jay-Robert et al. 1997).

In the Iberian Peninsula, separated from North Africa by

the Strait of Gibraltar, the settling of a more thermophilous

fauna would be probably more difficult and the shift in the

distribution of mesophilous species should induce a sig-

nificant drop in species richness. An early emergence in

spring could allow insects to maintain local populations

and limit the decline in species richness (Stefanescu et al.

2003), but nowadays adult Scarabaeidae have a typical

spring-summer period of activity everywhere in Europe

(Wassmer 1994; Jay-Robert et al. 2008a, b).

The extent of the faunistic changes forecasted in our

study (and the insignificance of differences between A2

and B1 scenarios) should force breeding industry and

conservationists to collaborate on a win–win strategy

which associates grazing network (e.g. Natura2000 areas)

and agroecological guidelines (which still remain to be

designed). In that case, Scarabaeidae should succeed in

adjusting their distribution to climate change and com-

pensate for the local risk of rarefaction of mesophilous

Aphodiidae and Geotrupidae species.

Acknowledgments We are very grateful to John Thompson (UMR

5175 CEFE, Montpellier, France) who revised the English version of

the manuscript. ED received financial support from the Agence Na-

tionale de la Recherche contract ANR-05-BDIV-014. WT received

support from European Commission’s FP6 ECOCHANGE (Chal-

lenges in assessing and forecasting biodiversity and ecosystem

changes in 17 Europe, No 066866 GOCE) project.

References

Aragón P et al (2010) The contribution of contemporary climate to

ectothermic and endothermic vertebrate distributions in a glacial

refuge. Glob Ecol Biogeogr 19:40–49

Araújo M et al (2008) Quaternary climate changes explain diversity

among reptiles and amphibians. Ecography 31(1):8–15

Barbero E et al (1999) Dung beetle conservation: effects of habitat

and resource selection (Coleoptera: Scarabaeoidea). J Insect

Conserv 3:75–84

Biström O et al (1991) Abundance and distribution of coprophilous

Histerini (Histeridae) and Onthophagus and Aphodius (Scara-

baeidae) in Finland (Coleoptera). Entomologica Fennica 2:53–66

J Insect Conserv

123

Page 11

Cabrero-San˜udo FJ, Lobo JM (2003) Estimating the number of

species not yet described and their characteristics: the case of

Western Palaearctic dung beetle species (Coleoptera, Scarabae-

oidea). Biodivers Conserv 12:147–166

Cardinale BJ et al (2006) Effects of biodiversity on the functioning of

trophic groups and ecosystems. Nature 443:989–992

Carpaneto GM et al (2007) Inferring species decline from collection

records: roller dung beetles in Italy (Coleoptera, Scarabaeidae).

Divers Distrib 13:903–919

Devictor V et al (2012) Differences in the climatic debts of birds and

butterflies at a continental scale. Nat Clim Change 2:121–124

Doube BM et al (1991) Native and introduced dung beetles in

Australia. In: Hanski I, Cambefort Y (eds) Dung beetle ecology.

Princeton University Press, Princeton, pp 255–282

Duncan RP et al (2009) Do climate envelope models transfer? A

manipulative test using dung beetle introductions. Proc R Soc B

267:1449–1457

Elith J et al (2006) Novel methods improve prediction of species’

distributions from occurrence data. Ecography 29:129–151

Errouissi F et al (2004) Composition and structure of dung beetle

assemblages in mountain grasslands of the Southern Alps. Ann

Entomol Soc Am 97(4):701–709

Friedlingstein P (2008) A steep road to climate stabilization. Nature

451:297–298

Guisan A (2003) Simuler la répartition géographique des espe`ces et

de la végétation (ou ‘‘Si DeCandolle avait eu on ordinateur …’’),

Synthe`se. Saussurea 33:79–99

Guisan A, Thuiller W (2005) Predicting species distribution: offering

more than simple habitat models. Ecol Lett 8:993–1009

Hanley JA, McNeil BJ (1982) The meaning and use of the area under a

receiver operating characteristic (ROC) curve. Radiology 143:29–36

Hanski I (1991) North temperate dung beetles. In: Hanski I,

Cambefort Y (eds) Dung beetle ecology. Princeton University

Press, Princeton, pp 75–96

Hanski I, Cambefort Y (1991) Dung beetle ecology. Princeton

University Press, Princeton

Hättenschwiler S et al (2005) Biodiversity and litter decomposition in

terrestrial ecosystems. Annu Rev Ecol Evol Syst 36:191–218

Heimann M, Reichstein M (2008) Terrestrial ecosystem carbon

dynamics and climate feedbacks. Nature 451:289–292

Hughes L (2000) Biological consequences of global warming: is the

signal already apparent? Trends Ecol Evol 15(2):56–61

Jay-Robert P et al (1997) Altitudinal turnover and species richness

variation in European montane dung beetle assemblages. Arct

Alp Res 29(2):196–205

Jay-Robert P et al (2008a) Temporal coexistence of dung-dweller and

soil-digger dung beetles (Coleoptera, Scarabaeoidea) in con-

trasting Mediterranean habitats. Bull Entomol Res 98:303–316

Jay-Robert P et al (2008b) Spatial and temporal variation of mountain

dung beetle assemblages and their relationships with environ-

mental factors (Aphodiinae: Geotrupinae: Scarabaeinae). Ann

Entomol Soc Am 101(1):58–69

Jay-Robert P et al (2008c) Relative efficiency of extensive grazing vs.

wild ungulates management for dung beetle conservation in a

heterogeneous landscape from Southern Europe (Scarabaeinae,

Aphodiinae, Geotrupinae). Biol Conserv 141:2879–2887

Jiménez-Valverde A, Hortal J (2003) Las curvas de acumulación de

especies y la necesidad de evaluar la calidad de los inventarios

biológicos. Revista Ibérica de Aracnologıa 8:151–161

Jiménez-Valverde A, Lobo JM (2007) Threshold criteria for conver-

sion of probability of species presence to either–or presence–

absence. Acata Oecologica 31:361–369

Kadiri N et al (1997) Conséquences de l’interaction entre préférences

pour l’habitat et quantité de ressources trophiques sur les

communautés d’insectes coprophages (Coleoptera : Scarabaeoi-

dea). Acta Oecologica 18(2):107–119

Keith DA et al (2008) Predicting extinction risks under climate

change: coupling stochastic population models with dynamic

bioclimatic habitat models. Biol Lett 4:560–563

Lobo JM (1993) The relationship between distribution and abundance

in a dung beetle community (Col., Scarabaeoidea). Acta

Oecologica 14(1):43–55

Lobo JM (2001) Decline of roller dung beetle (Scarabaeinae)

populations in the Iberian peninsula during the 20th century.

Biol Conserv 97:43–50

Lobo JM, Martın-Piera F (1991) La creación de un banco de datos

zoológico sobre los Scarabaeidae (Coleoptera: Scarabaeoidea)

Ibero-Baleares: una experiencia piloto. Elytron 5:31–37

Lobo JM, Martın-Piera F (1999) Between-group differences in the

Iberian dung beetle species-area relationship (Coleoptera: Scara-

baeidae). Acta Oecologica 20(6):587–597

Lobo JM, Martın-Piera F (2002) Searching for a predictive model for

species richness of Iberian dung beetle based on spatial and

environmental variables. Conserv Biol 16(1):158–173

Lobo JM et al (1997a) Les atlas faunistiques comme outils d’analyse

spatiale de la biodiversité. Annales de la Société Entomologique

de France (N.S.) 33(2):129–138

Lobo JM et al (1997b) Diversity and spatial turnover of dung beetle

(Coleoptera: Scarabaeoidea) communities in a protected area of

south Europe (Don˜ana National Park, Huelva, Spain). Elytron

11:71–88

Lobo JM et al (2001) Diversity, distinctiveness and conservation

status of the Mediterranean coastal dung beetle assemblage in

the Regional Natural Park of the Camargue (France). Divers

Distrib 7:257–270

Lobo JM et al (2002) Modelling the species richness distribution of

French dung beetles and delimiting the predictive capacity of

different groups of explanatory variables (Coleoptera: Scara-

baeidae). Glob Ecol Biogeogr 11:265–277

Lobo JM et al (2008) AUC: a misleading measure of the performance of

predictive distribution models. Glob Ecol Biogeogr 17:145–151

Lumaret JP (1990) Atlas des Scarabéides Laparosticti de France.

Secrétariat Faune—Flore/MNHN, Paris

Lumaret JP, Kirk AA (1991) South temperate dung beetles. In:

Hanski I, Cambefort Y (eds) Dung beetle ecology. Princeton

University Press, Princeton, pp 97–115

Lumaret JP, Lobo JM (1996) Geographic distribution of endemic

dung beetles (Coleoptera, Scarabaeoidea) in the Western Palae-

arctic region. Biodivers Lett 3:192–199

Luoto M, Heikkinen RK (2008) Disregarding topographical hetero-

geneity biases species turnover assessments based on bioclimatic

models. Glob Change Biol 14:483–494

McCarty JP (2001) Ecological consequences of recent climate

change. Conserv Biol 15(2):320–331

Nakicenovic N, Swart R (2000) Emissions scenarios: a special report

of Working Group III of the Intergovernmental Panel on Climate

Change. Cambridge University Press, Cambridge

Parmesan C (2006) Ecological and evolutionary responses to recent

climate change. Annu Rev Ecol Evol Syst 37:637–669

Parmesan C, Yohe G (2003) A globally coherent fingerprint of

climate change impacts across natural systems. Nature

421:37–42

Paulian R, Baraud J (1982) Faune des Coléopte`res de France. II.

Lucanoidea et Scarabaeoidea. Lechevalier, Paris

Real R et al (2010) Species distribution models in climate change

scenarios are still not useful for informing policy planning: an

uncertainty assessment using fuzzy logic. Ecography 33:304–314

Root TL et al (2003) Fingerprints of global warming on wild animals

and plants. Nature 421:57–60

Rosenlew H, Roslin T (2008) Habitat fragmentation and the

functional efficiency of temperate dung beetles. Oikos

117:1659–1666

J Insect Conserv

123

Page 12

Roslin T (1999) Spatial ecology of dung beetles. PhD thesis.

University of Helsinki

Roslin T (2000) Dung beetle movements at two spatial scales. Oikos

91:323–335

Roslin T, Koivunen A (2001) Distribution and abundance of dung

beetles in fragmented landscapes. Oecologia 127(1):69–77

Schröter D et al (2005) Ecosystem service supply and vulnerability to

global change in Europe. Science 310:1333–1337

Soberón M, Llorente B (1993) The use of species accumulation

functions for the prediction of species richness. Conserv Biol

7(3):480–488

Stat Soft (2001) Statistica 6. Tulsa, OK 74104, USA

Stefanescu C et al (2003) Effects of climatic change on the phenology

of butterflies in the northwest Mediterranean Basin. Glob

Change Biol 9:1494–1506

Thuiller W (2003) BIOMOD: optimising predictions of species

distributions and projecting potential future shifts under global

change. Glob Change Biol 9:1353–1362

Thuiller W (2007) Climate change and the ecologist. Nature

448:550–552

Thuiller W et al (2009) BIOMOD—a platform for ensemble

forecasting of species distributions. Ecography 32(3):369–373

Wassmer T (1994) Seasonality of coprophagous beetles in the

Kaiserstuhl area near Freiburg (SW-Germany) including the

winter months. Acta Oecologia 15(5):607–631

Wilson MV, Shmida A (1984) Measuring beta diversity with presence

absent data. J Ecol 72:1055–1064

J Insect Conserv

123