Phenotypic disparity in Iberian short-horned grasshoppers (Acrididae): the role of ecology and phylogeny

The combination of model-based comparative techniques, disparity analyses and ecomorphological correlations constitutes a powerful method to gain insight into the evolutionary mechanisms that shape morphological variation and speciation processes. In this study, we used a time-calibrated phylogeny of 70 Iberian species of short-horned grasshoppers (Acrididae) to test for patterns of morphological disparity in relation to their ecology and phylogenetic history. Specifically, we examined the role of substrate type and level of ecological specialization in driving different aspects of morphological evolution (locomotory traits, chemosensitive organs and cranial morphology) in this recent radiation. We found a bimodal distribution of locomotory attributes corresponding to the two main substrate type guilds (plant vs. ground); plant-perching species tend to exhibit larger wings and thicker femora than those that remain on the ground. This suggests that life form (i.e., substrate type) is an important driving force in the evolution of morphological traits in short-horned grasshoppers, irrespective of ancestry. Substrate type and ecological specialization had no significant influence on head shape, a trait that showed a strong phylogenetic conservatism. Finally, we also found a marginal significant association between the length of antennae and the level of ecological specialization, suggesting that the development of sensory organs may be favored in specialist species. Our results provide evidence that even in taxonomic groups showing limited morphological and ecological disparity, natural selection seems to play a more important role than genetic drift in driving the speciation process. Overall, this study suggests that morphostatic radiations should not necessarily be considered as “non-adaptive” and that the speciation process can bind both adaptive divergence mechanisms and neutral speciation processes related with allopatric and/or reproductive isolation.


Background
Adaptive radiations, groups that have rapidly diversified from a common ancestor to exploit a wide suite of ecological niches, have intrigued evolutionary biologists for over a century [1,2]. Considered as important biodiversity engines, these bursts of speciation are thought to arise from ecological opportunity in the form of vacant ecological niches that become available due to the colonization of new environments (spatial dispersal) or the acquisition of adaptive innovations that allow the access to novel niche dimensions (ecological dispersal) [3,4]. Spatial and/or ecological dispersal can be driven by modifications of existing environments via climatic changes [4,5]. When lineages first enter these new adaptive zones, morphological evolution should initially be rapid; as niche spaces become increasingly saturated, the rate of morphological evolution would be expected to slow down [6,7]. A core prediction based on the above scenario (the so-called early-burst model) is the existence of early rapid diversification followed by a slowdown in net diversification over time [8][9][10]. Thereby, lineage and morphological diversification are frequently positively correlated [11,12], but the opposite pattern (i.e., a negative relationship) has also been reported [13,14]. Nevertheless, recent studies have shown that processes underlying phenotypic disparity and those generating species diversity can be uncoupled, suggesting that ecological opportunity is not the only diversification force [15][16][17]. For example, speciation by simple geographic isolation can generate a pattern of declining speciation trough time without the intervention of niche-filling processes [18]. On the other hand, speciation bursts can occur in multiple pulses linked to environmental changes (e.g., pulses of orogenic uplift or glacial retreat), thus eroding the diversity-dependent signature of diversification. This may explain why a large number of studies performed on radiating clades have failed to detect early-bursts of phenotypic evolution, despite some of them represent the most classic examples of adaptive radiation, including Darwin's finches from the Galapagos Islands and Anolis lizards of the Caribbean Islands (reviewed in [6]).
Recently, it has been recognized that there are many different types of evolutionary radiations and not all adaptive radiations conform to the early-burst model [19,20]. For example, under a repeated radiation scenario, it is expected that subclades tend to resemble each other in morphological disparity and, contrary to that predicted under the early-burst paradigm, morphological evolution should not necessarily be initially rapid. This kind of adaptive radiation ("iterative radiation") may arise in systems dominated by constraints, such as when evolution within clades is driven by repeated adaptation to similar environments [21][22][23]. This model of radiation exemplifies that episodes of ecological opportunity can be recurring over the evolutionary history of a lineage [24]. Beyond the epithet "adaptive", non-adaptive radiations (or morphostatic radiations) constitute another form of radiation that remains largely unstudied [25][26][27]. Nonadaptive radiations arise through processes that are unrelated to niche exploitation and the resulting species usually exhibit low morphological disparity and allopatric distributions [28]. For example, non-adaptive radiations driven by sexual selection result in new species that are ecologically similar to their ancestors such are the case of Hawaiian Laupala crickets [29]. The boundaries of adaptive and non-adaptive radiations are sometimes diffuse and it has been suggested that both processes are extremes along a continuum [30]. Thus, some lineages may exhibit features consistent with these two concepts [31,32].
Striking radiations showing extraordinary phenotypic divergence, most of them comprising clades with restricted geographic distributions such as islands or ancient lakes, are over-represented in the literature on evolutionary radiations [33] (Fig. 1). Conversely, there is a paucity of studies analyzing the tempo and mode of trait evolution on species-rich clades that present little phenotypic disparity and that are more common among continental radiations. In this sense, recent radiations with limited morphological variability provide valuable case-studies for a comprehensive understanding of the evolutionary mechanisms that drive speciation (Fig. 1). However, despite of morphostatic radiations are likely not the exception but the rule within many groups such as invertebrates, this form of radiation remains largely unstudied. Short-horned grasshoppers (Acrididae), a superfamily that comprises over 6600 species [34], provide an excellent opportunity to study disparification dynamics in an evolutionary radiation exhibiting little apparent adaptive phenotypic disparity. According to a recent study, this group seems to have undergone a significant increase in diversification rate with little extinction, with the major diversification events occurring during the Cenozoic after the Cretaceous-Paleogene boundary [34]. Song and colleagues [34] suggested that the emergence of a new niche space (open grasslands) during the Cenozoic and the subsequent colonization of new habitats may have prompted an explosive adaptive radiation in short-horned grasshoppers. Later, during the Pleistocene, the most speciose acridid subfamily (Gomphocerinae) may have undergone one or several independent radiations due to the evolution of complex species-specific acoustic signals in different clades [35]. Contrary to that predicted under the adaptive radiation model, short-horned grasshopper species do not seem to exhibit remarkable morphological adaptations to different feeding habitats and environments. The existence of a low degree of morphological variation does not necessarily imply an absence of adaptation to different microhabitats. For example, plethodontid salamanders show limited morphological diversification in spite of presenting an extensive array of ecotypes ranging from arboreal to aquatic or fossorial species [36,37]. Nevertheless, acridid grasshoppers seem to be rather conservative in microhabitat usage and, probably, the main differentiating factor among species is the substrate type wherein they perch; some species can be found on the plant canopy whereas others remain on the ground. Although most studies have focused on the relationship between microhabitat and morphology [38][39][40], it is likely that, at a finer scale, substrate type can also impose different selective pressures leading to the evolution of more or less subtle morphological and behavioral adaptations. The ability of species to exploit a range of habitats (niche breadth) could be also correlated to key morphological traits like forelimb length or cranial morphology as evidenced by studies on other taxa (see [41] and references therein).
In the present study, we used a time-calibrated phylogeny of 70 species of Iberian short-horned grasshoppers to test for patterns of morphological disparity in relation to their ecology and phylogenetic history. Specifically, we aimed at testing whether the major morphological changes occurred early in the acridid's diversification history or if random-walk patterns ("morphological drift") or other processes (e.g., evolutionary constraints) have driven the diversification of phenotypic attributes in this group. We also examined the potential role of substrate use and the level of ecological specialization (i.e., niche breadth of a given species) in shaping morphological attributes (locomotory traits, chemosensitive organs and cranial morphology) in this invertebrate radiation. If substrate type and ecological specialization imposes a selective pressure across taxa, it would result in the same solution -phenotypic optima-multiple times, leading to convergence among species even between distantlyrelated clades.

