Skip to main content
Download PDF

Thermal niche evolution and geographical range expansion in a species complex of western Mediterranean diving beetles

BMC Evolutionary Biology volume 14, Article number: 187 (2014) Cite this article

Abstract

Background

Species thermal requirements are one of the principal determinants of their ecology and biogeography, although our understanding of the interplay between these factors is limited by the paucity of integrative empirical studies. Here we use empirically collected thermal tolerance data in combination with molecular phylogenetics/phylogeography and ecological niche modelling to study the evolution of a clade of three western Mediterranean diving beetles, the Agabus brunneus complex.

Results

The preferred mitochondrial DNA topology recovered A. ramblae (North Africa, east Iberia and Balearic islands) as paraphyletic, with A. brunneus (widespread in the southwestern Mediterranean) and A. rufulus (Corsica and Sardinia) nested within it, with an estimated origin between 0.60-0.25 Ma. All three species were, however, recovered as monophyletic using nuclear DNA markers. A Bayesian skyline plot suggested demographic expansion in the clade at the onset of the last glacial cycle. The species thermal tolerances differ significantly, with A. brunneus able to tolerate lower temperatures than the other taxa. The climatic niche of the three species also differs, with A. ramblae occupying more arid and seasonal areas, with a higher minimum temperature in the coldest month. The estimated potential distribution for both A. brunneus and A. ramblae was most restricted in the last interglacial, becoming increasingly wider through the last glacial and the Holocene.

Conclusions

The A. brunneus complex diversified in the late Pleistocene, most likely in south Iberia after colonization from Morocco. Insular forms did not differentiate substantially in morphology or ecology, but A. brunneus evolved a wider tolerance to cold, which appeared to have facilitated its geographic expansion. Both A. brunneus and A. ramblae expanded their ranges during the last glacial, although they have not occupied areas beyond their LGM potential distribution except for isolated populations of A. brunneus in France and England. On the islands and possibly Tunisia secondary contact between A. brunneus and A. ramblae or A. rufulus has resulted in introgression. Our work highlights the complex dynamics of speciation and range expansions within southern areas during the last glacial cycle, and points to the often neglected role of North Africa as a source of European biodiversity.

Background

Information on the thermal biology of a species is fundamental to understand its ecology, biogeography and evolution, as species are only capable of tolerating a limited range of climatic conditions. Ambient temperature affects all biological processes [1],[2], especially in ectotherms [3], and is usually assumed to be one of the main determinants of their spatial distribution [4]. However, in most biogeographical studies the thermal tolerance of species is extrapolated exclusively from their current distributions [5], and even when palaeoclimatic or genetic data are considered (as in e.g. [6]–[8]), it is rare for these to be combined with experimental data on the actual physiological tolerance of the study organisms [9],[10]. Despite this, the need for integrative approaches is increasingly being recognised [11]–[13], particularly given the limitations of current distributional data for inferring historical or ecological processes [14],[15].

Here we attempt such an integrative approach in a clade of diving beetles that has diversified and expanded its range in the western Mediterranean region during the late Pleistocene, the Agabus brunneus complex [16]. Previous work has revealed that thermal tolerance is a good predictor of geographical range extent in these beetles, in which more widespread species have wider thermal windows than their narrow-range relatives [17]. Two species of the complex have partly overlapping distributions in southwest Europe and North Africa, whilst the third is confined to the islands of Corsica and Sardinia [16].

Traditionally, the study of the recent evolutionary history of the European fauna and flora has largely considered the direct effect of the Pleistocene glaciations, particularly the recolonization of previously glaciated areas from unglaciated refugia and the genetic changes resulting from such range movements [18]–[20]. In most cases, unglaciated areas are simply seen as refugia for northern species, little attention being paid to evolutionary and biogeographical processes in them, other than those which affected these species [21],[22]. In contrast to this view, the current diversity of the Mediterranean area is increasingly seen to result from processes which are not directly related to the range movements of northern species during glacial-interglacial cycles [23]–[26], but our understanding of its origin remains fragmentary, particularly for highly speciose groups such as most insects.

In this study we integrate current and palaeoclimatic information with a molecular phylogeography, morphological analysis and experimentally derived thermal tolerance data to understand the role of thermal niche differences in shaping geographical expansion and speciation processes within the A. brunneus complex. Our specific goals are to: 1) test for climatic niche divergence among these species, and associate these differences with their current distribution; 2) test for differences in the estimated ecological niche of each species, and reconstruct the changes in their potential distributions through the last glacial cycle; and 3) evaluate species limits and reconstruct the speciation processes, demographic evolution and range expansion within the A. brunneus complex.

We use mitochondrial and nuclear sequence data from populations throughout the extant geographical ranges of all three species of the complex to reconstruct their demographic history and geographic expansion, and explore these within the context of changes in estimated potential distributions through the last glacial cycle. Using the current distribution of the species, we test for differences in climatic niche, and contrast these with experimental data obtained from range edge populations. By integrating these diverse data we are able to reconstruct the evolutionary history of the A. brunneus complex in the southwestern Mediterranean region, and illustrate how late Pleistocene climate changes may have shaped its current diversity by promoting ecological differentiation within a southern refuge.

Methods

Taxonomic background on the Agabus brunneus complex

The Agabus brunneus complex (Coleoptera, Dytiscidae) includes three recognised species of diving beetles with a distribution centred in the western Mediterranean area [16]: Agabus brunneus (Fabricius, 1798), A. rufulus Fairmaire, 1860 and A. ramblae Millán & Ribera [16]. Together with the more distantly related A. didymus (Olivier, 1795), which is more widely distributed in the western Palaearctic, they form the Agabus brunneus group [27],[28]. Agabus brunneus has a wide distribution through North Africa and western Europe, including the Iberian and Italian peninsulas, some Mediterranean islands (Elba, Sicily) and France, with isolated populations in southern England ([16]; Figure 1). Some old, isolated records in the eastern Mediterranean (Greece, Syria) with all probability refer to other species (e.g. A. dilatatus (Brullé, 1832), unpublished observations). Agabus ramblae was recognised based on external morphology and male genitalia as a distinct species previously confounded with A. brunneus, and has a disjunct distribution in the South and East of the Iberian Peninsula, the Balearic Islands, Central Morocco and Tunisia [29]. It is usually found in mineralized, temporary running waters [16], and despite its more recent description, exhaustive re-examination of material of A. brunneus sensu lato suggests that its apparently restricted range is genuine. The third species, A. rufulus, was traditionally considered a colour variety of A. brunneus and recorded from various localities in Italy (including the Tyrrhenian islands), Spain and North Africa [30],[31]. The revision of Millán & Ribera [16] revealed that it is, instead, a Corsico-Sardinian endemic. Agabus ramblae and A. rufulus have fully allopatric distributions, but the range of A. brunneus completely overlaps with that of A. ramblae, and it has also been recorded in Corsica and Sardinia. Agabus brunneus and A. ramblae are very rarely syntopic, with only 10 reported co-occurrences in the same locality (9 in the province of Albacete and 1 in the nearby area of Jaén, A. Millán, personal communication, 2014). Prior to this study very limited mitochondrial data were available for the three species [32],[33], and details of their phylogenetic/phylogeographic relationships, or age of divergence, were lacking.

Figure 1
figure 1

Known distribution of the species of the Agabus brunneus complex, including all localities used for DNA sequencing and ecological niche modelling. Blue circles, A. brunneus; red circles, A. ramblae; green circles A. rufulus.

Full size image

Morphological identification of species

The main diagnostic difference among species of the A. brunneus complex is the shape of the median lobe of the male aedeagus, in particular its relative size, the degree of asymmetry in ventral view and the shape of the apex in lateral view [16], and it is this character suite which was primarily used to assign material to species here. Female specimens were assigned to species by association with males and according to body dimensions and colouration. For 121 males of the three species we measured the maximum length of the median lobe of the aedeagus (AL) and the asymmetry of the aedeagus in ventral view (AD), as these are the main characters separating A. brunneus from A. ramblae, the two coexisting species [16]. We measured asymmetry (AD) as the difference between the width of the right (RD) and left (LD) sides in the point of its maximum width. We also measured the total body length from the anterior margin of the pronotum to the apex of the elytra (BL), as well as the maximum body width at the widest point (BW) (Additional file 1: Figure S1; Additional file 2: Table S1). Measured males included all specimens used for DNA sequencing (see below) and additional material from a number of sources, including areas for which no fresh material could be obtained (such as Mallorca). Species of the genus Agabus have a very uniform shape, and the length of a specimen is a good surrogate of total size [34]. For both A. brunneus and A. ramblae there is no difference between the body length of males and females [16]. We tested for differences in the length and shape of the aedeagus across the three species using MANOVA, and checked for possible intraspecific geographical variation within A. ramblae from Morocco, the Iberian Peninsula and the Balearic Islands, and A. rufulus from Corsica and Sardinia. All analyses were performed in IBM SPSS Statistics v. 20 (IBM, Armonk, NY, US).

