Contributions to Zoology, 70 (4) (2002)
Landmark-based morphometric analysis of two sibling species of the genus Asida (Coleoptera, Tenebrionidae)
Keywords: Geometric morphometry, Thin-Plate spline analysis, canonical variates analysis, species discrimination, allometry, Tenebrionidae (Coleoptera).
The case described here analyses morphological change at the boundary between ecological and evolutionary scales. The size and shape of 8 populations of two sibling species of tenebrionid beetles (Asida planipennis and A. moraguesi) are analysed using landmark-based methods. The two species differ in size, shape and in allometric trajectory. Thin-Plate Spline Analysis (TPSA) combined with Canonical Variate Analysis (CVA) reveal the specific shape changes that allow the best inter-species discrimination. These changes involve the outline of the posterior margin of the pronotum. Moreover, the landmark-based method provides useful tools for interpreting the intra-species variability of some continuously varying morphological characters. In the case of A. planipennis, size and shape are correlated at the inter-population level, but are independent at the intra-population level. Moreover, size and shape do not reflect any spatial (i.e., geographical) structure or phylogenetic inertia at the inter-population level. These facts favour site-specific environmental conditions as the main cause of shape and size variability in this species. One environmental variable is suggested to be the cause of the inter-population morphological differences detected.
Faunas from oceanic and oceanic-like islands show a high rate of endemicity. Oceanic-like islands are those that have been previously connected to continents but are faunistically similar to oceanic islands because the island-continent connection occurred in the deep past or did not promote significant faunal transfer (Alcover et al., 1998). The fauna from the Balearic islands, as an archetypal example of oceanic-like islands, is characterised by a relevant percentage of endemic species that are typically flightless and with reduced dispersal potential (Pons & Palmer, 1996; Palmer et al., 1999).
Sea-level of the Mediterranean has experienced notable changes during the glacial-interglacial pulses (Cuerda, 1975; Tuccimei et al., 2000). Some inter-island faunal transfers within the Balearic Archipelago were possible during marine regressions, and subsequent marine transgressions implied isolation events. The outcome is a puzzling pattern of inter-population differentiation. A well-known example is that of the lizard Podarcis lilfordi Günther 1874, with more than 20 subspecies described. Most of islets around Mallorca and Menorca Islands harbour a population of this species displaying some specific morphological differentiation (Mayol, 1985; Salvador, 1979).
Here I analyse the body size and shape of 8 Balearic insular populations of flightless endemic beetles of two sibling species of the genus Asida Schaufuss Latreille 1802 (Coleoptera, Tenebrionidae). They display some morphological differentiation at both intra- and inter-population level. The above-depicted biogeographical scenario seems suitable for analysing two relevant aspects of morphological patterns. First, I examine the inter-species discrimination considering both shape and allometric trajectories. Second, I check the significance of the effects of common history and geographical location on the pattern of inter-population shape differences. Inter-population shape similarity can arise because two populations are (phylogenetically) closely related or (geographically) closely located. Or it can result of the fact that populations share a similar set of site-specific environmental variables. The effects of spatial structure between sampling sites and inter-population phylogenetic inertia are here explicitly tested whereas the effects on body size of one site-specific variable (island size) are elucidated elsewhere (Palmer, submitted).
Material and methods
Collecting samples and measuring size and shape
I studied 105 males of two sibling endemic species of the genus Asida (A. planipennis Schaufuss 1869 and A. moraguesi Schaufuss 1879) from 7 islands of the Balearic Archipelago (W Mediterranean; see Fig. 1). Sample sizes and island sizes are detailed in Table 1. The specimens analysed are deposited in the Museu de la Naturalesa de les Illes Balears (Palma de Mallorca, Spain). I studied males only in order to simplify the problem and to focus only on inter-population differences (males and females show morphological differences). In addition, females are apparently scarce, and the few available individuals precluded performing a separate analysis. The taxonomy follows Español, 1954.
Body size and shape were recorded using the methods of Geometric Morphometrics (Rohlf & Marcus, 1993; http://life.bio.sunysb.edu/morph/). Conventional analyses of shape are characterised by the application of multivariate statistical methods to a set of variables (usually lengths). On the contrary, Geometric Morphometrics records shape as two- (or three-) dimensional coordinates. The relationships between landmarks (i.e., the shape) are preserved by fitting them into an appropriate function after rotation, translation and scaling. The parameters of the fitted function are then susceptible to be used as variables in standard multivariate statistical analysis (Rohlf & Marcus, 1993). The geometry (shape and size) of the pronotum was captured by recording two-dimensional coordinates of 8 (Fig. 2) landmark points (Rohlf & Marcus, 1993).
Two-dimensional coordinates of these landmarks were collected from every specimen using an image-analyser composed of a Wild dissecting microscope, a video camera, and VIDAS21 software (Kontron Elektronic, Munich). Specimens’ size was estimated by its centroid size of the pronotum (i.e., the squared root of the sum of the squared distances from each landmark to the centre of gravity, Rohlf & Marcus, 1993). Centroid size was determined using the TPSReg program (Rohlf, 1998a). Every specimen was measured three times and the three values of the centroid size obtained were then averaged. The centroid size of a population is determined by averaging the centroid size of every specimen measured.
Pronotum shape was captured with the same landmarks mentioned above (Fig. 2). Landmark configuration derived from the three replicate measures of the same specimen were translated, rotated and scaled (to a centroid size =1) in order to super-impose them (i.e., to overlay configurations so that homologous landmarks match as closely as possible). The super-imposition procedure implies the determination of a consensus landmark configuration by an iterative procedure based in the Procrustes metric (Bookstein, 1991). Procrustes distance is approximately the squared root of the sum of the squared differences between the positions of all landmarks of one object to the reference (Slice et al., 1996). Such a procedure minimizes shape differences between the consensus configuration and the configuration of all objects being superimposed. Symmetry axis was then used to reflect the left side of the three-replicate consensus on its right side. The consensus configuration for these two hemipronotums was then determined (Klingenberg & McIntyre, 1998).
The resulting 105 hemipronotum consensus configurations were the input to undertake a Thin-Plate Spline Relative Warps Analysis (Bookstein, 1991; Rohlf & Marcus, 1993; Rohlf, 1998b). This analysis fits between-landmark differences to a mathematical function and provides useful information at two levels: 1) displaying shape changes as grid deformations (i.e., splines), and 2) providing a parameter matrix (W matrix, i.e., the warp scores) ready to be analysed by standard multivariate statistical methods. Shape decomposition performed by TPSRW allows making comparisons of shape, size and shape taken together and with different types of shape, i.e., the non-uniform component or local-scale change, and the uniform component or general change. Here I focus in local-scale changes.
The obtained parameter matrix (W matrix) for 102 specimens from 7 populations (Na Plana population, with only three specimens, was excluded) was used to test the significance of inter-population shape differences. In that case, the W matrix was used as input for a Canonical Variate Analysis (NTSYS; Rohlf, 1993) to obtain a matrix of Generalized Mahalanobis Distances (i.e., distances measuring the inter-population shape differences). Testing for the significance these inter-population shape differences was made using unbiased Mahalanobis distances (an adjustment for sample sizes; Marcus, 1993).
The existence of allometric patterns was analysed by regressing shape on ln-transformed size (using TPSreg; Rohlf, 1998a). This procedure was completed including and excluding A. moraguesi in the analysis. In addition, separate analyses were made 1) considering the full set of specimens and 2) the consensus for every population only.
Spatial structure and phylogenetic inertia
Inter-population shape similarity can arise because populations are (geographically) closely located or (phylogenetically) closely related, or because they live under the same set of environmental features. The usual pattern of spatial autocorrelation arises
when the values for a single variable (e.g., size) are similar for populations located within a narrow range, whereas tend to be dissimilar for distant populations. The existence of this pattern can be evidenced, for example, by clustering the interpopulation (geographical) distances into a few classes. The Moran’s I statistic (e.g., Legendre & Legendre, 1998) measures similarity (i.e., correlation) between each observation (e.g., the size of one population) and all the other observations within the same distance class. Here, the observed interpopulation distances were clustered into three classes, and three Moran’s I values were obtained. Significance of the observed Moran’s I values was assessed by random permutation of the observed distance matrix (1000 iterations). Tests were one-tailed because low morphological similarity at short distance has no biological sense in the current scenario. Statistical testing of the whole set of Moran’s I values (known as correlogram; Legendre & Legendre, 1998) implies multiple comparisons. Therefore, progressive (Legendre & Legendre, 1998) Bonferroni correction was adopted.
Multivariate Mantel correlogram (Legendre & Legendre, 1998) was developed as a way to compute a correlogram for multivariate data. This is the case for shape that is here defined by the Mahalanobis distance matrix (derived from a canonical variate analysis of four shape variables corresponding to the two first partial warps; see above). Multivariate Mantel correlogram uses normalized Mantel Z statistic. Similarly to the case of size, geographic distance matrix was recoded in to three categories and three separate Mantel’s Z values were determined. As an example, the computational procedure for the first distance class implies to code as 1 all neighbouring populations and the remainder values are set to zero. Mantel’s Z statistic should be considered comparable to the Moran’s I statistic excepting for the fact that positive indices of the last statistic indicates positive autocorrelation. This is not the case for Mantel’s Z because a distance matrix generates a negative coefficient in the case of positive autocorrelation (i.e., inverting the signal).
Similarly to the case of geographic distance, the null hypothesis of phylogenetic independence among samples was tested by a correlogram of the estimated pair-wise inter-population distances (Gittleman & Kot, 1990; Diniz-Filho et al., 1998). I lack an explicit (e.g., DNA-based) phylogeny for the populations considered. However, the history of the inter-island isolation process can be deducted from paleogeographic scenarios. Thus, the sea level dropped below 100 m in the Western Mediterranean basin during the last Ice Age (»18.000 yr BP). At this time all the populations considered became connected (Cuerda, 1975). The northern mega-island of the Balearics (i.e., Mallorca, Menorca, Cabrera and the inshore islets; Fig. 1) was probably inhabited by a common biota (Palmer et al., 1999). The sea-level rise at the end of the last Ice Age determined a hierarchical pattern of isolation. Therefore, the last faunal transfer between two islands in the case of lack of dispersal since the isolation event (i.e., the inter-population phylogenetic distances) can be estimated from the bathymetry and the sea-level eustatic curve. Lack of dispersal could be a realistic assumption for Asida beetles because 1) they are flightless, 2) their large body size precludes passive dispersal as aerial plankton, 3) ecological constraints preclude ornithochory and other passive dispersal mechanisms, and 4) all Western Mediterranean islands (or island groups) that remained isolated during the sea-level drop at the last Ice Age are populated by specific endemic Asida species (Palmer et al., 1999). Under the above assumption, the existence of phylogenetic inertia in the variables of interest (centroid size and shape of the pronotum) have been tested using the bathymetric pair-wise matrix as an estimate of pair-wise phylogenetic distances.
Allometric trajectories and the taxonomic status of Asida moraguesi.
Canonical Variate Analysis of the W matrix of the whole data set shows that all populations of A. planipennis cluster well apart from the single population of A. moraguesi (Fig. 3). The two first canonical axes explain 91% of the cumulative between-groups variance. However, the first one explains alone up to 85% of this variance. Therefore, the first canonical axis defines the main trend of shape change. Note that each axis can be interpreted as a shape gradient between two extremely deformed splines. These splines are obtained by deformation of a regular grid that corresponds to the general consensus configuration (i.e., the “average” shape). The regions of the splines showing large deformations (between the extreme splines) indicate the shape changes that allow the best discrimination. These changes concern the outline of the posterior margin of the pronotum. This margin is extremely necklined in A. moraguesi. In the case A. planipennis, it ranges from moderately necklined to rounded (na Moltona population). In order to facilitate the interpretation of the shape gradient, the splines corresponding to the extreme populations are superimposed (by eye) on the image of a specimen from each of the two extreme populations (Fig. 3). The splines were obtained by regressing the shape on the first canonical axis using TPSReg (Rohlf, 1998a). Consistently with the pattern reflected in Fig. 3, unbiased pair-wise Mahalanobis distances within the canonical space indicates significant shape differences between A. moraguesi and all populations of A. planipennis (Table 2).
Fig. 3. Minimum spanning tree of population centroids superimposed over the plot of the first two canonical axes. Splines illustrating the extreme morphological change along the first canonical axis were obtained regressing shape on the canonical scores (first axis). They were superimposed on the image of one specimen from the na Moltona population (right), and A. moraguesi (left). Shape differences between A. moraguesi and A. planipennis are statistically significant (Table 2).
In addition, A. moraguesi and A. planipennis seem to display different allometric trajectories. This pattern is evidenced both when considering consensus configurations only and when considering all measured specimens together. In the two cases, shape appeared to be size-dependent, but A. moraguesi falls clearly below the trend outlined for A. planipennis (Fig. 4). In the case of the size-shape relationship, the splines on the shape axis were obtained by regressing shape on size (using TPSReg; Rohlf, 1998a).
Fig. 4. Scatter plot of shape (first partial warp only) against centroid size. Two related regression analyses are shown. A, consensus configurations from the 8 populations regressed against size (splines showing the shape changes explained by size are added). B, all specimen configurations are included in the regression analysis. In both cases, A. moraguesi falls below the trend outlined by A. planipennis.
Shape-size relationship within A. planipennis
The foregoing results justify removing the single population of A. moraguesi from the analysis in order to get the allometric pattern free of possible distortions due to species-specific effects. After removing A. moraguesi, size showed a significant effect on shape at inter-population level (i.e., when the shape of the consensus specimen from each population was regressed on centroid size). Size explains 54% of shape variability measured as procrustes distance between the general consensus configuration and the consensus configuration of every population. Generalized Goodall F-test (Rohlf, 1998a) based on such distances suggests a significant effect of size on shape (F = 5.8941, df = 6, 30, Prob. = 0.0004). In contrast, size and shape seem to be uncorrelated at intra-population level. Extensive results are not reported here, but the percentage of Goodall F values that are larger or equal to the observed ones (after 1000 permutations) falls around 10% for all populations.
Spatial structure and phylogenetic inertia
Evidence for spatial (i.e., geographical) patterns was not found in either pronotum centroid size or pronotum shape. Tests applied looking for phylogenetic inertia were based on the assumption of no-dispersal since the last inter-island connection that allowed faunal transfer. Under such an assumption, evidence for phylogenetic inertia was not found for either pronotum centroid size or pronotum shape (see Fig. 5).
Fig. 5. Correlograms of inter-population variations of A. planipennis. Fig. 5a refers to the case of pronotum size. Moran’s I statistic is used for testing both the relationship between size and phylogenetic inertia, and between size and geographical distance. Fig. 5b refers to the case of pronotum shape. Mantel’s Z statistic is used for testing both the relationship between shape and phylogenetic inertia, and between shape and geographical distance). Three distance classes were defined. Class marks from the pair-wise bathymetric matrix (used as estimate of phylogenetic distance) wee 0 to 30 m, 30 to 60 m and more than 60 m. Class marks from the pair-wise geographic distances were 0 to 20 km, 25 to 80 km and more than 80 km. Evidences for spatial structure or phylogenetic inertia were not found because none of the Moran’s I values were significantly larger (i.e., close to 1.0) than those expected to occur by chance. Evidences for spatial structure or phylogenetic inertia were not found because none of the Mantel’s Z values were significantly smaller (i.e., close to -1.0) than those expected to occur by chance.
This general assessment of the lack of evidences for either spatial patterns or phylogenetic inertia should be commented in detail. Firstly, sample size (both, number of populations and number of specimens of certain populations) is small. This fact could reduce the statistical power of the correlograms. However, it should be noted that there are some pairs of geographically close populations showing clear differences in shape and size (e.g., Illa dels Conills Island and Na Moltona Island). Similarly, some pairs of phylogenetically distant populations show similar size and shape (e.g., Dragonera Island and Menorca Island). All these cases suggest that it is probable that the correlograms remain qualitatively unchanged by a mere increase of the number of sampled populations. The case of size and phylogenetic inertia is especially clear because the existence of autocorrelation would imply to invert the signal currently observed for the phylogenetically closer populations.
Secondly, including or excluding the population of A. moraguesi implies the same qualitative results. In addition and for the case of phylogenetic inertia, comparable results were obtained considering this species within the same three classes or including it within an ad hoc fourth class. The exclusion of A. moraguesi from the calculations relative to the phylogenetic inertia is justified because the alternative hypothesis of an ancient (i.e., anterior to the last Ice Age) speciation between this specie and A. planipennis should be also considered. In that case bathymetry could not be considered a reliable estimate of phylogenetic distance between A. planipennis and A. moraguesi.
The population quartered at Artà (north-eastern Mallorca; Fig. 1) is considered to be a separate species (A. moraguesi; Español, 1954; Pons & Palmer, 1996). This taxonomic status is supported by the results presented here. A. moraguesi shows both specific shape differences and a size-shape relationship far from the allometric trajectory displayed by A. planipennis. I would like to emphasise the usefulness of Thin-Plate Spline Analysis combined with Canonical Variates Analysis for intuitive locating, understanding and describing the key shape differences that discriminate A. planipennis and A. moraguesi (Fig. 3).
Focusing now on A. planipennis, there are no-significant effects of either spatial structure or phylogenetic inertia. Therefore, the alternative hypothesis (i.e., shape and size have evolved by selection related to site-specific environmental factors) should be considered (see also the next paragraph). In addition, size and shape seem to be significantly correlated at inter-population scale, but appear to be independent at intra-population level. Size and shape remain uncorrelated at intra-population level probably due to the small morphological variability occurring at such level (site-specific selective forces could constraint size and shape). However, size and shape are correlated at inter-population level. This correlation became observable only at this scale probably because of the large environmental changes linked to shifting from islands ranging from 0.06 to more than 1000 km2.
The relationship between body size and some environmental variables has been recently analysed for the same populations of A. planipennis considered here (Palmer, submitted). The relationship between island size (as a global variable that summarises some environmental features) and body size delineates an interesting bell-shaped pattern (Palmer, submitted). Populations from island with an specific area would display the largest body size, and departures in both directions from such an island size (i.e., larger or smaller islands) would imply body size reduction. This is an implicit prediction for isolated populations of the same species into a more general theoretical framework known as the island rule. The island rule (Foster, 1964) initially claimed that different groups of flightless non-marine mammals experience different trends in body size changes (i.e., carnivores tend to dwarfism and rodents to gigantism) but has been generally reformulated as a graded trend from gigantism in small-bodied species and to dwarfism in large-bodied species (Brown & Lomolino, 1998). Lomolino (1985) proposed that ecological release (e.g., decreased competition or predation) would promote increasing body size and, conversely, resource limitation would decrease it. In the case of A. planipennis, body size increases from tiny islands to medium-sized islands, reaching a theoretical maximum at 38 km2, and additional increases of island area are linked to body size decreases (Palmer, submitted).
Alcover JA, Sans A, Palmer M. 1998. The extend of extinctions of mammals on islands. Journal of Biogeography 25: 913-918.
Bookstein FL. 1991. Morphometric tools for landmark data: Geometry and Biology. Cambridge: Cambridge Univ. Press.
Brown JH, Lomolino MV. 1998. Biogeography. Sunderland: Sinauer Assoc.
Cuerda J. 1975. Los Tiempos Cuaternarios. Palma de Mallorca: Institut d’Estudis Balearics.
Diniz-Filho JAF, de Sant’Ana CER, Bini LM. 1998. An eigenvector method for estimating phylogenetic inertia. Evolution 52: 1247-1262.
Español F. 1954. Los Tenebriónidos (Col.) de Baleares. Trabajos del Museo de Ciencias Naturales de Barcelona (Nueva serie Zoológica) 1: 1-93.
Foster J. 1964. The evolution of mammals on islands. Nature 202: 234-235.
Gittleman JL, Kot M. 1990. Adaptation: statistic and null model for estimating phylogenetic effects. Systematic Zoology 39: 227-241.
Klingenberg CK, McIntyre GS. 1998. Geometric morphometrics of developmental instability: Analyzing patterns of fluctuating asymetry with procrustes methods. Evolution 52: 1363-1375.
Legendre P, Legendre L. 1998. Numerical Ecology. Amsterdam: Elsevier Science B.V.
Lomolino MV. 1985. Body size of mammals on islands: The island rule reexamined. The American Naturalist 125: 310-316.
Marcus LF. 1993. Some aspects of multivariate statistics for morphometrics. In: Marcus LF, Bello E, Garcia-Valdecasas A,eds. Contributions to morphometrics. Madrid: CSIC-Monografias del Museo Nacional de Ciencias Naturales, 96-130.
Mayol J. 1985. Rèptils i amfibis de les Balears. Palma de Mallorca: Ed. Moll.
Palmer M, Pons GX, Cambefort Y, Alcover JA. 1999. Historical processes and environmental factors as determinants of inter-island differences in endemic faunas: the case of the Balearic Islands. Journal of Biogeography 26: 813-823.
Pons GX, Palmer M. 1996. Fauna endèmica de les Illes Balears. Palma de Mallorca: Inst. Est. Balearics,Conselleria d’Obres Públiques, Ordenació del Teritori i Medi Ambient, Soc. Hist. Nat. Balears.
Rohlf FJ. 1993. NTSYS-pc. Numerical taxonomy and multivarite analysis system, version 1.18. New York: Exeter: Setauket.
Rohlf FJ. 1998a. TPSRegr: a program for regresion of partial warps scores. New York: Dept. Evolutionary Biology. Univ. New York. Stony Brook.
Rohlf FJ. 1998b. TPSRelW.Thin-plate spline relative warp . New York.: Dept. Evolutionary Biology. State Univ. of New York. Stony Brook.
Rohlf FJ, Marcus LF. 1993. A revolution in morphometrics. Trends in Ecology and Evolution 8: 129-133.
Salvador A. 1979. Materiales para una “Herpetofauna Balearica”. Taxonomia de las lagartijas baleares del archipiélago de Cabrera. Zool. Beitr. Bonn 30: 176-191.
Slice DE, Bookstein FL, Marcus LF, Rohlf FJ. 1996. A glossary for Geometric Morphometrics. In: Marcus LF, Corti M, Loy A, Naylor GJP, Slice DE,eds. Advances in Morphometrics. New York: Plenum Press, 531-551.
Tuccimei P, Ginés J, Delirtala C, Pazzelli L, Taddeucci A, Clamor B, Fornós J, Ginés A, Gràcia F. 2000. Dataciones Th/U de espeleotemas freáticos recolectados a cotas inferiores al actual nivel marino en cuevas costeras de Mallorca (España): aportaciones a la construcción de un acurva eustática detallada en los últimos 300 ka para el Mediterráneo Occidental. Endins 23: 59-71.
I thank Les Marcus, Tonyo Alcover, Txus Gómez-Zurita, Carlos Juan, and José Aleandre Felizola Diniz-Filho for their help and suggestions. I also thank James F. Rohlf and all other people that make morphometric software freely available on the web. This research was granted by the scientific project DGES PB96-0090, Ministry of Education and Sciences (Spain). This paper is an hommage to the late Les Marcus.