Taxon sampling
The Iberian Peninsula is renowned for its high level of biodiversity and number of endemic species, and also as one of the main refugial areas in Europe during the Pleistocene Ice Ages [42,43]. Most of species included in this study are endemic to the Iberian Peninsula or have a distribution restricted to Iberia, France and North Africa. Our dataset included 70 taxa of short-horned grasshoppers belonging to four different subfamilies; slant-faced grasshoppers (subfamily Gomphocerinae, 43 spp.), band-winged grasshoppers (subfamily Oedipodinae, 17 spp.), spur-throated grasshoppers (Catantopinae, 5 spp.) and 5 spp. belonging to different subfamilies (Calliptaminae ×3 spp., Dericorythinae and Eyprepocnemidinae) (Additional file 1: Table S1). Our dataset accounts for around three quarters of all extant species of acridid grasshoppers that have been recorded in the Iberian Peninsula [44]. Iberian acridids do not descend from a single common ancestor and, thus, they do not constitute a local radiation in this region, but a paraphyletic group [34]. However, it should be noted that several recent studies have pointed out that analyses of adaptive radiation do not require to include all descendants of a particular common ancestor, as such requisite is overly restrictive (see e.g. [19]). In fact, no theory of adaptive radiation predicts monophyly [4].

Phylogeny reconstruction
We sequenced four mitochondrial gene fragments: cytochrome c oxidase subunit 1 (COI), NADH dehydrogenase subunit 5 (ND5), 12S rRNA (12S) and a fragment Fig. 1 The axes of evolutionary radiation. Clades can show (a) strong phenotypic and ecological disparity like Anolis lizards or Hawaiian honeycreepers. Some clades (b) exhibit low morphological disparity despite of they have a wide variety of ecotypes (e.g., Plethodon salamanders) whereas others (c) present considerable phenotypic disparity, but little ecological disparity (e.g. African lake cichlids) suggesting a prominent role of sexual selection (disruptive or diversifying) as driving force in the speciation process. Clades (d) with limited morphological diversification and low ecological disparity (e.g., Acridid grasshoppers, Bythinella spring snails or Muroid rodents) constitute the most extreme cases along these axes and therefore, they have been frequently regarded as non-adaptive radiations containing parts of 16S rRNA (16S). For some taxa we failed to obtain reliable sequences, so we complemented our dataset with sequences obtained from previous studies [45,46]. Sequences were aligned in MAFFT online version 7 [47] and concatenated using Sequencematrix 1.7.8. [48]. We calculated the best-fit models of nucleotide substitution for each of the four genes using jModelTest 0.1.1 [49]. The TIM2 + I + Γ substitution model was selected for 12S, GTR + I+ Γ for 16S, TrN + I+ Γ model for ND5 and TPM3uf + I+ Γ was selected for COI. More details about these procedures are given in [50].
Phylogenetic inference was carried out under both maximum likelihood (ML) and Bayesian frameworks. Maximum likelihood analyses were conducted with two search replicates and 1000 bootstrap replicates using GARLI version 2.0 [51]. Bayesian inference analyses were conducted in BEAST 1.8.0 [52] in order to estimate a time-calibrated phylogeny. We used an uncorrelated lognormal relaxed-clock model and applied a Yule process as tree prior. Two calibration points were used in order to calibrate the tree based on absolute times. We employed as a first calibration point the split between Gomphocerinae and Oedipodinae, estimated to have occurred~100 Mya ago. This estimate is based on dated ancient cockroach fossils [53]. As a second calibration point, we used the divergence between Sphingonotus azurescens (mainland species) and S. guanchus (endemic to La Gomera Island, Canary Islands), whose divergence is estimated to have occurred around 3.5 Mya [54]. Sphingonotus guanchus was only included in BEAST analyses for calibration purposes. Note that dating should be treated with caution as our calibrations are based on previous estimates. We ran analyses for a total of 100 × 10 5 generations, with a sampling frequency of 1000 generations. We ensured that replicated analyses converged (effective sample size values >200) using Tracer 1.4.1. Tree and log files (9500 trees after a 5% burn-in) of the two runs were combined with Log-Combiner 1.4.7 [52], and the maximum clade credibility (MCC) tree was compiled with TreeAnnotator 1.4.7. The obtained tree topologies were rather similar and consistent with previous studies [46]. See [50] for more details.