Thermal tolerance

We determined the thermal tolerance of two populations of different species (Additional file 3: Table S5), one of A. brunneus from NE Spain (Girona, Ser river 42°08'48''N 2°34'48''E) and one A. ramblae from central Morocco (Tinghir, Toudgha river 31°33'25''N 5°34'49''W). The experimental procedure followed that used by Calosi et al. [17], from which we extracted data for three additional populations of three different species (A. brunneus, A. ramblae and A. rufulus). Individuals were acclimated for 7 days at 14.5°C in controlled conditions of pH, water composition, light regime and food (red chironomid larvae). Specimens were then separated into sub-groups and thermally ramped (±1°C min–1) in a computer-controlled water bath (Grant Instruments Ltd, Herts, UK) to obtain measures of their Upper Thermal Limit (UTL) and Lower Thermal Limit (LTL). Temperature was directly measured in one of the wells where individuals were placed for the experiment with a calibrated digital thermometer (Omega HH11; Omega Engineering Inc., Stamford, CT, USA) (see [17] for details). Data were analysed with an ANOVA with species and population as factors and DHS or Tukey post-hoc tests using IBM SPSS Statistics v. 20.

Taxon sampling, DNA extraction and sequencing

Specimens were collected in the field and directly preserved and stored in absolute ethanol. We included molecular data from 68 populations of A. brunneus, 22 A. ramblae and 6 A. rufulus, with up to five individuals per location when available, giving a total of 203 sequenced individuals covering the entire geographical ranges of the three species (Figure 1; Additional file 4: Table S2). As outgroups we used A. didymus (the sister of the A. brunneus complex [16],[27]) together with other published sequences representing a wide range of genera/species of Agabini (Additional file 4: Table S2). Trees were rooted in the genus Platynectes, which is within the Agabini but clearly outside the Agabus group of genera [32],[35],[36].

For DNA isolation we employed commercial DNA tissue kits (Additional file 4: Table S2) following the manufacturer instructions. Voucher specimens and DNA aliquots are deposited in the Natural History Museum (NHM, London), Museo Nacional de Ciencias Naturales (MNCN, Madrid) and Institute of Evolutionary Biology (IBE, Barcelona) (Additional file 4: Table S2).

To define the closest outgroups and the general time frame of diversification we used a combination of mitochondrial (a fragment of 827 nucleotides at the 3’ end of cox1, a continuous fragment between 798–803 (Agabini) or 802 nucleotides (A. brunneus complex) including the 3’ end of rrnL + full trnL + 5’ end of nad1) and two fragments of the nuclear genes SSU and H3 of 608 and 327 nucleotides respectively (see Additional file 5: Table S3 for the primers used and the general PCR conditions). For some specimens in which the cox1 fragment could not be amplified in a single PCR we used internal primers to obtain two non-overlapping fragments of ca. 400 bp each.

For the detailed phylogeographic and coalescent analyses we sequenced two gene fragments, one mitochondrial (3’-cox1) and one nuclear (H3). Sequence errors/ambiguities were edited using Geneious Pro 5.3.6 (http://www.geneious.com). New sequences have been deposited in GeneBank with accession numbers LM654767-LM655064 and LM655068-LM655168 (Additional file 4: Table S2).

Phylogenetic analyses

Length-variable sequences were aligned with the on-line version of MAFFT v.7 [37] using the Q-INS-i algorithm, which considers the secondary structure of RNA, and default values for other parameters. The final aligned matrix for the analyses of Agabini was 2579 nucleotide long.

General phylogenetic relationships

To determine the relationships among the main lineages within the A. brunneus complex and its phylogenetic relationships within Agabini we used the combined mitochondrial and nuclear sequence from a selection of specimens of the three species. Specimens were selected to cover the geographic range of all species, with a particular focus on potential contact areas (identified through preliminary analyses). Analyses used Bayesian probabilities as implemented in BEAST v1.7 [38] with a partition by genes (except for the trnL, pooled with the rrnL fragment) and a GTR + I + G evolutionary model for each partition. BEAST was run for 100 million generations, with 10% considered as the burn-in fraction after checking convergence of all parameters with the effective sample size (ESS) as measured in TRACER v1.5 [38].

To establish a temporal framework for the origin and evolution of the A. brunneus complex we used the mitochondrial genes only, for which there are recent calibrations for different families of Coleoptera with very homogeneous estimations for the rate of a combination of protein coding and ribosomal mitochondrial regions, calibrated with fossils and different biogeographic events [39]–[41]. As a prior we used a normal distribution with an average combined rate of 0.01 substitutions/site/million years (MY) and a standard deviation of 0.001, with other settings identical to the above analysis. To ensure that we obtained the same topology as in the analysis employing the full sequence, we constrained the monophyly of ingroup and outgroup, the genera and the A. brunneus complex. We used TRACER to calculate the mean and 95% highest posterior density interval for divergence times. We tested different alternative topologies for the relationships amongst species of the A. brunneus complex via the use of Bayes factors as estimated with the stepping-stone (SS) and the path-sampling (PS) algorithms in BEAST [42], and with the harmonic mean estimator (HME) in TRACER 1.5 [43] for comparison, in this case requiring an improvement in marginal likelihood of 10 units per additional parameter before accepting a more complex model [44],[45]. We tested three different topologies: 1) unconstrained (C0); 2) respective monophyly of A. brunneus and A. ramblae (C1); and 3) monophyly of A. rufulus (C2). We also analysed the H3 sequences separately in BEAST and RAxML, using a range of outgroups (Additional file 4: Table S2), a single partition with a GRT + G evolutionary model and the previously estimated age of the genus Agabus as prior to calibrate the tree in BEAST.

To test for alternative demographic models and to establish haplotype distribution and relationships we used the mitochondrial gene cox1 of all sequenced specimens of the A. brunneus complex, with the same settings as described for the analyses above, but with a mean rate of 0.02 (following [39],[40]), with a standard deviation of 0.001. Demographic models were tested first with the haplotypes of the A. brunneus complex only (i.e. excluding A. didymus), without topological constraints, and then with the haplotypes of A. brunneus only (i.e. also excluding A. ramblae and A. rufulus) (Additional file 4: Table S2).

Four models were tested: 1) constant population size; 2) exponential growth; 3) expansion; and 4) logistic growth. Models were compared through Bayes factors as above, i.e. using the HME, PS and SS. We also constructed a Bayesian skyline plot (BSP, [46]) with the combined results of two independent analyses of the A. brunneus complex in TRACER v1.5. The BSP constructs a model of demographic history based on how the number of coalescent events over a given interval differs from that expected under a neutral model for a panmictic population, then summarizes all possible genealogies and provides confidence intervals for all parameters in the model. It estimates changes in effective population size to analyse the population expansion of a species.

Population genetic analysis

We estimated some measures of haplotype diversity and analysed raggedness indices for demographic expansion with Arlequin 3.5 [47] using the cox1 gene. We tested the validity of the estimated stepwise expansion model using a parametric bootstrap approach.

Ecological niche modelling

We used ecological niche modelling (ENM) based on large-scale climatic data sets and known occurrence points to characterise the environmental niche of all three species, and to test for niche divergence amongst them. Climatic data were obtained at a spatial resolution of approximately 0.08° from WORLDCLIM version 1.3 (19 bioclimatic variables from http://www.worldclim.org; [48]; Additional file 6: Table S4). As records of species occurrences, we employed all known localities for the A. brunneus species complex identified according to the morphology of the males (Additional file 7: Table S6) at the same resolution than the bioclimatic variables. Most of the localities used in the ENM were also represented in the molecular analyses.

Bioclimatic values were first subjected to a principal component analysis (PCA) to obtain uncorrelated environmental factors (Varimax rotation). We used the values of the first two PCA factors to represent the climatic space of the whole study area. We plotted the occurrences of the three species in this same space, to visualize the section of the climatic space occupied by each species.

As the results of the thermal tolerance experiments clearly pointed to the importance of lower thermal limits, we used a Kruskal-Wallis ANOVA to compare the minimum temperatures of the coldest month between the localities of the three species. Multiple comparison tests were used to detect significant differences between means. The PCA and the Kruskal-Wallis ANOVA were conducted in Statistica version 8.0 (www.statsoft.com, 2007).

To compare the climatic niche of species we generated ecological niche models using MaxEnt [49], with Schoener’s D [50] as a measure of niche similarity between each pair of species as calculated by ENMTOOLS [51]. These values were calculated by comparing the climatic suitability of each grid cell in the study area obtained with MaxEnt. As niche differences may be simply a result of the spatial autocorrelation of the explanatory environmental variables [52], we conducted a background similarity test, also implemented in ENMTOOLS. This test uses randomization to determine whether two species are more or less similar than expected based on the differences in the environmental background in which they occur. A null distribution of 100 niche similarity values was generated by comparing the model suitability values of one species to those generated from random cells drawn from the distribution of the other species. The observed D value of niche similarity between the two species was then compared with the null distribution generated for each of them. The background area of each species should be adjusted to the habitat available, and should be biologically realistic [51]. For the insular A. rufulus, the background area was geographically restricted to Corsica and Sardinia, and for A. ramblae and A. brunneus we defined the background area as the Freshwater Ecoregions in which each species occurs following the classification in Abell et al. [53] (Additional file 8: Figure S2). This method has been used in a number of studies (e.g. [54],[55]), including aquatic Coleoptera [56],[57].

To provide a climatic context for the interpretation of the demographic models and the changes in distribution of species, we estimated the potential distribution of A. brunneus and A. ramblae for current climatic conditions, the reconstructed conditions during the last glacial maximum (LGM, 21,000 years before present, YBP) and the last interglacial (LIG, ca. 120,000-140,000 YBP). To estimate potential (not realized) distributions we used a multidimensional-envelope, as it provides a better estimate from observed occurrences (see [58] for details, Additional file 9: Text S1). For both past scenarios we used the same 19 bioclimatic variables at the same resolution as for current climate (see above). For the LGM we used a simulation of the general circulation model (GCM) from the Community Climate System Model (CCSM, http://www.ccsm.ucar.edu/, [59]). The original GCM data were downloaded from the PMIP2 website (http://www.pmip2.cnrs-gif.fr/). For the LIG we used data provided by Otto-Bliesner et al. [60], available at www.worldclim.org.

Results

Morphology

The three species differed in body size, as measured with BL and BW (MANOVA, Roy’s greatest root F = 292.134, P < 0.001), with pairwise differences being significant for all comparisons (Table 1). Agabus brunneus was the largest of the three species, and A. ramblae the smallest. The measures of size and asymmetry of the aedeagus (AL and AD respectively) were also significantly different (MANOVA, Roy’s greatest root F = 513.120, P < 0.001, Table 1). Agabus ramblae and A. brunneus were fully separated by AD, with no intermediate specimens, whilst the specimens identified as A. rufulus had an intermediate, overlapping shape (Additional file 10: Figure S3). For both A. ramblae (comparing Morocco, Iberian peninsula and the Balearic Islands) and A. rufulus (Corsica and Sardinia) there were no significant differences between major geographical areas.

Table 1 Comparison of the measurements of body and male genitalia (aedeagus) among the three species of the A. brunneus complex
Full size table

Thermal tolerance

Overall differences in thermal tolerance were highly significant between populations for both LTL (F = 9.87, d.f. = 4, P < 0.001) and UTL (F = 27.2, d.f. = 4, P < 0.001; Figure 2, Table 2, Additional file 3: Table S5). For UTL the highest differences were between the two populations of A. brunneus (north and south Iberian peninsula, post-hoc Tukey P < 0.001; higher limit in the southern population), which were also different for LTL (post-hoc Tukey P < 0.01; lower limit in the northern population). Differences between the two populations of A. ramblae (Morocco and SE Spain) were significant for UTL (post-hoc Tukey P < 0.05; higher limit in the northern population) but not for LTL (Figure 2; Table 2).

Figure 2
figure 2

Results of thermal tolerance experiments on the five populations tested. A) upper thermal limit (UTL), B) lower thermal limit (LTL). CAD, Cádiz; GIR, Girona; MUR, Murcia; MOR, Morocco; SAR, Sardinia.

Full size image
Table 2 Differences in thermal tolerance among the species and populations of the A. brunneus complex
Full size table

Between species there were significant differences in LTL (F = 10.9, d.f. = 2, P < 0.001; Figure 2; Table 2; lower limit in A. brunneus, higher in A. ramblae), but not UTL. Post-hoc analyses for LTL were highly significant in the case of A. brunneus vs. A. ramblae (P < 0.001) and A. brunneus vs. A. rufulus (P < 0.05), but not for A. ramblae vs A. rufulus (Table 2).

Phylogeny and phylogeography

Phylogenetic placement of the A. brunneus complex

There were no length differences in the mitochondrial and nuclear protein coding genes (nad1, cox1 and H3), and amongst the ribosomal genes length differences were restricted to outgroups, with a maximum difference of 2–5 bp. There were few amino-acid changes in the protein coding genes within the A. brunneus complex, and only two in the H3 fragment: one shared by all A. rufulus with the exception of one specimen from Sardinia (AH223), and another shared by all A rufulus and all A. brunneus.

In the Bayesian analysis of the combined data (mitochondrial plus nuclear) the genus Agabus was recovered as monophyletic with strong support (Figure 3), and the A. brunneus group also monophyletic but with lower support (posterior probability, pp = 0.61) and with a long stem branch. The monophyly of the A. brunneus complex was strongly supported (pp = 1), with A. ramblae basal and paraphyletic and a polyphyletic A. rufulus: two Corsican specimens were sister to Menorcan A. ramblae, while two Sardinian specimens were nested within A. brunneus.

Figure 3
figure 3

Phylogenetic relationships of the Agabus brunneus complex, obtained with BEAST and the combined mitochondrial and nuclear sequence. Numbers in black, node posterior probabilities; numbers in red, estimated age of the nodes (Ma) and 95% confidence intervals, obtained from the analysis of the mitochondrial data only (see Figure 4 and Additional file 11: Figure S4).

Full size image

Internal phylogeny of the A. brunneus complex, divergence age estimation

The estimation of a time window for the diversification of the A. brunneus complex, using the mitochondrial sequence from the above specimens and an a priori rate obtained from related groups, gave an age for the stem group of approximately 10.4 Ma (million years ago). The estimated origin of the sampled diversity within the complex was much more recent (0.6 Ma, Figure 3). Models using the two constraints tested (monophyly of A. ramblae, C1, and A. rufulus, C2) preformed significantly worse than the unconstrained model (C0), with the monophly of A. rufulus being the worst of all for all three Bayes factor estimators (HME, PS and SS; with more than 10 units in the difference in –lnLH for the HME, and more than four units in the case of PS and SS, [61], Table 3). The preferred mitochondrial topology had the Moroccan haplotypes of A. ramblae sister to the rest of the complex (range of uncorrected mitochondrial p distances, p = 0.004-0.012), Corsican A. rufulus sister to Menorcan A. ramblae with a divergence of ca. 0.34 Ma (p = 0.008), and Iberian A. ramblae sister to A. brunneus + Sardinian A. rufulus, with a divergence date of ca. 0.35 Ma (p = 0.002; Figure 4, Additional file 11: Figure S4).

Table 3 Bayes factors for the topological comparisons
Full size table
Figure 4
figure 4

Summary trees obtained from A) combined mitochondrial data ( cox1 + rrnl + trnl + nad1 ), B) cox1 and C) H3 . Node supports in A) and B), Bayesian posterior probabilities; C) Bayesian pp/ML bootstrap values. In red, estimated age of the nodes (see Additional file 11, Additional file 12 and Additional file 13 for the complete topologies of the trees).

Full size image

We also analysed the cox1 and H3 sequences independently. The analysis of the cox1 gene for all sequenced specimens (203, Additional file 4: Table S2), using A. didymus as an outgroup, resulted in a topology very similar to that described above, but with some additional Sardinian A. rufulus grouped with Corsican specimens and additional Iberian A. ramblae grouped with the Moroccan specimens and nested within A. brunneus (Figure 4, Additional file 12: Figure S5), in some cases with identical haplotypes. Two specimens of A. ramblae from Menorca (AH348, AH352, the later a male with typical A. ramblae aedeagus, Additional file 2: Table S1) and one from Morocco (AH311, female) were also nested within A. brunneus.

The analysis of the H3 gene recovered the three species as respectively monophyletic with good support (A. ramblae and A. rufulus pp = 1, ML bootstrap = 89 and 91 respectively; A. brunneus pp = 0.8, MLb = 74; Figure 4, Additional file 13: Figure S6), with only one exception, one female from Albacete identified as A. ramblae (AH224) was grouped with A. brunneus. Within A. brunneus and A. rufulus there was some variation (one position), without geographical structure (Figure 4, Additional file 13: Figure S6).

Demographic models of the expansion of the A. brunneus complex