Morphometric analyses i) Locomotory morphology
Our dataset included 316 preserved specimens representing 3 to 5 male adults for 70 species of short-horned grasshoppers. We gathered a dataset of five morphological attributes of known correlation with locomotor performance and ecology [55]. Morphological measurements included: (1) structural (i.e., head + thorax) body length, (2) tibia length, (3) femur width, (4) femur length, and (5) forewing (tegmina) length. All measurements were taken by the same observer (VGN) using a ZEISS stereomicroscope (SteREO Discovery V.8; Carl Zeiss Microscopy GmbH, Germany). Because disparity in body size and level of sexual size dimorphism has been addressed in a previous study [50] and the former trait (body size) is strongly correlated with other morphological variables, we sizecorrected the remaining four variables using structural body length as size measurement. The relative length and relative width of the femur are inversely correlated (see Additional file 1: Figure S1), so we calculated a width/length femur ratio. As a result, our final dataset was composed of three variables related with locomotory morphology: tibia length, femur width/length ratio and tegmina length. Tegmina are modified leathery forewings present in some insects like grasshoppers. The major role of this structure is that of protecting the hindwings (and therefore, hind-and forewing length are highly correlated; p < 0.001, R 2 = 0.98). These also have an aerodynamic function and thus can be considered an informative trait about the dispersal potential of a species [56]. On the other hand, leg length is thought to determine the jumping performance of insects like grasshoppers, leafhoppers and froghoppers [57][58][59]. Thus, relative femur length can provide information on jumping ability in insects that use a catapult-like mechanism to jump. We performed a phylogenetic principal component analysis (pPCA) on these three variables using the R package phytools [60]. The phylogenetic principal component analysis yielded a single principal component (hereafter pPClm) that explained 53.2% of the variance of our three morphological variables (loadings; tibia length: −0.518, width/length femur ratio: −0.737, tegmina length: −0.886). The positive extreme of pPClm represents species with short wings and tibiae and thicker femora (e.g., Pezotettix giornae, Podisma carpetana), and the negative extreme of this axis represents species with large wings and tibiae and more stylized femora (e.g., Chorthippus jucundus, Stethophyma grossum). We also implemented a nonphylogenetic principal component analysis (PCA) as a recent study has suggested that using phylogenetic PC axes as trait data could bias results [61]. However, we found similar results by using both approaches (correlation pPCA vs. PCA; p = 0.028, r = 0.33).

ii) Antenna length
In conjunction with variables directly related to locomotor performance we also measured antenna length in a similar way as described above. Antenna length was corrected for body size, so we refer to relative antenna length hereafter. Members of the family Acrididae show relatively short and stout antennae in comparison with crickets and katydids. However, the length of antennae largely differs among taxa and these differences may be linked to the diet and degree of ecological specialization of each species [55]. Acridid grasshoppers recognize their host plant by sensorial stimuli, which are perceived through chemoreceptors located on the labrum and antennae [62]. Chemosensory sensilla on antennae express two classes of proteins involved in the recognition of odors and taste, and thus the antennae are thought to play an important role in host orientation and food selection [63]. In this vein, we would expect a greater development of this structure in taxa with a narrower ecological niche (i.e., specialist species) and species preferring to perch on plants, which may display higher sensitivity to plant odour perception [64].