We compared four coalescent demographic models in BEAST using only the cox1 data of the three species within the A. brunneus complex. The best model according to Bayes factors using the PS and SS estimators was logistic growth (with a difference of more than four units of –lnLH), followed by exponential growth (Table 3). On the contrary, marginal likelihood (HME) suggested a constant population size model as optimal (with a difference of more than 20 units of –lnLH), followed by logistic growth. The constant population model performed worst for both PS and SS estimators (Table 3).

The combined two independent runs of the Bayesian skyline plot gave a good convergence when the burn-in fraction was extended to 40 million generations (40%), and suggested exponential population growth at ca. -0.03 Ma, followed by a levelling off at ca. -0.01 Ma with a slight decrease towards the present, i.e. a sigmoidal curve of population growth (Figure 5).

Figure 5
figure 5

Bayesian skyline plot of the cox1 sequence data of the 203 studied specimens of the Agabus brunneus complex. Blue lines, 95% highest probability density; horizontal axis, time before present (Ma); vertical axis, effective population size (NeT).

Full size image

The expansion raggedness indices were lower for populations of A. brunneus and A. rufulus (from Corsica) than for A. ramblae, indicating an expansion in the former two species (Table 4). Agabus ramblae also had a higher molecular diversity, as measured with the Theta S and Theta Pi indexes (Table 4).

Table 4 Measures of raggedness and molecular diversity of the three species of the A. brunneus complex
Full size table

Ecological niche modelling data

The two first axes of the PCA of climatic variables for all localities of the A. brunneus complex jointly accounted for 82.4% of the total variance, and were interpreted as representing `aridity’ and `seasonality’ gradients respectively. The first axis was positively correlated to maximum temperature of the warmest month, and negatively to precipitation of the driest month, whilst the second axis was negatively correlated with temperature seasonality. The environmental space of A. brunneus encompassed almost completely that of the other two species: Agabus ramblae occupied the more seasonal and arid extreme of the climatic space of A. brunneus complex, and A. rufulus was climatically close to A. ramblae, although in areas with lower aridity and seasonality (Additional file 14: Figure S7).

The values of the minimum temperature of the coldest month of the three species (A. brunneus: –8.8°C, A. ramblae: –4.8°C, A. rufulus: –1.9°C) were significantly different, as measured with a Kruskal-Wallis test (N = 686; H = 73.32, P < 0.05). All pairwise comparisons were also significantly different (at P < 0.05) except that for A. brunneus and A. ramblae, which was close to significance (P = 0.07).

Results of ENM of each species, estimated with MaxEnt, had low spatial overlap. The degree of niche overlap estimated by Schoener’s D statistic was lower than 0.373 for all pairwise comparisons between the three species, suggesting differences in the climatic niche among them ([52]; Figure 6 and Figure 7). The null hypothesis of no differences in ecological niches explained by environmental differences between areas of occupancy alone was accepted (lower p-value = 0.17) for the comparison between A. ramblae and A. rufulus (despite a low D value, 0.20, this was not significantly more different than expected by chance, Figure 7), and rejected (p < 0.05) for comparisons between A. brunneus and the other two species.

Figure 6
figure 6

Environmental niche models (ENMs) using MaxEnt for A) A. brunneus ; B) A. ramblae ; and C) A. rufulus . Higher MaxEnt probabilities (red colours) represent areas more suitable for the species according to the MaxEnt models, lower values (blue) represent less suitable areas. These continuous measures of habitat suitability were used in ENM-based tests of niche similarity.

Full size image
Figure 7
figure 7

Background similarity test results for species within the Agabus brunneus species complex. Observed niche overlap values (arrows) were compared with null distributions (100 replicates) generated by comparing model suitability values of one species to those generated from random cells drawn from the background area of the other species.

Full size image

The potential distributions of A. brunneus and A. ramblae during the last interglacial (LIG), reconstructed using a climatic envelope, were similar but far more restricted than their current distributions, with most of the Iberian peninsula and North Africa considered unsuitable (Figure 8). This was especially the case for A. ramblae, for which, in terms of current range, only the Balearic Islands and some areas in the High Atlas in Morocco appeared appropriate (Figure 8). The largest difference between the LIG reconstruction and the current climate for the studied variables was seasonality, with the areas considered to be unsuitable having values beyond those of their current ranges (Additional file 15: Figure S8).

Figure 8
figure 8

Estimated potential distribution for A. brunneus and A. ramblae during the last interglacial (LIG), last Glacial Maximum (LGM) and the present. In red, areas considered to be climatically suitable for the species (according to the ecological conditions of the localities in which they are currently found). Blue dots, current localities of the species.

Full size image

During the last glacial maximum (LGM) the potential distribution of the two species increased to cover most of their current ranges: for A. ramblae only the northernmost known locality was outside the potential LGM distribution (Alcampell, in the province of Huesca), which was also an outlier in the representation of the climatic niche of the species (Additional file 14: Figure S7). For A. brunneus some Pyrenean localities, those on the north coast of France and the south coast of Britain, north Italy and the French Massif Central were outside the LGM potential distribution (Figure 8), but most of its current distribution corresponded to its potential range in the LGM. On the contrary, for both species the estimated potential present day distribution was much wider than the actual range: for A. brunneus it included the whole Mediterranean area and most of Europe, and for A. ramblae most of the Mediterranean and central France (Figure 8).

Discussion

Species limits within the A. brunneus complex – morphology, molecules and physiology

Our initial criterion for species recognition was the morphology of the aedeagus, in agreement with current taxonomy [16]. The simple measures used were able to unambiguously discriminate between A. brunneus and A. ramblae, but A. rufulus, the insular species, had an intermediate morphology for these characters, although the shape of the aedeagus in lateral view allows the unequivocal identification of this species [16]. The three species were clearly recovered as monophyletic with the nuclear marker (H3), with the exception of one female from the southeast of the Iberian Peninsula, characterised as A. ramblae but nested with other peninsular A. brunneus, which may represent a misclassified individual (in the same area both A. brunneus and A. ramblae can be found, Additional file 4: Table S2). In contrast to nuclear genes, the mitochondrial markers recovered a paraphyletic A. ramblae as ancestral to A. brunneus and A. rufulus. This paraphyly could be due to incomplete lineage sorting resulting from the recent evolution of the group [62]–[64], which is estimated to have diverged mostly within the last 0.5 MY –an insufficient time to reach reciprocal monophyly. However, in our uncalibrated tree for Agabini the estimated rate of cox1 was approximately eight times higher than that of H3, which being nuclear, should have a longer coalescent time than mitochondrial genes [65]. A possible explanation could be a complete replacement of the A. ramblae mitochondrial genome by that of A. brunneus in the Iberian peninsula, due to an early introgression event which did not result in phenotypic change between North African and Iberian populations of A. ramblae (either morphological, ecological or in thermal tolerance). The clear morphological separation between the two species argues against continued events of gene flow between them, which if present should have produced a higher frequency of specimens with intermediate morphologies. But given the low variability among the species of the A. brunneus complex in the H3 fragment, this difference in coalescent times could also be simply due to random effects.

There is some additional evidence of a mismatch between morphology and some of the genetic markers used that suggests occasional introgression, but this is only seen in geographically marginal areas: Sardinia (between A. rufulus and A. brunneus), Menorca (between A. ramblae and A. brunneus), and possibly Tunisia (also between A ramblae and A. brunneus). Mitochondrial haplotypes of A. brunneus were found in individuals identified by morphology and nDNA as A. ramblae in Menorca, and A. rufulus in Sardinia. This could be due to a secondary colonisation of the islands by continental A. brunneus, similarly to what has been described in a related genus of diving beetle (Meladema) in the Canary Islands [66]. In Tunisia, southern populations were characterized as A. ramblae by both morphology and nDNA, but the mtDNA clustered with A. brunneus and Iberian haplotypes of A. ramblae. Two possible scenarios may account for such results: Tunisian A. ramblae could have arrived directly from the Iberian peninsula, which seems unlikely given the geographical distance and the presence of sea barriers; or they could have arrived from elsewhere in North Africa, and then hybridised with northern A. brunneus either also from North Africa or from Sicily. A Moroccan origin is supported by the presence of other water beetles typical of arid or saline habitats with a similar distribution through central and south Morocco to south Tunisia, such as Enochrus risii Arribas et al. [67] or Ochthebius salinator Peyerimhoff [68]. The situation of A. ramblae in the Iberian peninsula is more complex and difficult to interpret: some specimens had mitochondrial haplotypes clustering with those of A. ramblae from Morocco, suggesting the persistence within Iberia of some of the ancestral Moroccan haplotypes and a derived origin of most of both the Iberian A. ramblae and A. brunneus. But again, the replacement of the A. ramblae mitochondrial genome by that of A. brunneus through introgression in secondary contact zones cannot be discarded.