iii) Head shape
Different selective pressures linked to a plethora of factors both intrinsic (e.g., genetic) and extrinsic (e.g., predation risk, habitat structure) can play a role in driving the evolution of head shape. The functional trade-offs that shape head morphology may be related to foraging, locomotory and anti-predator strategies, which may differ among species that move in different microhabitats, for example, between species that climb vertical structures and ground-dwelling taxa [65]. Here, we aimed at examining head shape variation in relation to ecological factors and phylogenetic history. To that end, we photographed a total of 221 specimens (3-4 individuals per species) using a ZEISS stereomicroscope and the ZEISS image analysis software (ZEN2). We established the scale and digitalized landmarks from each photograph using the function 'digitize2D' in the R package geomorph [66]. We selected 14 homologous landmarks which capture the outline of the head and the size and relative position of the eyes (see inset in Fig. 2). These landmarks were subjected to a generalized Procrustes analysis (GPA) wherein all specimens are translated to the origin, scaled to unit centroid size and optimally rotated until the coordinates of corresponding points align as closely as possible [67]. The resulting coordinates in the tangent space represent the head shape (side view) of each specimen. From these coordinates, we calculated an average set of landmark coordinates for each of the 70 species. We averaged each xand y-coordinate, also calculating the species' average centroid size. We then performed a PCA on the aligned Procrustes coordinates using the function 'plotTangentSpace' in geomorph. The first two PC (PChs and PChs2) accounted for 87.1% of shape variation in the sample (60.0 and 27.1%, respectively). We inspected how head shape was related to the log of centroid size (a size estimate) using the function 'proc-D.allometry'. Head shape and centroid size did not show a linear relationship and thus, the allometric component of our data can be considered negligible.

Ecological data
We aimed to examine the evolutionary relationship between grasshopper morphology (locomotory morphology, relative antenna length, and head shape) and substrate type. For this purpose, we categorized each species on the basis of substrate preference (ground vs. plant) from the literature (e.g. [68,69]) and our own personal observations. In addition, we examined the association between our morphological variables and the level of ecological specialization (i.e., the niche breadth of a given species). Ecological specialization refers to the ability of a species to exploit a range of resources and its capacity to use each one as result of evolutionary trade-offs (the "jack-of-alltrades is master of none" hypothesis; [70,71]). Here, we used the 'Paired Difference Index' (PDI) as an estimator of ecological specialization [72,73]. PDI values were calculated from a species-habitat matrix in which we rated the level of association of each species (from complete generalist, 0, to complete specialist, 3) with the nine most common habitats in which the studied species can be found. The PDI is a robust specialization index which takes into account not only the number of resources used by a species, but also the strength of the association between the species and its resources [73]. Scores of specieshabitat association were obtained directly from the literature and our own personal observations. PDI values were computed using the R package bipartite [74].

Phylogenetic comparative analyses
We assessed the phylogenetic signal of our focal variables (locomotory morphology, relative antenna length, and head shape) in order to test if these traits tracked the evolutionary history of the group and were not randomly distributed across taxa. We computed Pagel's lambda (λ; Fig. 2 Scatterplot of the two first principal components of head shape variation (PChs vs. PChs2). Green and brown circles indicate plant-perching and ground-perching species, respectively. Thin-plate spline deformation grids for the two most extreme cases relatives to the overall reference shape (C. wattenwylianus and B. tryxalicerus) are shown. A picture illustrating the position of the 14 landmarks used to characterize head shape variation is also represented in the right-bottom corner of figure [75]) for the first two principal components of head shape (PChs and PChs2), and for the two remaining variables, locomotory morphology (pPClm) and relative antenna length, in phytools [76]. Pagel's λ can take a value from 0, which means that phylogeny has no impact on the distribution of the trait and the values can be treated as independent, to a value of 1, which means that phylogeny fully predicts the distribution of the trait. We then mapped the evolution of these continuous traits on our phylogeny by using the function 'contMap' in phytools, which estimates states at internal nodes using ML and interpolates the states along each edge following [77].
We examined the association between the two ecological factors (substrate type and degree of ecological specialization) and our morphological variables (relative antenna length, PChs and pPClm). First, in order to test if substrate type ("ground-perching" vs. "plant-perching"; categorical variable) has influence on these variables, we performed phylogenetic ANOVAs using the function 'phy.anova' (10,000 simulations) in the R package geiger [78]. Because substrate type largely varies among subfamilies (most gomphocerine species are plant-perching species whereas most of the remaining species stay on the ground) we performed a second phylogenetic ANOVA only including Gomphocerinae species (n = 43). Thereby, we assessed the influence of substrate type on morphology while controlling for the non-independence of the data because of shared ancestry. Secondly, we analyzed the relationship between our morphological variables and the level of ecological specialization (PDI index; continuous variable) using phylogenetic generalized least squares (PGLS λ ). PGLS regression analyses were performed using the caper package [79] and graphically visualized by means of phylogenetically independent contrasts (PICs) computed using the PDAP:PDTREE module in Mesquite v.3.04 [80].