The two species with broadly overlapping ranges, A. brunneus and A. ramblae, were also ecologically different as measured through ecological niche modelling and the background test. Our experimental results show that the greater resistance of A. brunneus to lower temperatures may have been a key feature to allow its range expansion during the LGM (see below). This difference was reflected in the significantly lower minimum temperature of the coldest month of the places in which A. brunneus is currently found. This species had also significant differences in thermotolerance between populations, with the northern one (with an average lowest temperature of the coldest month of 2.8 °C) being more resistant to cold than the southern population (average lowest temperature of the coldest month of 6.0 °C). With our data it is not possible to discriminate whether this difference results from local adaptation or phenotypic plasticity, although Calosi et al. [17] found that members of the group did not significantly adjust their LTL after a short period of acclimation in the laboratory. This suggests that A. brunneus populations may instead adapt to local temperature conditions, these evolutionary changes possibly facilitating range expansions.

Evolutionary history of the A. brunneus complex

We found strong support for the monophyly of the A. brunneus complex, and also recovered a monophyletic A. brunneus group, albeit with lower support. The long stem branch of the complex is atypical within the rest of Agabini, and as there are no known species worldwide that could be more closely related to it [28], it appears that the A. brunneus complex has a very isolated position amongst extant Agabini.

Although the support for the internal relationships of the A. brunneus complex was low, the selected mitochondrial topology recovered the Moroccan populations of A. ramblae as paraphyletic, suggesting an origin of the complex in western North Africa, with colonization of the Iberian peninsula ca. 0.5 Ma. As seen above, much of the likely recent introgression between the species of the group is restricted to a few populations in secondary contact zones, so we do not expect this to affect the mitochondrial phylogeny. Even if some paraphyly was due to early introgression within the Iberian peninsula this would not affect our biogeographic scenario. The colonization of the Balearic Islands and Corsica and Sardinia happened in a narrow temporal window, possibly from Iberia. Although we cannot discard a direct colonization from North Africa, this seems less plausible due to the longer geographical distances involved. According to our estimation, A. brunneus split from SE Iberian A. ramblae ca. 0.25 Ma. An alternative scenario is that A. brunneus originated in Morocco and colonized the Iberian Peninsula in parallel with A. ramblae, but as well as being less parsimonious this hypothesis seems less plausible since some Iberian A. ramblae have mitochondrial haplotypes clustering with those from Morocco.

The ecological and physiological differences between A. brunneus and A. ramblae may have originated during the speciation process or evolved later, with an initial separation only due to isolation. Either way, at some point A. brunneus acquired the capacity to resist colder temperatures.

The demographic analyses estimated a population expansion of the complex at the start of the last glaciation 30–40,000 YBP, in agreement with the extension of their potential distributions during the last glacial maximum (LGM, 21,000 YBP). For both widespread species, potential distributions during the LGM covered practically the totality of their current ranges, and were mostly determined by minimum temperatures and climatic seasonality. It is remarkable that only very few current known localities (mostly for A. brunneus) are outside the reconstructed potential range of both species during the LGM, despite a large increase in apparently suitable geographical areas both in Europe and north Africa. During the LGM sea levels could have been up to 200 m lower, probably extending the suitable surface and potentially favouring the expansion of the continental species of the group to areas now isolated by sea barriers, such as south Britain.

The current absence of A. brunneus from central and northern Europe cannot be attributed to the effect of anthropogenic habitat modification, as there are no historic or Quaternary fossil records of the species outside its current range [24]. The external appearance of species of the A. brunneus complex is very characteristic, so that even incomplete remains would be recognizable, and in central and northern Europe the fossil record from the LGM and the Holocene is very complete [24]. Possible explanations for the absence of range expansion in A. brunneus and A. ramblae after the LGM are the presence of undetected climatic, biotic or ecological limiting factors, or simply a lack of sufficient time for these species to arrive at equilibrium with their potential ranges. All species of the A. brunneus complex are exclusively found in running waters, and such lotic taxa have, in general, weaker dispersal abilities than their lentic relatives, leading to a stronger mismatch between their realized and potential distributions in central and north Europe [66],[69].

Conclusions

Using a combination of morphology, genetics, ecological niche modelling of current and paleoecological data and physiological experiments we have reconstructed the surprisingly complex evolutionary history of this diving beetle clade in the western Mediterranean. The A. brunneus complex diversified ca. 0.6-0.25 Ma, most likely in the south of the Iberian peninsula after the colonization of A. ramblae from north Morocco. Whilst insular populations (A. ramblae in the Balearic Islands and A. rufulus in Corsica and Sardinia) did not apparently differentiate substantially in either morphology or ecology, continental A. brunneus evolved the most distinctive morphology within the complex, as well as wider tolerance to cold habitats, something that seems to have facilitated range expansion.

From a reduced potential distribution during the LIG, A. brunneus and A. ramblae appear to have expanded their ranges during the last glacial (0.03-0.01 Ma) (A. brunneus to a much wider area), covering most of their LGM potential rages in the western Mediterranean. This expansion was accompanied by a population expansion, as identified through demographic models. However, despite much wider current potential distributions, both species have not occupied areas beyond their LGM potential distribution except for some isolated populations of A. brunneus in France and England. In Sardinia, the Balearic Islands and possibly Tunisia, secondary contact between species of the complex has resulted in introgression, with some specimens showing discordance between mitochondrial haplotypes typical of A. brunneus and nuclear sequences and morphology typical of A. rufulus or A. ramblae respectively.

Our work highlights the complex dynamics of speciation and range expansions within refugia during the last glacial cycle, and the fact that the biota of southern Europe, in addition to being a source of colonisers of formerly glaciated areas in the north, experienced much evolutionary change during this time period. It also highlights the fundamental but often neglected role of North Africa as source of biodiversity in Europe [70]–[74].

Availability of supporting data

All raw data are included in the Supplementary files with the exception of the sequences, deposited in the EMBL database with accession numbers LM654767-LM655064 and LM655068-LM655168.

Authors’ contributions

AHG, AC and IR conceived the study. AHG and IR coordinated the sampling. AHG obtained most the sequences and the morphometric and distribution data. AHG and DTB conducted the physiological experiments. AHG, DSF and IR analysed the data. All authors contributed to the writing and improving the manuscript, and approved the final version.

Additional files