Disparity analyses
In order to visualize the relationship between phylogeny and taxon distribution in the morphological space, we built a 'phylomorphospace' , which projects the branches of a phylogenetic tree into a 2D morphospace defining trait variation among species [81]. When subclades occupy limited morphospace it means that most disparity is accounted for by early divergence (early-burst of trait evolution), while if subclades occupy large regions of the morphospace (i.e., stronger overlap) it means that morphological diversity ('disparity') is largely explained by recent divergence. We generated a phylomorphospace by plotting head shape against locomotion morphology (PChs vs. pPClm) using phytools [60]. We also generated a second morphospace from the two first principal components for head shape (PChs vs. PChs2).
We investigated how morphological disparity has accumulated over the group' evolutionary history. We evaluated the rate of morphological evolution in relation to lineage diversification by means of disparitythrough-time (DTT) analyses [82], as implemented in the geiger package [78]. We compared observed relative disparity with average phenotypic disparity simulated under Brownian Motion (1000 simulations) in order to test whether disparity was accumulated during the early or recent history of the group. When the rate of morphological evolution is constant, as expected under BM, the DTT curve is expected to decline linearly toward zero through time. Under an early-burst model, disparity will decline sharply and much earlier, while if evolution within subclades is fast the observed DTT curve will fall above the Brownian profile. In order to quantify the magnitude of the difference between the null disparity profile computed under BM, and the observed disparity profile from our data, we computed the so-called morphological disparity index (MDI) statistic using the geiger package [78]. A positive MDI value is indicative of greater than expected subclade disparity (i.e., disparity distributed primarily within subclades) whereas negative MDI values indicate greater than expected subclade disparity (i.e., disparity distributed primarily among subclades). A negative MDI is characteristic of adaptive radiations, as rapidly diversifying taxa are expected to evolve distinct morphologies in response to new adaptive zones and slow once niches are filled. DTT analyses were performed twice; by including all taxa and only including species of the subfamily Gomphocerinae. We then ran a node height test [83], which tests for accelerations or decelerations in trait evolution, by comparing the independent contrasts for a trait with the respective node height (i.e., relative age). A significant negative relationship between absolute contrast value and node age supports the hypothesis of adaptive radiations as it would imply that species are dividing niche space more finely through time, consistent with a niche-filling model. Node height tests were performed using the R package caper [79].

Tempo and mode of morphological evolution
Evolutionary models are conceived to infer the different possible processes shaping phenotypic evolution and provide a method of testing different predictions regarding the tempo and mode of evolution. We assessed the fit of four evolutionary models to our phenotypic variables using the functions 'fitContinuous' and 'OUwie' implemented in the R packages geiger [78] and OUwie [84], respectively. We fitted three different single-rate models; Brownian-Motion (i.e., diffusive drift), Early-Burst (i.e., exponential declining of evolutionary rates) and single-optimum Ornstein-Uhlenbeck (i.e., bounded evolution around a single phenotypic optimum), and a multi-peak Ornstein-Uhlenbeck model. Brownian-motion (BM) operates under the assumption that trait evolution proceeds a random walk wherein trait variance across lineages accumulates proportional to time (single-rate model; [77]). The Early-Burst (EB) model predicts rapid morphological disparity early in the radiation, followed by a slowdown in the diversification rate as ecological niches are filled over time [82]. Support for the single-optimum OU model (OU1) would imply that there is a single phenotypic optima (θ) for all taxa. This model is often associated with a process of stabilizing selection in which variation of phenotypic traits revolves around one or more stationary peaks [85]. Lastly, we assessed the fit of a multivariate OU model (OUVM) with separate morphological optima and separate random walk variances (σ 2 ) for each substrate guild (1: plant, 2: ground), and one global parameter (α), which determines the strength of selection towards those optima. For the multi-peak OU model we determined the possible ancestral substrate regimes along the internal branches of the phylogenetic tree using the Stochastic Mutational Mapping on Phylogenies (SIMMAP) tool in phytools, and sampled 500 character histories in order to incorporate evolutionary uncertainty. We compared the fit of the BM, EB, OU1 and OUMV models using the Akaike information criterion corrected for small sample size (AICc), which can be employed to compare models that differ in the number of parameters and therefore have non-comparable likelihoods.

Ecology-morphology association
We found a marginally significant correlation between the relative length of antennae and the level of ecological specialization (PGLS; estimate: 0.304 ± 0.157, t = 1.92, p = 0.057); specialist species tend to have larger antennae than generalist species (Fig. 3). There was no significant association between either head shape (PChs) or locomotory morphology (pPClm) and the level of ecological specialization (t = 1.54, p = 0.128, and t = −1.37, p = 0.174, respectively).
We found a significant effect of substrate type on head shape (PChs) (phylogenetic ANOVA, F 1,68 = 20.81, phylop = 0.029) indicating that plant-perching species exhibit a more conical head in comparison with ground species. However, when restricting our analyses to the Gomphocerinae subset, we did not find significant differences in head shape between both groups (F 1,42 = 3.81, phylop = 0.25). This suggests that differences in head morphology between substrate types are mainly due to most of Oedipodinae being "ground-dwelling" taxa whereas most Gomphocerinae species are adept climbers that prefer to perch on plants instead of remaining on the ground. Regarding locomotory morphology (pPClm), we detected marginally significant differences between plant-perching and ground-dwelling species (F 1,68 = 14.60, phylop = 0.077); the former tend to exhibit larger tibiae and larger and thinner femora than ground-dwelling species, which possess a more stout morphology (Fig. 4). This difference remained marginally significant even when restricting our analyses to the Gomphocerinae subfamily (F 1,42 = 9.00, phylo-p = 0.074); gomphocerine species that prefer to stay on the ground exhibit a slightly different locomotory morphology (short wings, thicker legs) in comparison with gomphocerine perching on plants. There were no significant differences in relative antenna length between ground and plant-perching species for either the entire dataset (F 1,68 = 2.93, phylo-p = 0.31) or the Gomphocerinae subset (F 1,42 = 0.33, phylo-p = 0.25).

Disparity analyses
Phylomorphospace plots of PChs versus pPClm and of PChs versus PChs2 indicate a moderate phylogenetic structuring of taxon distribution in the phenotypic space (see Additional file 1: Figure S4). We projected our phylogeny onto a plot defined by the first axis of head shape (PChs) and the PC axis representing variation in locomotory morphology (pPClm). We found that Catantopinae and Calliptaminae occupy a unique region of morphospace at the positive extreme of pPClm, characterized by relatively small and stout limbs, whereas Gomphocerinae and Oedipodinae species showed a broader distribution along the morphospace (Fig. 5). Within the 2-dimensional space defined by PChs versus PChs2, we observed that the first axis, which represents variation between an elongated (slant-faced) to a straight head shape (60% of variation), separates most of gomphocerine species from the rest of taxa (Fig. 2). When discerning between both substrate types, we observed that plant-perching species (as most Gomphocerinae are) are more abundant in the right-side of the morphospace while ground species (23% of Oedipodinae falls within this category) tend to occupy the left-side (Fig. 4).
The DTT analyses yielded positive MDI statistics for all morphological variables (relative antenna length: 0.166, head shape PChs: 0.070, locomotory morphology pPClm: 0.137; Fig. 6), however the differences between both profiles were not statistically significant in either case (all p-values >0.05). It means that trait evolution did not significantly deviate from a BM model, which prevents us to state that disparity is concentrated within subclades (i.e., that closely related species differ considerably in morphology). For head shape, we found a higher pulse of trait disparification during the Paleogene, which is in agreement with previous studies in which it has been suggested that major lineages of Acrididae went through a major radiation in the Cenozoic [34]. In all cases we observed a departure of the observed DTT curves from BM expectations beyond the 95% confidence interval during the Pleistocene suggesting rapid acceleration in morphological diversification during this period. When restricting our analyses to this period (i.e., by including only the Gomphocerinae subfamily) we obtained again positive MDI values in all cases (relative antenna length: 0.355, head shape PChs: 0.841, locomotory morphology pPClm: 0.456), which implies that due to niche evolution subclades overlap (i.e., disparity is within subclades), and all contain a significant proportion of variation found through the Gomphocerinae group at a given time. The node-height test resulted in a positive but nonsignificant relationship between the absolute values of standardized length contrasts and node age for all morphological traits (PChs: b = 0.011 ± 0.008, t = 1.36, p = 0.18; relative antenna length: b = 0.013 ± 0.007, t = 1.72, p = 0.09; pPClm: b = 0.018 ± 0.01, t = 1.86, p = 0.07).

Tempo and mode of morphological evolution
Results from the multivariate model fitting analysis support stabilizing selection (OU1, OUVM) as the most plausible evolutionary scenario for our morphological data (summarized in Table 1). According to the obtained AICc values, the OUMV model provides the best fit to head shape evolution (PChs) followed by the BM model and lastly the OU1 and the EB models (Table 1). This means that head shape diversification may be driven by adaptive constraints, although a constant rate model could not be discarded. Indeed, when models of head shape evolution were tested exclusively on the Gomphocerinae, BM    Table 1 Comparison of four evolutionary model fits, random-walk variances (σ 2 ) and primary trait optima (in parentheses) for the three morphological variables: head shape (PChs), locomotory morphology (pPClm) and relative antenna length ΔAICc is the model's mean AICc minus the minimum AICc between models. Bolded rows represent the best fit model as indicated by the lowest AICc score. Estimated phenotypic optima (θ) for ground and plant-perching regimes are shown where applicable. Analyses were performed for (a) all species and (b) only including the Gomphocerinae subset provided the best fit for the observed pattern ( Table 1). Divergence of locomotion attributes was also best described by an OUMV model (Table 1), which suggests that acridid body morphology is tied to substrate type. This result was confirmed when restricting our analyses to the Gomphocerinae dataset, which allowed us to distinguish ancestral from derived similarity (Table 1). However, a single-peak OU model could not be excluded in this latter analysis. Regarding relative antenna length, the evolution of this trait was better described by a single-peak OU model either considering all species or restricting the analyses to Gomphocerinae (Table 1).

Discussion
Like some examples of non-adaptive radiations (Plethodon: [36], Rattus: [32]), acridid grasshoppers tend to retain a rather conserved body plan with little overt ecomorphological specialization among taxa (Fig. 6). However, we observed a greater morphological disparity within Oedipodinae in comparison with that observed for the Gomphocerinae, suggesting that the level of morphological resemblance among species differ between the two subfamilies. Gomphocerinae (slant-faced grasshoppers) constitute a very recent radiation, which is suspected to be mostly the result of divergent evolution of isolated populations induced by climatic oscillations during the Pleistocene [42]. Some gomphocerine species spread out from southern refugia and expanded their ranges northwards during interglacials, whereas many others (montane species) became restricted to high altitudes ("sky-islands") due to their limited dispersal capacity [35,86]. Hence, many Iberian species probably arose in association with these mountain ranges (Sistema Central, Betic ranges, Picos de Europa; [87]). At a smaller scale, the evolution of elaborate acoustic signals as premating reproductive isolation mechanisms (via sexual selection) may have allowed gomphocerine grasshoppers to diversify in regions in which ecologically similar and phylogenetically related species coexist [46,86,87]. Thus, it is likely that both divergence in allopatry [42,86] and prezygotic isolation mechanisms [46,88,89] have played a crucial role in promoting speciation in this subfamily. However, our results should be interpreted cautiously because the present study assesses the dynamic of morphological evolution among a small fraction of the Acrididae' overall diversity. Thus, our findings may reflect local processes (i.e. factors that have affected which species can occur together in the same regional theatre) rather than global mechanisms that operate across the entire clade.