References

  1. Spicer JI, Gaston KJ: Physiological Diversity and its Ecological Implications. 1999, Blackwell Science, Oxford

    Google Scholar 

  2. Portner HO, Bennett AF, Bozinovic F, Clarke A, Lardies MA, Lucassen M, Pelster B, Schiemer F, Stillman JH: Trade-offs in thermal adaptation: the need for a molecular to ecological integration. Physiol Biochem Zool. 2006, 79: 295-313. 10.1086/499986.

    PubMed  Article  Google Scholar 

  3. Angilletta MJ, Niewiarowski PH, Navas CA: The evolution of thermal physiology in ectotherms. J Therm Biol. 2002, 27: 249-268. 10.1016/S0306-4565(01)00094-8.

    Article  Google Scholar 

  4. Calosi P, Bilton DT, Spicer JI, Votier SC, Atfield A: What determines a species’ geographical range? Thermal biology and latitudinal range size relationships in European diving beetles (Coleoptera: Dytiscidae). J Anim Ecol. 2010, 79: 194-204. 10.1111/j.1365-2656.2009.01611.x.

    PubMed  Article  Google Scholar 

  5. Peterson AT, Soberón J, Pearson RG, Anderson RP, Martínez-Meyer E, Nakamura M, Araújo MB: Ecological Niches and Geographic Distributions: a Modelling Perspective. 2011, Princeton, Princeton University Press

    Google Scholar 

  6. Carnaval AC, Moritz C: Historical climate modelling predicts patterns of current biodiversity in the Brazilian Atlantic forest. J Biogeogr. 2008, 35: 1187-1201. 10.1111/j.1365-2699.2007.01870.x.

    Article  Google Scholar 

  7. Kozak KH, Wiens JJ: Niche conservatism drives elevational diversity patterns in Appalachian salamanders. Am Nat. 2010, 176: 40-54. 10.1086/653031.

    PubMed  Article  Google Scholar 

  8. Pepper M, Fujita MK, Moritz C, Keogh JS: Palaeoclimate change drove diversification among isolated mountain refugia in the Australian arid zone. Mol Ecol. 2011, 20: 1529-45. 10.1111/j.1365-294X.2011.05036.x.

    PubMed  Article  Google Scholar 

  9. Moritz C, Langham G, Kearney M, Krockenberger A, VanDerWal J, Williams S: Integrating phylogeography and physiology reveals divergence of thermal traits between central and peripheral lineages of tropical rainforest lizards. Philos Trans R Soc B Biol Sci. 2010, 367: 1680-1687. 10.1098/rstb.2012.0018.

    Article  Google Scholar 

  10. Sánchez-Fernández D, Aragón P, Bilton DT, Lobo JM: Assessing the congruence of termal niche estimations derived from distribution and physiological data. A test using diving beetles. Plos One. 2012, 7 (10): e48163-10.1371/journal.pone.0048163.

    PubMed  PubMed Central  Article  Google Scholar 

  11. Weber MG, Agrawal A: Phylogeny, ecology, and the coupling of comparative and experimental approaches. Trends Ecol Evol. 2012, 27: 394-403. 10.1016/j.tree.2012.04.010.

    PubMed  Article  Google Scholar 

  12. Nesse RM: Tinbergen’s four questions, organized: a response to Bateson and Laland. Trends Ecol Evol. 2013, 28: 681-2. 10.1016/j.tree.2013.10.008.

    PubMed  Article  Google Scholar 

  13. Eme D, Malard F, Colson-Proch C, Jean P, Calvignac S, Konecny-Dupré L, Hervant F, Douady CJ: Integrating phylogeography, physiology and habitat modelling to explore species range determinants. J Biogeogr. 2014, 41: 687-699. 10.1111/jbi.12237.

    Article  Google Scholar 

  14. Beale CM, Lennon JJ, Gimona A: Opening the climate envelope reveals no macroscale associations with climate in European birds. Proc Natl Acad Sci U S A. 2008, 105: 14908-14912. 10.1073/pnas.0803506105.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  15. Dawson TP, Jackson ST, House JI, Prentice IC, Mace GM: Beyond predictions: biodiversity conservation in a changing climate. Science. 2011, 332: 53-58. 10.1126/science.1200303.

    PubMed  CAS  Article  Google Scholar 

  16. Millán A, Ribera I: The Agabus (Gaurodytes) brunneus group, with description of a new species from the western Mediterranean (Coleoptera: Dytiscidae). Coleopt Bull. 2001, 55: 107-112. 10.1649/0010-065X(2001)055[0107:TAGBGW]2.0.CO;2.

    Article  Google Scholar 

  17. Calosi P, Bilton DT, Spicer JI, Atfield A: Thermal tolerance and geographical range size in the Agabus brunneus group of European diving beetles (Coleoptera: Dytiscidae). J Biogeogr. 2008, 35: 295-305.

    Google Scholar 

  18. Hewitt GM: Some genetic consequences of ice ages, and their role, in divergence and speciation. Biol J Linn Soc. 1996, 58: 247-276. 10.1111/j.1095-8312.1996.tb01434.x.

    Article  Google Scholar 

  19. Hewitt GM: Genetic consequences of climatic oscillations in the Quaternary. Philos Trans R Soc Lond B Biol Sci. 2004, 359: 183-195. 10.1098/rstb.2003.1388.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  20. Schmitt T: Molecular biogeography of Europe: Pleistocene cycles and postglacial trends. Front Zool. 2007, 4: 13-10.1186/1742-9994-4-13.

    Article  Google Scholar 

  21. Gómez A, Lunt DH: Refugia within refugia: patterns of phylogeographic concordance in the Iberian Peninsula. Phylogeography of Southern European Refugia. Edited by: Weiss S, Ferrand N. 2004, Kluwer, Dordrecht, The Netherlands, 155-188.

    Google Scholar 

  22. Husemann M, Schmitt T, Zachos FE, Ulrich W, Habel JC: Palaearctic biogeography revisited: evidence for the existence of a North African refugium for Western Palaearctic biota. J Biogeogr. 2014, 41: 81-94. 10.1111/jbi.12180.

    Article  Google Scholar 

  23. Bilton DT, Mirol PM, Mascheretti S, Fredga K, Zima J, Searle JB: Mediterranean Europe as an area of endemism for small mammals rather than a source for northwards postglacial colonization. Proc R Soc B. 1998, 265: 1219-1226. 10.1098/rspb.1998.0423.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  24. Abellán P, Benetti CJ, Angus RB, Ribera I: A review of Quaternary range shifts in European aquatic Coleoptera. Glob Ecol Biogeogr. 2011, 20: 87-100. 10.1111/j.1466-8238.2010.00572.x.

    Article  Google Scholar 

  25. Hortal J, Diniz-Filho JAF, Bini LM, Rodríguez MÁ, Baselga A, Nogués-Bravo D, Rangel TF, Hawkins BA, Lobo JM: Ice age climate, evolutionary constraints and diversity patterns of European dung beetles. Ecol Lett. 2011, 14: 741-748. 10.1111/j.1461-0248.2011.01634.x.

    PubMed  Article  Google Scholar 

  26. Ribera I, Castro A, Díaz JA, Garrido J, Izquierdo A, Jäch MA, Valladares LF: The geography of speciation in narrow-range endemics of the “Haenydra” lineage (Coleoptera, Hydraenidae, Hydraena). J Biogeogr. 2011, 38: 502-516. 10.1111/j.1365-2699.2010.02417.x.

    Article  Google Scholar 

  27. Nilsson AN: World Catalogue of Insects, Vol. 3: Dytiscidae (Coleoptera). 2001, Apollo Books, Stenstrup, Denmark

    Google Scholar 

  28. Nilsson AN: A World Catalogue of the Family Dytiscidae, or the Diving Beetles (Coleoptera, Adephaga). Version 1.I.2013. Available at ; last accessed 30 March 2014., [http://www2.emg.umu.se/projects/biginst/andersn/]

  29. Przewoźny M, Jaskuła R, Rewicz T: First African records ofAgabus(Gaurodytes)ramblaeMillán & Ribera, 2001 with an updated checklist ofAgabusin Maghreb (Coleoptera: Dytiscidae).Afr Ent 2014, in press.,

  30. Guignot F: Les Hydrocanthares de France. 1933, Miscelanea Zoologica, Toulouse

    Google Scholar 

  31. Franciscolo ME: Fauna d’Italia, Vol. XIV: Coleoptera Haliplidae, Hygrobiidae, Gyrinidae, Dytiscidae. 1979, Edizioni Calderini, Bologna

    Google Scholar 

  32. Ribera I, Nilsson AN, Vogler AP: Phylogeny and historical biogeography of Agabinae diving beetles (Coleoptera) inferred from mitochondrial DNA sequences. Mol Phylogenet Evol. 2004, 30: 545-562. 10.1016/S1055-7903(03)00224-0.

    PubMed  CAS  Article  Google Scholar 

  33. Abellán P, Sánchez-Fernández D, Picazo F, Millán A, Lobo JM, Ribera I: Preserving the evolutionary history of freshwater biota in Iberian National Parks. Biol Conserv. 2013, 162: 116-126. 10.1016/j.biocon.2013.04.001.

    Article  Google Scholar 

  34. Ribera I, Nilsson AN: Morphometric patterns among diving beetles (Coleoptera: Noteridae, Hygrobiidae, and Dytiscidae). Can J Zool. 1995, 73: 2343-2360. 10.1139/z95-275.

    Article  Google Scholar 

  35. Nilsson AN: A new view on the generic classification of the Agabus-group of genera of the Agabini, aimed at solving the problem with a paraphyletic Agabus (Coleoptera: Dytiscidae ). Koleopterol Rundsch. 2000, 70: 17-36.

    Google Scholar 

  36. Ribera I, Vogler AP, Balke M: Phylogeny and diversification of diving beetles (Coleoptera: Dytiscidae). Cladistics. 2008, 24: 563-590. 10.1111/j.1096-0031.2007.00192.x.

    Article  Google Scholar 

  37. Katoh K, Misawa K, Kuma K, Miyata T: MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002, 30: 3059-3066. 10.1093/nar/gkf436.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  38. Drummond AJ, Suchard MA, Xie D, Rambaut A: Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012, 29: 1969-1973. 10.1093/molbev/mss075.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  39. Papadopoulou A, Anastasiou I, Vogler AP: Revisiting the insect mitochondrial molecular clock: the mid-Aegean trench calibration. Mol Biol Evol. 2010, 27: 1659-1672. 10.1093/molbev/msq051.

    PubMed  CAS  Article  Google Scholar 

  40. Ribera I, Fresneda J, Bucur R, Izquierdo A, Vogler AP, Salgado JM, Cieslak A: Ancient origin of a Western Mediterranean radiation of subterranean beetles. BMC Evol Biol. 2010, 10: 29-10.1186/1471-2148-10-29.

    PubMed  PubMed Central  Article  Google Scholar 

  41. Andújar C, Serrano J, Gómez-Zurita J: Winding up the molecular clock in the genus Carabus (Coleoptera: Carabidae): assessment of methodological decisions on rate and node age estimation. BMC Evol Biol. 2012, 12: 40-10.1186/1471-2148-12-40.

    PubMed  PubMed Central  Article  Google Scholar 

  42. Baele G, Lemey P, Bedford T, Rambaut A, Suchard MA, Alekseyenko AV: Improving the accuracy of demographic and molecular clock model comparison while accommodating phylogenetic uncertainty. Mol Biol Evol. 2012, 29: 2157-67. 10.1093/molbev/mss084.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  43. Drummond AJ, Rambaut A: BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007, 7: 214-10.1186/1471-2148-7-214.

    PubMed  PubMed Central  Article  Google Scholar 

  44. Pagel M, Meade A, Barker D: Bayesian estimation of ancestral character states on phylogenies. Syst Biol. 2004, 53: 673-684. 10.1080/10635150490522232.

    PubMed  Article  Google Scholar 

  45. Miller KB, Bergsten J, Whiting MF: Phylogeny and classification of the tribe Hydaticini (Coleoptera: Dytiscidae): partition choice for Bayesian analysis with multiple nuclear and mitochondrial protein-coding genes. Zool Scr. 2009, 38: 591-615. 10.1111/j.1463-6409.2009.00393.x.

    Article  Google Scholar 

  46. Drummond AJ, Rambaut A, Shapiro B, Pybus OG: Bayesian coalescent inference of past population dynamics from molecular sequences. Mol Biol Evol. 2005, 22: 1185-92. 10.1093/molbev/msi103.

    PubMed  CAS  Article  Google Scholar 

  47. Excoffier L, Lischer HEL: Arlequin suite ver. 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010, 10: 564-567. 10.1111/j.1755-0998.2010.02847.x.

    PubMed  Article  Google Scholar 

  48. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A: Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005, 25: 1965-1978. 10.1002/joc.1276.

    Article  Google Scholar 

  49. Phillips SJ, Anderson RP, Schapire RE: Maximum entropy modelling of species geographic distributions. Ecol Modell. 2006, 190: 231-259. 10.1016/j.ecolmodel.2005.03.026.

    Article  Google Scholar 

  50. Schoener TW: The Anolis lizards of Bimini: resource partitioning in a complex fauna. Ecology. 1968, 49: 704-726. 10.2307/1935534.

    Article  Google Scholar 

  51. Warren DL, Glor RE, Turelli M: ENMTools: a toolbox for comparative studies of environmental niche models. Ecography. 2010, 33: 607-611. 10.1111/j.1600-0587.2009.06041.x.

    Article  Google Scholar 

  52. Warren DL, Glor RE, Turelli M: Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution. Evolution. 2008, 62: 2868-2883. 10.1111/j.1558-5646.2008.00482.x.

    PubMed  Article  Google Scholar 

  53. Abell R, Thieme ML, Revenga C, Bryer M, Kottelat M, Bogutskaya N, Coad B, Mandrak N, Balderas SC, Bussing W, Stiassny MLJ, Skelton P, Allen GR, Unmack P, Naseka A, Ng R, Sindorf N, Robertson J, Armijo E, Higgins JV, Heibel TJ, Wikramanayake E, López HL, Reis RE, Lundberg JG, Pérez MHS, Petry P: Freshwater ecoregions of the World: a new map of biogeographic units for freshwater biodiversity conservation. Bioscience. 2008, 58: 403-414. 10.1641/B580507.

    Article  Google Scholar 

  54. Martínez-Gordillo D, Rojas-Soto O, Espinosa de los Monteros A: Ecological niche modelling as an exploratory tool for identifying species limits: an example based on Mexican muroid rodents. J Evol Biol. 2010, 23: 259-270. 10.1111/j.1420-9101.2009.01897.x.

    PubMed  Article  Google Scholar 

  55. McCormack JE, Zellmer AJ, Knowles LL: Does niche divergence accompany allopatric divergence in Aphelocoma jays as predicted under ecological speciation? Insights from tests with niche models. Evolution. 2010, 64: 1231-1244.

    PubMed  Google Scholar 

  56. Hawlitschek O, Porch N, Hendrich L, Balke M: Ecological niche modelling and nDNA sequencing support a new, morphologically cryptic beetle species unveiled by DNA barcoding. PLoS One. 2011, 6: e16662-10.1371/journal.pone.0016662.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  57. Sánchez-Fernández D, Lobo JM, Abellán P, Millán A: Environmental niche divergence between genetically distant lineages of an endangered water beetle. Biol J Linn Soc. 2011, 103: 891-903. 10.1111/j.1095-8312.2011.01668.x.

    Article  Google Scholar 

  58. Jiménez-Valverde A, Lobo JM, Hortal J: Not as good as they seem: the importance of concepts in species distribution modelling. Divers Distrib. 2008, 14: 885-890. 10.1111/j.1472-4642.2008.00496.x.

    Article  Google Scholar 

  59. Kiehl JT, Gent PR: The Community Climate System Model, Version 2. J Clim. 2004, 17: 3666-3682. 10.1175/1520-0442(2004)017<3666:TCCSMV>2.0.CO;2.

    Article  Google Scholar 

  60. Otto-Bliesner BL, Brady EC, Clauzet G, Tomas R, Levis S, Kothavala Z: Last Glacial Maximum and Holocene climate in CCSM3. J Clim. 2006, 19: 2526-2544. 10.1175/JCLI3748.1.

    Article  Google Scholar 

  61. Kass RE, Raftery AE: Bayes factors. J Am Stat Assoc. 1995, 90: 773-795. 10.1080/01621459.1995.10476572.

    Article  Google Scholar 

  62. Funk DJ, Omland KE: Species-level paraphyly and polyphyly: frequency, causes, and consequences, with insights from animal mitochondrial DNA. Annu Rev Ecol Evol Syst. 2003, 34: 397-423. 10.1146/annurev.ecolsys.34.011802.132421.

    Article  Google Scholar 

  63. Kizirian D, Donnelly MA: The criterion of reciprocal monophyly and classification of nested diversity at the species level. Mol Phylogenet Evol. 2004, 32: 1072-1076. 10.1016/j.ympev.2004.05.001.

    PubMed  Article  Google Scholar 

  64. Hörandl E: Paraphyletic versus monophyletic taxa — evolutionary versus cladistic classifications. Taxon. 2006, 55: 564-570. 10.2307/25065631.

    Article  Google Scholar 

  65. Palumbi SR, Cipriano F, Hare MP: Predicting nuclear gene coalescence from mitochondrial data: the three-times rule. Evolution. 2001, 55: 859-868. 10.1554/0014-3820(2001)055[0859:PNGCFM]2.0.CO;2.

    PubMed  CAS  Article  Google Scholar 

  66. Ribera I, Foster GN, Vogler AP: Does habitat use explain large scale diversity patterns in European water beetles?. Ecography. 2003, 26: 145-152. 10.1034/j.1600-0587.2003.03271.x.

    Article  Google Scholar 

  67. Arribas P, Andújar C, Sánchez-Fernández D, Abellán P, Millán A: Integrative taxonomy and conservation of cryptic beetles in the Mediterranean region (Hydrophilidae). Zool Scr. 2013, 42: 182-200. 10.1111/zsc.12000.

    Article  Google Scholar 

  68. Abellán P, Millán A, Ribera I: Parallel habitat-driven differences in the phylogeographical structure of two independent lineages of Mediterranean saline water beetles. Mol Ecol. 2009, 18: 3885-3902. 10.1111/j.1365-294X.2009.04319.x.

    PubMed  Article  Google Scholar 

  69. Sánchez-Fernández D, Lobo JM, Millán A, Ribera I: Habitat persistence mediates time to equilibrium in the geographical distribution of Iberian diving beetles. Global Ecol Biogeogr. 2012, 21: 988-997. 10.1111/j.1466-8238.2011.00743.x.

    Article  Google Scholar 

  70. Cosson J-F, Hutterer R, Libois R, Sarà M, Taberlet P, Vogel P: Phylogeographical footprints of the Strait of Gibraltar and Quaternary climatic fluctuations in the western Mediterranean: a case study with the greater white-toothed shrew, Crocidura russula (Mammalia: Soricidae). Mol Ecol. 2005, 14: 1151-1162. 10.1111/j.1365-294X.2005.02476.x.

    PubMed  CAS  Article  Google Scholar 

  71. Fritz U, Barata M, Busack SD, Fritzsch G, Castilho R: Impact of mountain chains, sea straits and peripheral populations on genetic and taxonomic structure of a freshwater turtle, Mauremys leprosa (Reptilia, Testudines, Geoemydidae). Zool Scr. 2006, 35: 97-108. 10.1111/j.1463-6409.2005.00218.x.

    Article  Google Scholar 

  72. Horn A, Roux-Morabito G, Lieutier F, Kerdelhue C: Phylogeographic structure and past history of the circum-Mediterranean species Tomicus destruens Woll. (Coleoptera: Scolytinae). Mol Ecol. 2006, 15: 1603-1615. 10.1111/j.1365-294X.2006.02872.x.

    PubMed  CAS  Article  Google Scholar 

  73. Habel J, Meyer M, Mousadik A, Schmitt T: Africa goes Europe: The complete phylogeography of the marbled white butterfly species complex Melanargia galathea/M. lachesis (Lepidoptera: Satyridae). Org Divers Evol. 2008, 8: 121-129. 10.1016/j.ode.2007.04.002.

    Article  Google Scholar 

  74. Habel JC, Dieker P, Schmitt T: Biogeographical connections between the Maghreb and the Mediterranean peninsulas of southern Europe. Biol J Linn Soc. 2009, 98: 693-703. 10.1111/j.1095-8312.2009.01300.x.

    Article  Google Scholar 

Download references

Acknowledgements

We especially thank all the collectors mentioned in Additional file 4: Table S2 for sending material for study; J.A. Carbonell for the habitus photographs of A. brunneus and A. ramblae; A. Faille, H. Fery and P. Queney for unpublished information on the distribution of A. brunneus and A. rufulus; and A. Izquierdo, R. Alonso and A. Faille for help in the laboratory. We also thank A. Millán and P. Abellán for comments and unpublished information, and two anonymous Referees for comments and suggestions. Permits to collect A. brunneus and A. ramblae were obtained from the different Spanish Regional governments. No permits or ethical approval were required for the experimental procedures. This work was supported by an FPI grant to AH-G and projects CGL2007-61665 and CGL2010-15755 from the Spanish government to IR. We acknowledge support of the publication fee by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).