Patterns of morphological disparity through time
Our results do not support the existence of an early-burst of phenotypic diversification, adding to the number of study-cases showing no support for this evolutionary model [6,7,90,91]. Obtained MDI values indicate that substantial amount of variance is clustered within subclades, pointing out to the existence of low phylogenetic niche conservatism. That is, there was no tendency of closely related species to be more similar to each other in ecomorphological traits than they are to more distant relatives as expected if some processes had constrained niche divergence among phylogenetically close species [92,93]. This is in contrast to that predicted for taxa that experience early-bursts of diversification, where the partitioning of morphological disparity through time (disparification sensu [94]) is expected to be low [10,82,95]. Overall, DTT plots showed a general tendency for relative disparity in phenotypic attributes to decrease over time, but with two main pulses of increases in disparification during the Cenozoic and the Pleistocene. The observed pattern is not consistent with that predicted under an iterative radiation scenario wherein it should be expected the appearance of repeated diversification peaks [23]. The fact that acridids partitioned morphological disparity within rather than among clades suggests that acridid lineages did not evolve along distinct morphological trajectories through time.
Thus, it is unlikely that acridid lineages explored different adaptive zones, which makes sense taking into account the limited morphological variability they exhibit.

Evolutionary models of morphological divergence
Overall, results from the multivariate model fitting analysis summarized in Table 1 favored an Ornstein-Uhlenbeck model, indicating bounded evolution around a single (OU1) or two phenotypic optima (OUVM). Specifically, in relation to head shape, we found that a single model does not uniquely explain the evolution of this trait for the entire dataset; BM and OUMV were similarly informative. However, when models of evolution were tested exclusively on the Gomphocerinae subset, BM was the most informative model. It suggests that head shape evolution displays an idiosyncratic component and that, irrespective of substrate on which they rest, Gomphocerinae species have a more conical head shape in comparison with band-winged grasshoppers due to their phylogenetic legacy. On the other hand, our results indicate that locomotory morphology in Acrididae is inconsistent with a BM process, and has not evolved under a constant rate over time. A multivariate-peak OU (OUMV) model provided the best fit in both cases (when considering all taxa and after including only Gomphocerinae), which means that the evolution of this trait revolves around adaptive peaks. In this context, it is possible that biomechanical constraints, which underlie much of locomotory morphology variation, impose restrictions on morphological diversification [96][97][98]. Specifically, our results suggest that stabilizing selection pulls the locomotory morphology towards two convergent adaptive optima (θ plant , θ ground ) during the group's evolutionary history. A key difference between the two substrate type regimes in terms of locomotory attributes relies on the different degree of development of femora (which has influence on jumping performance) and wings (which determines flight capability) that each guild of species exhibit (see more below). Interestingly, we observed a negative relationship between relative forewing length and relative femur width/length ratio in Gomphocerinae (see Additional file 1: Figure S5) pointing out to the existence of a trade-off between flying-and jumping-based dispersal investments in this group. With regard to antenna length, a single-peak OU model was the model that provided the best fit to the observed data, suggesting that antennae evolution has remained constrained by directional selection.

Ecology-morphology association
We investigated the role of substrate type and level of ecological specialization in driving morphological evolution in our study group. The variation in structural characteristics within the habitat (substrate type, size and incline) has direct effects on animal locomotion and can lead to the evolution of morphological and behavioral adaptations (e.g., anti-predator strategies) [99]. For example, stick insects (order Phasmatodea) of New Caledonia and New Guinea ("land lobsters") are grounddwelling species and show a distinct ecomorph; they are flightless and exhibit a stocky body form and thick hind legs. This is in contrast with the majority of stick insects, which are solitary canopy-dwellers and show an instantly recognizable morphology (stick-like bodies with large and thin legs) [100]. In particular, substrate type may be especially important for small animals that jump using a catapult-like system such as grasshoppers and leafhoppers [101]. Here, we found that plant-perching grasshoppers tend to exhibit a slender morphology in comparison with those species that prefer to stay on the ground, which show a more compressed body (shorter and thicker femora). This result is consistent with that reported by [58] in a study with leafhoppers. These authors suggested that insect species that jump from plants have longer legs because it reduces the amount of energy lost to bending the leaf, whereas species that jump from a stiffer substrate can have shorter-legs [58]. Hence, in species living on plants, natural selection would have favored long hind legs because these allow short preparation times when an emergency jump is necessary [102]. The existence of moderate differences in morphology between ground-dwelling and plant-perching species was confirmed when restricting our analyses to the Gomphocerinae subset. Despite of members of the subfamily Gomphocerinae exhibit low variability in morphology (and as a result of this phenotypic conservatism some species can be almost only distinguished by song attributes), gomphocerine species perching on plants tend to present a less rotund morphology -larger and thinner legs-in comparison with gomphocerine species that prefer to stay on the ground, which show squat bodies and short limbs. Our results suggest that there are two possible ways to achieve superior jumping capabilities: to evolve longer limbs or to evolve more muscular limbs (as stronger legs will produce more acceleration) [102]. This finding is in agreement with that reported by [96] in a study with anole species, which also seem to face trade-offs that prevent them from simultaneously optimizing different skills of jumping ability (larger vs. thicker legs).
Head shape differed significantly between ground-and plant-species across taxa. Plant-perching species show a sharper head than ground-species, which may be favored by natural selection as this head shape facilitates camouflage against potential predators by making it difficult to see the mimic against the surroundings (i.e., background matching) [103]. However, this difference between both categories become non-significant within the Gomphocerinae, which implies that differences in head shape between ground-and plant-species have to do with its different phylogenetic history per se rather than with the possible existence of a substrate type effect. In bandwinged grasshoppers (subfamily Oedipodinae), the orientation of the face is usually nearly vertical (straight head shape) whereas species of the subfamily Gomphocerinae (also known as slant-faced grasshoppers) show a sharper profile, that is, a more conical head.
Lastly, we also found a marginally significant association between the relative length of antennae and the level of ecological specialization. Larger antennae can bear a higher number of chemoreceptors, which would be favorable for grasshoppers in the search for shelter or food [104]. Hence, in those species that can thrive only in a limited range of habitats or have a limited diet (specialist species), the development of sensory organs may be favored, whereas in species with a broader ecological niche breadth (generalist species) selection for larger antennae may be weaker. Thus, the length of antennae in short-horned grasshoppers could be employed as a proxy for the degree of ecological specialization in terms of habitat requirements of a given species.

Conclusions
Short-horned grasshoppers exploit a variety of habitats, but lack the extreme specializations ("ecomorphs") observed in other taxa (see Fig. 1). Accordingly, at a first glance, it seems that acridid grasshoppers exhibit a conservative and "all-purpose" morphology allowing them to perform equally well under different environments, which is contrary to one of the premises that are frequently used to catalogue a radiation as adaptive [19,105]. However, in this study we provide evidence that natural selection could lead to a convergent body plan in very distant clades in this recent radiation [106]. Specifically, we found substrate type-related macroevolutionary variation in locomotory attributes and head morphology. The observed pattern was consistent with an evolutionary scenario with two adaptive peaks (selective agents) related to life-form (i.e., substrate type; plant vs. ground). Thus, these findings suggest that, although this radiation does not meet the conditions to be classified as adaptive sensu stricto (at least in comparison with textbook examples of rapid adaptive radiations such as Madagascan vangas and tetragnathid spiders on the Hawaiian Islands), it provides evidence that natural selection can act on a very small scale, and therefore it is not always detectable. So, even in a priori morphostatic radiations, adaptive processes can lead to subtle phenotypic variation. Hence, our results support the notion that ecological factors play a key role in evolutionary processes by acting as selective mechanisms and by imposing constraints that shape morphological traits [107]. In sum, this study highlights that natural selection, and not genetic drift, constitutes the main factor driving the disparification process in this recent radiation despite of the lack of apparent phenotypic and ecological variability. This calls into the question the use of term "non-adaptive" as synonymous with "morphostatic" since the total absence of an adaptive component seems to be unlikely in any radiation.

Additional file
Additional file 1: Table S1. Information on substrate type (plant-or ground-perching) and niche breadth estimated in form of categorical (generalist vs. specialist species) and continuous (PDI: 'Paired Difference Index' values) variables for the 70 grasshopper taxa included in the present study. Figure S1. Relationship between relative femur length and relative femur width. Dots show mean values for each acridid species and colors indicate clade membership (blue: Gomphocerinae, red: Oedipodinae, yellow: Calliptaminae-Dericorythinae-Eyprepocnemidinae, grey: Catantopinae). The arrow denotes a case (B. tryxalicerus) that deviates remarkably (i.e., an outlier) from the general trend. Figure S2. Maximum likelihood ancestral reconstruction of head shape variation (PChs) in Iberian short-horned grasshoppers. Figure S3. Maximum likelihood ancestral reconstruction of locomotory morphology (pPClm) variation. Figure S4. Phylomorphospace projection of acridid species on the first two principal components of head shape variation, which account for 87% of the variance. For illustrative purposes, B. tryxalicerus, the most extreme case in both axes, was not represented. The inset shows the phylogenetic relationships among species including B. tryxalicerus, which is highlighted in color blue. Figure S5. Relationship between relative forewing length and femur width/length ratio. Dots show mean values for each acridid species and colors indicate clade membership (blue: Gomphocerinae, red: Oedipodinae, yellow: Calliptaminae-Dericorythinae-Eyprepocnemidinae, grey: Catantopinae). Figure S6. Head shape variation (PChs) plotted against locomotory morphology variation (pPClm). Dots show mean values for each acridid species and colors indicate clade membership (blue: Gomphocerinae, red: Oedipodinae, yellow: Calliptaminae-Dericorythinae-Eyprepocnemidinae, grey: Catantopinae). Figure S7. The time-calibrated phylogeny for 70 species of short-horned grasshoppers used in the present study. (DOCX 517 kb)