Author information

Affiliations

  1. Institute of Evolutionary Biology (CSIC-Universitat Pompeu Fabra), Passeig Maritim de la Barceloneta 37-49, Barcelona, 08003, Spain

    Amparo Hidalgo-Galiana, David Sánchez-Fernández, Alexandra Cieslak & Ignacio Ribera

  2. Marine Biology and Ecology Research Centre, Plymouth University, Drake Circus, Plymouth, PL4 8AA, UK

    David T Bilton

Authors
  1. Amparo Hidalgo-Galiana
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. David Sánchez-Fernández
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. David T Bilton
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Alexandra Cieslak
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Ignacio Ribera
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding authors

Correspondence to Amparo Hidalgo-Galiana or Ignacio Ribera.

Additional information

Competing interest

The authors declare that they have no competing interest.

Electronic supplementary material

12862_2014_187_MOESM1_ESM.pdf

Additional file 1: Figure S1.: Measures used for the identification of the specimens. A) Median lobe of the aedeagus of A. ramblae in ventral view, with the measurements used. The global measure of asymmetry used was AD = RD-LD. B) Maximum body length (BL, excluding head) and width (BW). (PDF 768 KB)

12862_2014_187_MOESM2_ESM.xlsx

Additional file 2: Table S1.: List of specimens with morphometric measurements. Ref., reference of the specimen: with AH and AI, voucher numbers of extracted specimens (see Additional file 4: Table S2); BM, other material not used for DNA extraction. BL, body length (from anterior side of prototum to apex of elytra); BW, maximum body width; AL, length of the aedeagus; AD, asymmetry of the aedeagus (see Additional file 1: Figure S1). (XLSX 41 KB)

12862_2014_187_MOESM3_ESM.xlsx

Additional file 3: Table S5.: Thermal tolerance data. Population: SPA-Cádiz: Spain, Cádiz (data from ref.[17]); SPA-Girona: Spain, Girona (42°08'48''N 2°34'48''E); SPA-Murcia: Spain, Murcia (data from ref.[17]); MOR-Tinghir: Morocco, Tinghir (31°33'25''N 5°34'49''W); ITA-Sardinia: Italy, Sardinia (data from ref.[17]). (XLSX 48 KB)

12862_2014_187_MOESM4_ESM.xlsx

Additional file 4: Table S2.: List of the specimens used for DNA extraction, with accession numbers of the sequences. Extraction method: Charge, Charge Switch gDNA Tissue Kit (Invitrogen, Carlsbad, USA); Invisorb, Invisorb Spin DNA Extraction Kit (Invitek GmbH, Berlin, Germany); Qiagen, Qiagen DNeasy Tissue Kit (Qiagen GmbH, Hilden Germany). In the column 16S, with asterisk: single fragment including the 3′ end of rrnL, the full trnL and the 5′ end of nad1; in the rest either the fragment is divided in two sequences or includes only the rrnL fragment. (XLSX 72 KB)

Additional file 5: Table S3.: List of primers used and PCR conditions. (XLSX 42 KB)

Additional file 6: Table S4.: Climatic variables used in the niche modelling. (XLSX 36 KB)

12862_2014_187_MOESM7_ESM.xlsx

Additional file 7: Table S6.: Localities used for the niche modelling of the three species. BD (BIODIV), unpublished database of the Department of Ecology and Hydrology, University of Murcia (Spain); CKmap, Checklist and distribution of Italian fauna (http://www.faunaitalia.it/ckmap/); ESACIB, database “Escarabajos Acuáticos Ibéricos” (Sánchez-Fernández et al. 2008); GBIF, Global Biodiversity Information Facility (http://www.gbif.org); NBN Gateway, National Biodiversity Network of United Kingdom (http://www.nbn.org.uk). (XLSX 76 KB)

12862_2014_187_MOESM8_ESM.jpeg

Additional file 8: Figure S2.: Background area used for each species in the background similarity test. A) Agabus brunneus; B) A. ramblae; C) A. rufulus. (JPEG 2 MB)

Additional file 9 Text S1.: Method used to obtain the potential distribution of species. (DOCX 17 KB)

12862_2014_187_MOESM10_ESM.pdf

Additional file 10: Figure S3.: Bivariant plot of the measures of the median lobe of the aedeagus. Open circles, Agabus ramblae; grey diamonds, A. brunneus; black triangles, A. rufulus. Open square: male A. rufulus from Sardinia with an A. brunneus mitochondrial haplotype. (PDF 6 MB)

12862_2014_187_MOESM11_ESM.pdf

Additional file 11: Figure S4.: Ultrametric tree of Agabini obtained in BEAST using only mDNA data, constraining the monophyly of the ingroup and outgroup (genus Platynectes), the genera and the A. brunneus complex. To calibrate the tree we used an a-priori rate of 0.01 substitutions/site/MY (see text for details). Numbers in nodes in black font, posterior probabilities (above 0.5); in red, estimated age. (PDF 1019 KB)

12862_2014_187_MOESM12_ESM.pdf

Additional file 12: Figure S5.: Phylogenetic analyses of the cox1 data. Calibrated tree obtained in BEAST, using a mean rate of 0.02+/–0.001 substitutions/site/MY. Small numbers in nodes, estimated age (Ma), large numbers in nodes, posterior Bayesian probabilities (pp). Negative branches collapsed in polytomies. In red, specimens likely to have introgressed mitochondrial DNA from A. brunneus. (PDF 2 MB)

12862_2014_187_MOESM13_ESM.pdf

Additional file 13: Figure S6.: Phylogram obtained in RAxML with the H3 sequences. Numbers in nodes, posterior probabilities obtained in BEAST (if above 0.5) / bootstrap support (if above 50%). With an asterisk, female from Albacete (SE Spain) of uncertain identity. (PDF 266 KB)

12862_2014_187_MOESM14_ESM.jpeg

Additional file 14: Figure S7.: Representation of the scores of the two first axis of the PCA with all climatic variables. Grey surface, climatic space of the western Palaeartic (background). In colours, climatic space occupied by the three species. (JPEG 320 KB)

12862_2014_187_MOESM15_ESM.jpeg

Additional file 15: Figure S8.: Reconstructed seasonality and minimum temperature of the coldest month during the last glacial interval (upper row) and the last glacial maximum (lower row). Blue circles, current distribution of A. ramblae. (JPEG 384 KB)

Authors’ original submitted files for images

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Hidalgo-Galiana, A., Sánchez-Fernández, D., Bilton, D.T. et al. Thermal niche evolution and geographical range expansion in a species complex of western Mediterranean diving beetles. BMC Evol Biol 14, 187 (2014). https://doi.org/10.1186/s12862-014-0187-y

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12862-014-0187-y

Keywords

  • Thermal niche evolution
  • Cold tolerance
  • Demographic expansion
  • Dytiscidae
  • Western Mediterranean

BMC Ecology and Evolution

ISSN: 2730-7182

Contact us