Geometric morphometrics in entomology

Entomological Science (2018) 21, 164–184

doi: 10.1111/ens.12293

INVITED REVIEW

Geometric morphometrics in entomology: Basics and applications Haruki TATSUTA1 , Kazuo H. TAKAHASHI2 and Yositaka SAKAMAKI3 1

Department of Ecology and Environmental Sciences, Faculty of Agriculture, University of the Ryukyus, Nishihara, Japan, Graduate School of Environmental and Life Science, Okayama University, Okayama, Japan and 3Entomological Laboratory, Faculty of Agriculture, Kagoshima University, Kagoshima, Japan 2

Abstract The recent expansion of a variety of morphometric tools has brought about a revolution in the comparison of morphology in the context of the size and shape in various fields including entomology. First, an overview of the theoretical issues of geometric morphometrics is presented with a caution about the usage of traditional morphometric measurements. Second, focus is then placed on two broad approaches as tools for geometric morphometrics; that is, the landmark-based and the outline-based approaches. A brief outline of the two methodologies is provided with some important cautions. The increasing trend of entomological studies in using the procedures of geometric morphometrics is then summarized. Finally, information is provided on useful toolkits such as computer software as well as codes and packages of the R statistical software that could be used in geometric morphometrics. Key words: elliptical Fourier analysis, Procrustes analysis, thin-plate spline.

INTRODUCTION In entomology, much attention has been focused on the morphology of various body parts for identifying, naming and classifying organisms. Traditional taxonomy and systematics use the specific features of external morphology as diagnostic characters for the identification of taxa. In general, there are two categories that represent these morphological characters, “metric (or continuous) traits” and “meristic (or discrete) traits”. Morphometric comparisons are generally made using metric features of the body parts such as head and thorax. The traditional and most convenient way is to measure the linear distances between predetermined points on each body part and then compare these measurements using an appropriate statistical approach. Generally, in univariate or multivariate comparisons of morphological features, a linear distance-based approach is being followed; however, it provides less information on the morphological features. Furthermore, these metrics tend to be correlated with each other to a certain extent, and thus, may Correspondence: Haruki Tatsuta, Department of Ecology and Environmental Sciences, Faculty of Agriculture, University of the Ryukyus, Nishihara, Okinawa 903-0213, Japan. Email: [email protected] Received 14 July 2017; accepted 14 October 2017; first published 26 December 2017.

lead to considerable redundancy in the measurements (Zelditch et al. 2004). Hence, such traditional methods of measurement are not very effective for describing overall features of the body parts of interest. Another measurement approach, the use of ratios, also suffers from limitations. The ratio of independent standard normal random variables, such as the ratio of head length to the width, is sometimes regarded as a morphological feature. Such ratios are often considered as indicators of shape differences in the body parts of interest. However, such ratios vary with the sum of proportional and disproportional changes in the measured distances. Moreover, the mathematical property of ratios is quite ambiguous (Atchley et al. 1976). Furthermore, the ratio is sometimes compared under the a priori assumption of already known probability distributions (e.g. the Gaussian distribution), but such a comparison is inappropriate. It is known that the ratio of two independent standard normal random variables conforms to a Cauchy distribution, where the mean and higher moments do not exist (Mood et al. 1974), and thereby it is not recommended to compare the ratio of variables under the assumption of the conventional statistical distribution, such as the multivariate Gaussian distribution (Atchley et al. 1976; Atchley & Anderson 1978). Of course, the above condition postulates a rather unrealistic assumption, and in most of

© 2017 The Entomological Society of Japan

Geometric morphometrics in entomology

cases, the measurements are more or less correlated with each other. Under the relaxed conditions, however, it should be kept in mind that the ratios may induce undesirable statistical artifacts (Atchley & Anderson 1978); thus, better to avoid statistical comparisons without proper prudence. Penrose (1954) defined the total form as the summation of “size” and “shape”. In fact, there are numerous definitions in terms of size and shape (Rohlf & Bookstein 1987; Lestrel 2000) and these two notions are sometimes confounded in traditional morphometrics for some of the reasons outlined earlier. However, one can define size and shape in a geometrical manner. The definition of “shape” by Kendall (1977) is “all the geometrical information that remains when location, scale and rotational effects are filtered out from an object”. The word “scale” in Kendall’s definition of shape refers to “size”, which can be regarded as positive scalar a satisfying the equation g(aX) = ag(X), where g(X) represents any positive real-valued function of the configuration matrix (Mosimann & James 1979; Dryden & Mardia 1998). In other words, size change can be attributed to the proportional change when the two forms are compared. In traditional morphometrics, the size is modeled as a latent variable in a path diagram, named as “general size”, which co-varies intrinsically with shape. The approximate effect of size on the total measurements can be assessed using a path diagram (Humphries et al. 1981; Rohlf & Bookstein 1987; Bookstein 1991). Several methods for the adjustment of the size effect of measurements have also been developed (e.g. Burnaby 1966; Humphries et al. 1981; Rohlf & Bookstein 1987), but what one needs to consider here is which measure is suitable as a representative of general size. When a set of linear measurements of body parts is available, one sometimes regards the first principal component as general size because the loadings tend to be all positive and become similar with each other (e.g. Lande & Arnold 1983; Bolnick & Lau 2008). However, this conventional procedure has been shown to be sometimes quite erroneous, and thus, is inappropriate for size correction (for more details, see Berner 2011). Over the last few decades, the mathematical theory and geometric morphometric tools have been intensively developed. However, a major concern in geometric morphometrics has been how to analyze 2-D or 3-D Cartesian coordinates of a certain set of end-points across the specimen (Rohlf & Marcus 1993). The theory is rather involved for most researchers (and probably also for most entomologists!) who are not familiar with linear algebra; nevertheless, this approach definitely outweighs traditional morphometrics with reference to partitioning explicitly the geometric information into components of

Entomological Science (2018) 21, 164–184 © 2017 The Entomological Society of Japan

size and shape (Zelditch et al. 2004). To start, a size measure, “centroid size”, is introduced, which is calculated as the square root of the summed squared distances of each end-point from the centroid of a form (Zelditch et al. 2004). In contrast to general size, the measure of centroid size is uncorrelated with shape as long as shape varies isometrically (i.e. in the absence of allometry) and is used in the definition of the Procrustes distance (see subsequent section) (Bookstein 1991). In geometric morphometrics, the centroid size can be referred to as “size”. Also of importance is that the traditional approach loses the information of the coordinate positions (point locations) when geometric information is converted into linear distances; thus, geometric morphometrics enables one to draw intelligible pictures as the information of the end-point positions still remains after appropriate removal of irrelevant information, such as the position and orientation of the specimens. This would be particularly useful for biologists because concurrent local changes, in particular, regions of morphological structures, would then be easy to understand. In this review, a brief outline will be presented of the methods that are used within geometric morphometrics and several comments provided on technical issues regarding Procrustes superimposition and the use of elliptical Fourier (EF) analysis. Although both methods are feasible for comparing shape features after removing all information unrelated to shape, Procrustes superimposition relies on explicit landmarks and aims to minimize differences between configurations based on Procrustes distance (i.e. the summed squared distances between corresponding landmarks) (Zelditch et al. 2004). Elliptical Fourier analysis, in contrast, does not need explicit landmarks and was developed for fitting curves to complex closed contours (and even for open curves) (Kuhl & Giardina 1982). Elliptical Fourier analysis is superior to other Fourier methods in some respects (Rohlf & Archie 1984 and see later sections). In order to review the entomological studies where the geometric morphometric approach was applied during the past 15 years (2000–2015), focus is placed on those studies with respect to body parts of interest and the used methods of geometric analysis. In the final section, information is provided on some useful toolkits such as ready-to-use software and computer codes, and future perspectives, especially for 3-D data acquisition in entomology.

OVERVIEW OF GEOMETRIC MORPHOMETRICS Prior to the explanation of methods used in geometric morphometrics, one needs to explain the word

165

H. Tatsuta et al.

landmarks involved so as not to confound these different types of landmarks.

Landmark-based approaches

Figure 1 Male Dorcus titanus palawanicus (left) and magnified left mandible (right). Gray dots, landmarks; small white dots, semi-landmarks (pseudo-landmarks). Photograph courtesy of K. Goka.

“landmark”. Here the well-developed mandible of male in the stag beetle Dorcus titanus palawanicus (Lacroix, 1984) is used as an example. Five predefined end-points along the outline of left mandible (gray circles) are displayed in Figure 1. These are called “landmarks” and are defined as points of correspondence on each object that matches between and within populations (Dryden & Mardia 1998). These points consist in juxtaposition of tissues (type I), maxima of curvature (type II) and external points (type III) (Bookstein 1991). In addition to the above five points, one might also consider a few other points to be taken as landmarks (e.g. on the external points of small teeth along a row between a tip of mandible and the first inner tooth). Although it is possible to take at least some of these points into account for comparisons made on intraspecific levels (i.e. comparison between different local populations), it would sometimes happen that one could take only a few homologous landmarks as points with confidence when the shape of the mandible varies greatly among different species or higher taxa (e.g. jaggy inner teeth often disappear in some species) (see “Outline-based approaches” section below). These remaining landmarks, therefore, would not be sufficient to describe the fine-scale features of the mandible. In such cases, one may want to include an arbitrary number of points between landmarks as seen along a curvature between the first inner tooth and proximal point (Fig. 1, white circles indicated by arrow). These points are called “semi-landmarks” or “pseudo-landmarks”, where the point-to-point homology is not guaranteed among the different samples to be compared, and thus they cannot be considered identical with homologous landmarks as defined above. Therefore, one has to choose appropriate methods depending on the type of

166

It is not an exaggeration to say that landmark-based methods are at the center of geometric morphometrics because the definition of size and shape is reasonably clear and the mathematics now well established. Although there are plenty of choices when one examines the shape of biological configurations, it is essential to use the most appropriate method available with respect to the focus of the investigator’s research. Initially, the discussion will be on: (i) superimposition of landmarks; and (ii) the thin-plate spline method to familiarize readers with some of the special terms. However, the explanation will be brief here as one needs to consult the details of geometric morphometrics from some of the prominent textbooks that are available, e.g. the “Orange Book” by Bookstein (1991), the “Black Book” by Marcus et al. (1993), the “White Book” by Marcus et al. (1996), the “Blue Book” by Dryden and Mardia (1998, 2016), the “Green Book” by Zelditch et al. (2004), and more. It is especially encouraged that general readers study Zelditch et al. (2004) as this book provides comprehensive tutorials as well as various morphometric methods; however, some explanation on the statistics may be misleading (for details, see Rohlf 2005). The Procrustes superimposition method based on generalized least-squares (GLS) is first introduced for the comparison of configurations. The name “Procrustes” derives from a Greek mythical tavern owner who was a thief and murderer and tied his victims to a bed by either stretching or cutting their legs off to fit the bed. The GLS method has been extensively used in the biological sciences including entomology (e.g. Eberle et al. 2015; Laparie et al. 2016) and can be implemented using software and R library (e.g. tps series software, Rohlf 2015), or the shapes package in R library (Dryden & Mardia 2016), although it must be admitted that the underlying mathematical theory is complex. As indicated by Kendall (1977), the GLS method relies on translation, scaling and rotation to eliminate all information unrelated to shape (Zelditch et al. 2004). Prior to using the optimization procedure, an overview will be provided of the relationship of different “morphospaces” where the configurations are to be positioned and aligned. To understand these “spaces”, the following terminologies will be referred to: “figure space”, “pre-form space”, “pre-shape space”, “conformation space” and “shape space” as illustrated in



Figure 2 Relationships between various morphospaces. Dots in the center of triangles indicate centroid. Modified after Dryden and Mardia (2016).

Figure 2 (for the further details of these morphospaces, see Goodall 1991; Minaka 1999; Klingenberg 2016). Consider a k × m matrix of coordinates of an object representing k landmarks in m-dimensional space and another different form with the same points existing in the same dimensional space but arbitrarily positioned in that plane or space. Both of the objects constitute a “figure space” where both size and shape are still retained. After appropriate centering and translation, these figures are then arrayed into a “pre-form space”. There are two different steps that are available for further optimization of figures to achieve the resulting superimposition of the shapes (Fig. 2). One step is to standardize the size of figures first and then to rotate the figures to minimize the sum of distance between corresponding landmarks (= centroid size). The other step is to rotate figures first and then to standardize the size of figures. Either of the approaches achieves the concurrent figures shown in the spherical hyperspace (= “shape space” hereafter). After size standardization using centroid size for the two separate objects, the new configurations now meet the restrictions of the “pre-shape space” and each of the new configurations is a centered pre-shape which constitutes a km − m − 1 dimensional hypersphere (Dryden & Mardia 1998, 2016). For a set of triangles, the pre-shape space is 3-D, and therefore has to be visualized on a hypersphere. Although the hypersphere for quadrilaterals or polygons cannot be easily visualized, the configurations are arrayed in the same way. One then rotates the target figure so as to minimize the square root of the sum of squared distances between corresponding landmarks of the reference and target figure. After the treatment, the smallest great circle ρ and chordal (Euclidian) distance Dp (Fig. 3) between reference and target figure can be consequently found;


the latter corresponds to the partial Procrustes distance. This chordal distance Dp may not be the shortest possible “distance” between the two shapes, and then, one can find the shortest distance by changing the constraint of the centroid size of one shape (= target shape) so as to minimize further the distance from the other shape (= reference shape). The shortest point can be achieved with a line that is perpendicular to the target’s radius and passes through the reference’s position on the surface (Fig. 3). The corresponding centroid size of the target becomes cos(ρ) and the distance from the reference to the now re-aligned configuration on a new shape space is sin(ρ), that is, the full Procrustes distance Df. This “rescaling” produces a new shape-space hypersphere with a radius of 1/2, tangential to the previous shape space at the reference shape (Fig. 3). This rescaling process is equivalent to a Riemannian submersion, where a configuration on the tangent space to the pre-shape sphere is isometrically transformed onto the tangent space to the shape space (for further details, see Kendall et al. 1999); that is, the shape itself is invariant after the transformation. This new space is specifically called “Kendall’s shape space”, where the dimension is decreased after the optimization of rotation is made, resulting in km − m − 1 − m(m − 1)/2 dimensions. The other optimization for obtaining figures in shape space, rotation of figures, precedes the size standardization. After the rotation, the figures are arrayed on “conformation space” of km − m − m(m − 1)/2 dimensions (Rohlf 1996; Klingenberg 2016; but see Lele 1993 for a different number of dimensions). The same space has been called equivalently as “form space” (Rohlf 1996) or “size-and-shape space” (Dryden & Mardia 1998, 2016) or “allometric space” (Langlade et al. 2005), but Klingenberg (2016) proposed the term

167

H. Tatsuta et al.

Figure 3 Photograph of a section of the tangent space, aligned pre-shape space and Kendall’s shape space for the triangle case, after Zelditch et al. (2004) and Claude (2008). O corresponds to the mean reference shape, and Tp and Tf are preshape and Kendall position of a shape in the sphere, respectively. Dp and Df indicate partial Procrustes distance and full Procrustes distance, respectively.

“conformation” in order to distinguish the terms according to these intrinsic meanings. It should be noted that the conformation space is not a sphere, but a cone with a warped-product metric and the hemisphere of Procrustes-aligned pre-shapes is a special subspace in the conformation space because pre-shapes are defined as configurations that have been scaled to centroid size 1.0. Conformations with centroid sizes greater than 1.0 exist outside the hemisphere and smaller configurations fall into inside the sphere (for further geometric characteristics on this space, see Kendall et al. 1999; Klingenberg 2016). Because all of the shape spaces defined above are curved and thus non-linear, one cannot adopt statistical procedures that postulate a Euclidian space for comparing shapes (note that the most precise comparisons of aligned shapes can be conducted under non-linear space without any projection onto a Euclidian subspace). The surface of the hypersphere can be projected onto a tangent space where the data can be analyzed in Euclidian space (Fig. 3). The orthogonal or stereographic projection onto the tangent space might induce some non-trivial biases, but this would not happen in most of the cases in a real biological dataset (Klingenberg 2016). The GLS procedure mentioned above enables one to compare shapes on this space, and thus various multivariate statistical methods, such as principal component analysis and canonical discriminant analysis, can be adopted under the assumption of linearity. Here is the summary of the GLS fitting procedure (after Rohlf 1990a): 1. Center each configuration of landmarks at the origin by subtracting the coordinates of its centroid from

168

the corresponding coordinates of each landmark; it enables translation of each centroid to the origin. 2. Scale the landmark configurations to unit centroid size by dividing each coordinate of each landmark by the centroid size for that configuration. 3. Choose one configuration as a reference (base), and rotate the second configuration (target) to minimize the summed squared distances between corresponding landmarks (for all the landmarks) between the forms. If the number of forms is three or more, the reference shape is the average shape of different forms, which can be obtained by an iterative fitting procedure (Zelditch et al. 2004; Dryden & Mardia 2016). This procedure is called Procrustes analysis (Dryden & Mardia 1998, 2016) and the overall procedure is illustrated in box 2 in Klingenberg (2010). The prefix term partial corresponds to the optimization suing rotation and full includes rescaling after rotation for the optimization (see also the above explanation on shape space and Fig. 2). The bewildering terms “ordinary” and “generalized” in Procrustes analysis also need to be distinguished; the former refers to the case where there are only two configurations to be matched, whereas the latter involves a number of configurations greater than two. For most cases, the generalized Procrustes analysis is preferred, especially for estimating the “average” shape from the samples. The second approach to be discussed here is the thinplate spline. The spline is particularly useful for visualizing the localized deformation of a target shape from a reference shape by using deformation grids as first presented by D’Arcy Wentworth Thompson (Thompson 1917). Recall that the effects of translation along



Figure 4 Alternative models of uniform (affine) shape deformation: (a) translation along the vertical axis; (b) translation along the horizontal axis; (c) scaling; (d) rotation; (e) compression/dilation; (f ) shearing. The original (reference) square is shown with solid lines and the deformed shape is indicated by dotted lines (after Zelditch et al. 2004).

vertical and horizontal axes, rotation and scaling prior to the comparisons of shapes can be removed using the GLS procedure. According to Zelditch et al. (2004), these adjustments can be defined as an “implicit uniform (or affine) transformation”. In addition to these adjustments, compression/dilation and shearing as additional uniform transformations, which are called “explicit uniform deformations” and distinguished from implicit ones with respect to the accompanying shape change, also need to be considered; thus, in total, there are six uniform transformations (Fig. 4). The overall deformation from a reference to a target shape can then be expressed using these six uniform components and non-uniform (or non-affine) components. The name thin-plate spline represents a physical analogy involving the bending of a thin plate of metal. An interpolation function can be used for visualizing the transformation (termed as a “warping”) from one object (= reference) into another (= target) as a Cartesian deformation (Bookstein 1991). The function is formulated as zðx, yÞ = − UðrÞ = − r2 logr2 where r is the distance between a pair of landmarks in the reference configuration (scaled unit centroid size). The function U(r) satisfies the equation ∂2 ∂2 Δ U= + ∂x2 ∂y2 2

!2 U / δð0, 0Þ

where δ(0, 0) is delta function that is zero everywhere except at the origin (x = y = 0) but has an integral equal


to 1. U is the fundamental solution of the biharmonic equation Δ2U = 0, the equation for the shape of a thin steel plate raised as a function z(x, y) above the xy plane (Bookstein 1991). When the thin-plate spline deformation function f for two dimensions is considered, one needs to describe the deformations in the x and y directions separately so as to fit the reference configuration into the target configuration, i.e. fx(x, y) and fy(x, y), both of which are decomposed into uniform and non-uniform components (this framework can be easily extended to 3-D data, see Rohlf 1996). Here, a (k + 3) × 2 matrix Xtarget is defined to describe the coordination of the target configuration and k denotes the number of landmarks as defined above. Then the coordinates of landmarks in target configuration can be expressed as the product of (k + 3) × (k + 3) matrix L consisting of the set of the deformation functions f described above and (k + 3) × 2 weight matrix W such as Xtarget = LW: To solve for W, one obtains W = L − 1 Xtarget : The k × k upper left submatrix of L−1 is called the bending energy matrix L−1k (Bookstein 1991). The principal warps are a set of eigenvectors that span the tangent space defined by the non-uniform (or non-affine) components of the thin-plate spline and the first k − 3 column vectors of k × k eigenvectors estimated from the singular value decomposition of L−1k, and the principal warps are thus derived from the reference configuration.

169

H. Tatsuta et al.

The resulting eigenvalues correspond to bending energies that can be regarded as an inverse measure of scale; note that large eigenvalues correspond to eigenvectors representing small-scale features. The inner product of each row of V = Xtarget − Xreference (Xreference: coordinates of reference configuration) with the k − 3 columns of the principal warp eigenvectors yields partial warp scores (Rohlf 1996). The partial warp scores are defined from the optimization of the bending energy required to deform the grid, which is a function of the relative distance among the landmarks in the reference configuration. Therefore, analyzing the distribution of objects in the bivariate space defined by each partial warp is not recommended because a portion of partial warps describes a portion of the total shape change; the array of partial warp scores are determined by the magnitude of bending energy and not by a certain biological background (Bookstein 1996). Therefore, all the partial warp scores should be used together with two uniform components such as dilation/compression and shearing (Fig. 4) for the comparison of shape features; other uniform components do not differ because they have been already eliminated. It should also be noted that the coefficients of splines could be used to perform statistical analyses in the Euclidian space tangent to Kendall’s shape space (Bookstein 1991; Loy 2007). The relative warp analysis refers to the principal component analysis (PCA) of the partial warp scores for a set of n specimens when a weighting scalar α is set to 0 (Loy 2007). The value of α can be changed to other than 0, making the output a weighted PCA. The principal component vectors of the newly integrated matrix consisting of the partial warp scores of each specimen are called relative warp scores. They can be visualized as a transformation grid that shows deformations of the physical space of the reference configuration (Rohlf 1996). These calculations can be carried out using several software packages that are available as well as the R libraries, e.g. tpsRelw by Rohlf (2015) and shapes by Dryden and Mardia (2016). Understanding the fundamentals of shape spaces and splines is also useful for studying more advanced analyses such as morphological integration and functional aspects of evolution (Klingenberg 2010; Suzuki 2013). These studies include morphological integration and modularity (note that “modularity” is used in various contexts in ecology, but in geometric morphometrics and evolutionary developmental biology it refers to the pattern of covariation among traits). For readers who are interested in these crossover fields, it is recommended to read some good reviews (e.g. Klingenberg 2008, 2010), plus recent enthusiastic papers.

170

Effect of localized shape differences on shape analysis As discussed above, superimposition methods are generally based on the superimposition of one morphology onto another using landmark positions. Among them, GLS is the most frequently used method. The GLS method averages and allocates the change in shape to all the landmarks rather than to the ones with large displacements. Therefore, the GLS method is most appropriate when shape change is globally distributed and not localized to a few landmarks. However, use of the GLS method may be problematic if there is the presence of a Pinocchio effect (see below). Consequently, the question may be raised as to what happens if the GLS method is applied to a shape comparison where the shape change is strongly localized to a few landmarks. A hypothetical example of such a highly localized shape change in Drosophila wings is shown in Figure 5. In this example, the landmark-based approach is applied to quantify wing shape using eight homologous landmarks on the intersection of the wing veins (Fig. 5a). The shape difference between the two configurations (Fig. 5b,c) is localized to the two landmarks at the proximal side of the wings (Fig. 6a, arrows). The effect of such an extremely localized deformation on superimposition is called a “Pinocchio effect” (Zelditch et al. 2004). When the two configurations are superimposed using GLS, the Procrustes distance is then used as the criterion for optimal superimposition and the sum of the squared deviations of all the corresponding landmarks is minimized. As a result, the localized displacement of the two landmarks is averaged among all the landmarks (Fig. 6b) and the displacement of the two landmarks at the proximal position becomes understated (Fig. 6b). Generally, the Pinocchio effect is not a serious problem when the purpose of shape quantification is to describe the morphological differences between groups (species, populations, and communities etc.) or to test whether or not there exists a difference in shape between the groups, because geometric morphometric analyses tend to work appropriately regardless of the presence of such effects. However, in some cases, such as morphogenesis studies, where the localization of the shape deformation itself is of interest, the Pinocchio effect may be an issue. One of the factors governing the morphogenesis is morphogens, which are secreted signaling molecules that organize the spatial pattern of surrounding cells. The release of a morphogen from a source forms a morphogen gradient, which leads to the correct spatial arrangement of different cell types (Gurdon & Bourillot 2001). These morphogen gradients imply that morphogenesis is regulated by a specific



Figure 5 (a) Homologous landmarks on a wing of Drosophila melanogaster; (b,c) shape deformation between two hypothetical wings.

directionality and at a specific spatial scale, suggesting that the effect of morphogenic genes on development and morphogenesis is spatially limited and may result in local deformations. The wings of Drosophila melanogaster is one of the best-studied model systems in insects for elucidating the genetic basis of development and morphogenesis (Held 2002). During wing development, a set of cells is predestined at an embryonic stage to form a wing and wing axes are broadly established at the early larval


stage (Cohen et al. 1991). In this process, the posterior region of the wing imaginal disc is patterned by the protein Engrailed that activates the short range paracrine signaling ligand Hedgehog at the boundary between the anterior and posterior regions of the disc (Garcia-Bellido & Santamaria 1972; Lowrence & Morata 1976; Brower 1986; Hidalgo 1994; Tabata & Kornberg 1994; Sanicola et al. 1995). Hedgehog induces the expression of Decapentaplegic (Dpp), the canonical ligand of the transforming growth factor-β

171

H. Tatsuta et al.

signaling pathway that helps to establish the anterior– posterior axis of the wing (Zecca et al. 1995). Dworkin and Gibson (2006) tested the effect of 50 insertional mutations in several signal transduction pathways, including the transforming growth factor-β signaling pathway, on quantitative wing shape variation with the landmark-based geometric morphometric approach. They found that the majority of those mutations had a profound effect on wing shape, but the interpretation of the shape deformation might not be clear in biological terms because the landmark configurations were superimposed with GLS. The candidate genes were expected to have a spatially limited expression pattern and may be responsible for extreme localized deformations. Therefore, deformation in geometric shape may not fully reflect biological shape deformation in such cases. The interpretation of geometric shape deformation in a biologically relevant way is a common challenge in studies of genetic and developmental processes leading to wing shape in D. melanogaster (Palsson & Gibson 2004; Takahashi et al. 2011a). One approach in dealing with localized shape deformation is to use resistant-fit superimposition. A distinct difference between GLS and the resistant-fit superimposition methods is that the latter does not use Procrustes distance as the criterion for the optimization of the superimposition (Zelditch et al. 2004). The oldest and most well-known resistant-fit superimposition method is resistant-fit theta-rho analysis (RFTRA) which utilizes “repeated medians” to determine the optimal scaling and rotation for the superimposition (Zelditch et al. 2004). The RFTRA is robust or insensitive to a small number of landmarks with relatively large displacements because the local deformation does not affect the median scaling factor or the median rotation angles used for the superimposition (for details of RFTRA, see Siegel & Benson 1982). Figure 6(c) illustrates the result of the superimposition of the two configurations in Figure 5(b,c) using RFTRA. Compared to the superimposition using GLS (Fig. 6b), the RFTRA superimposition addresses the Pinocchio effect better than expected and does not induce the covariance of the landmark displacements that are clearly visible in the GLS superimposition. Nevertheless, there are several alternative resistant-fit methods based on specific distance matrices and they can produce different superimpositions, resulting in possibly different biological inferences (Zelditch et al. 2004). Although resistant-fit superimposition methods look promising, especially when dealing with extremely localized deformations, they have some disadvantages as well. As the resistant-fit superimposition does not use the Procrustes distance metric, it is not consistent

172

Figure 6 Comparison of superimposition methods. (a) Simple superimposition with arrows pointing at two landmarks with large displacement; (b) generalized least-squares superimposition; (c) resistant-fit theta-rho analysis superimposition.

with the general theory of shape (Zelditch et al. 2004). In contrast, GLS uses Procrustes distance metric for superimposition, which does not alter the shape as defined by Kendall (1984). When GLS is applied, the shape distance or geometrical similarity of two configurations can be expressed as a distance between points in a mathematically well-defined shape space (a similar drawback arises in outline analysis: see the next section). In contrast, the similarity of shape can also be evaluated by performing multivariate analyses on the coordinates obtained from the resistant-fit superimposition. Given the advantages and disadvantages of resistant-fit superimposition methods, based on the above analysis it is suggested that the GLS approach is generally recommended unless extremely localized deformations are present and biological inferences are more important than geometric ones.



Outline-based approaches Researchers often compare curves that lack explicit landmarks, such as short segments between edges or ridges between landmarks or line segments along body parts to be of interest. In such cases, one of the alternative options is to use an outline-based approach such as Fourier analysis (Rohlf & Archie 1984) or eigenshape analysis (Lohman & Schweitzer 1990). In these methods, particular functions such as the Fourier series or tangent angles are fitted to the digitized coordinates. Then a set of coefficients representing features of the curves are used as variables in subsequent analyses. However, these coefficients do not correspond to the definition of shape by Kendall (1977), and thus, cannot be easily compared with shape descriptions based on configurations of landmarks. In contrast to the methods based on landmarks, homologous landmarks are not needed in outline-based approaches. This would be greatly advantageous to compare curvilinear features in which even a few landmarks are very difficult to locate. Instead of explicit landmarks, a number of points can be put, as semilandmarks (or pseudo-landmarks), on the outline of interest. However, it should be noted that these semilandmarks just provide information on the “bowing” of the curve (i.e. the degree of deviation from a straight line connecting two explicit landmarks). Although several alternative approaches such as differential weighting and sliding semi-landmarks have been proposed, they still contain unsolved problems with respect to shape space and degrees of freedom (see chapter 15 in Zelditch et al. 2004 for this topic). The non-landmark-based approach to describe the curve is to fit appropriate descriptors to it. Although this procedure is placed outside of the rigid mathematical framework of shape space and is sometimes criticized by the members of the “landmark school”, the outline-based approach apparently outweighs the conventional landmark-based approach in approximation of concave and convex segments (e.g. MacLeod 1999). For this, different methods have been developed for two classes of outlines: open contour curves and closed contour ones. To estimate open contour curves, a simple polynomial equation, cubic splines (de Boor 1978; Evans et al. 1985) or Bezier curves (Engel 1986) can be adopted. In addition, Fourier analysis (described later) has been devised in an open curve approach (Lestrel & Wolfe 2017). The most appropriate method should be chosen depending on the complexity of the curves in question. For comparing the coefficients between two specimens, it is necessary to rotate and scale the objects


appropriately in such a way as placing the x-axis through the two reference points (e.g. starting and ending points of the curves), which are analogous to the Bookstein shape coordinates (Bookstein 1991). For closed contour curves, Fourier series decomposition with discrete Fourier transform is widely used for fitting the observed points along the outline (Lu 1965; Rohlf & Archie 1984; Rohlf 1990b). Given data consisting of the lengths of equally or unequally spaced radii emanating from a center of the biological form, one simple solution for fitting points is to adopt appropriate harmonic functions such as the Fourier function (Kuhl & Giardina 1982). The least-square estimates of the Fourier coefficients are given in Rohlf (1990b). However, it should be noted that the traditional Fourier transform cannot be applied to rectilinear coordinates as it stands, which is the default output of some software routines used for data acquisition (e.g. tpsDIG by F. J. Rohlf, also see later section in this review). In this case, the data must be transformed to polar coordinates relative to the centroid of the contour (Rohlf 1990b). In order to compare the Fourier coefficients of one specimen with those of others, one must make sure that the angle of the starting radius and digitizing direction (clockwise or counterclockwise) are the same for all specimens to be digitized. If this requirement is not satisfied, then the outlines need to be oriented so that the intersection of the first radius with the outline corresponds to the starting point of the measurements. Alternatively, if the digitizing direction is clockwise (instead of counterclockwise), one can remedy this easily by taking a mirror reverse of the raw data. The Fourier analysis of equally or unequally spaced radii needs a larger number of points closer to the central point (often the centroid) within the object. The presence of unequally spaced, sparse points may produce a bias that will need to be interpolated with more radii as the shape of the object increasingly deviates from a circular form (Rohlf & Archie 1984). In this case, another Fourier analysis can be adopted for contours with high complexity. This method expresses the contour by the equation ϕðt Þ* = ϕðt Þ− ϕð0Þ− t where t is the cumulative chordal distance along the curve, ϕ(t) is the angle of the tangent vector at the distance t, and ϕ(0) is the angle of the tangent vector for the starting point, which is omitted from the standardized data (Zahn & Roskies 1972). The standardization of t is required for comparing different specimens. Again, the starting point of measuring angles and the

173

H. Tatsuta et al.

Figure 7 Transformation series of reconstructed contours given different numbers (3, 5, 10, 15 and 30) of Fourier harmonics, shown in comparison with the original outline (rightmost) based on a picture of the Prosopocoilus inclinatus male mandible (left-most).

direction of taking points must be the same among specimens. Fourier coefficients are then calculated for the array of ϕ(t)* by the least-squares. The shape function of Zahn and Roskies (1972) has been extensively applied in eigenshape analysis (Lohman 1983; Lohman & Schweitzer 1990; MacLeod 1999), but a drawback in this method is the poor convergence of the tangent angle Fourier descriptors (see Rohlf & Archie 1984; Claude 2008). This is due to the “jagged” tangent angles of a contour in a digital image. Elliptical Fourier descriptors (see the next section) can avoid this problem because the 2-D coordinate system is relatively robust against the “jagged” curves (see the explanation of curve matching in Persoon & Fu 1986). Another drawback is that the reconstructed curve sometimes fails to close (for an example, see chapter 2 in Lestrel 1997). These drawbacks may be avoided by contriving to take semi-landmarks. Probably the most prominent method for the approximation of complex closed contours is the EF method (Kuhl & Giardina 1982; Ferson et al. 1985). The name is derived from the fact that each harmonic in the series is depicted by an ellipse with the first harmonic displaying the largest ellipse. The reconstructed contour approaches the original configuration as the number of harmonics increases (Fig. 7). The question is how to determine the optimal number of harmonics. There is no objective criterion for doing this, but one can infer the number of harmonics by computing the deviation (or residual) between the reconstructed and original outlines (for details, see Lestrel 1997; Claude 2008; Tatsuta et al. 2013). However, the number of harmonics should not exceed half the number of points, as required by the Nyquist frequency criterion (Lestrel 1997; Tatsuta et al. 2013). The EF method does not require that the data points be equally spaced on the boundary outline as other Fourier analytic approaches postulate (including the Fourier transform). Furthermore, the procedure proposed by Kuhl and Giardina (1982) for calculating Fourier descriptors can be applied to any type of descriptors based on discrete Fourier transformation. This method will fit outlines of

174

any complexity, even those with intersections and crossings, given that enough harmonics are taken (Rohlf & Archie 1984). Although empirical studies using the EF function are overwhelming for the 2-D outline, such studies can be easily extended to 3-D curves (Lestrel 1997; Claude 2008). Furthermore, this approach can also be adopted for open curves after an appropriate conversion of raw data, e.g. by taking a mirror image of the original configuration and combining it with the original so as to form a closed contour (see Lestrel & Wolfe 2017). Because the number of coefficients tends to be so large, the coefficients are usually summarized by ordination methods such as principal component analysis or canonical variate analysis for multiple groups to be compared. In such cases, one has to make sure that the statistical properties of Fourier coefficients satisfy a prerequisite of analysis to be adopted. For example, linear discriminant analysis usually requires multivariate normality, as demanded in ANOVA or MANOVA, and the conclusion based on the data that violate the assumption might be erroneous. However, these problems can be avoided by resampling procedures such as bootstrapping or cross-validation (see Zelditch et al. 2004). This issue will be addressed in more detail in the next section. Recently, the discrete wavelet transform has been applied to the form of biological outlines and has been reported to be superior to the EF method for identifying localized features. The reason for this is due to the property of periodicity; that is, while Fourier transforms range over the period −π to +π for the continuous case, wavelets are of limited duration and rapidly converge to zero at both ends of the interval (Lestrel 2000). Wavelet analysis can be done relatively easily in a ShapeR package in the R library (Libungan & Pálsson 2015).

Empirical notes for statistical use of EF coefficients After the calculation of numerous EF coefficients, one can choose either of the following two ways for



Figure 8 (a) Frequency distribution of correlation coefficients (r) between 77 standardized elliptical Fourier descriptors and (b) a part of the correlation matrix. The correlations of filled squares are significant after Bonferroni correction (P < 0.01): gray, 0.33 ≤ r < 0.8; black, 0.8 ≤ r.

statistical treatments: (i) directly treating the EF coefficients as variables for some multivariate analyses; or (ii) summarizing the EF coefficients by PCA first and then treating principal component scores as variables for some statistical analyses. In some cases, conventional multivariate statistical analyses (MANOVA, canonical variate analysis, discriminant analysis, etc) are implemented for numerous Fourier coefficients (e.g. Medal et al. 2003; Granier et al. 2005). Some of these statistical analyses require assumptions of normality, homogeneity of variance and independency of the variables utilized. For EF coefficients, however, there is no guarantee that they conform to normality, independency and homogeneity of variance. For example, the variance of EF coefficients in higher frequency tends to be smaller than that in lower frequency when the contours studied are smooth. This is obvious from the fact that the smoother the function is, the smaller the contribution of the highfrequency components in the Fourier series (Beerends


et al. 2003). In this case, the homogeneity of variance will be violated when the estimated EF coefficients are used as variables without appropriate transformation. For the independency of Fourier coefficients in mathematics, the Fourier coefficients are calculated from orthogonal axes so that they are independent of each other (Jacobshagen 1997); however, under the practical situation, some are found to be correlated, as seen in Jacobshagen (1997). Similar significant correlations among coefficients were found in the data of Fukudome and Sakamaki (2011) and Iwata and Ukai (2002). Fukudome and Sakamaki (2011) examined the variation of female genitalia shapes in three Trichogramma species (egg parasitoid wasps) by an EF analysis using the first 20 harmonics, and estimated 77 standardized Fourier coefficients (for the procedure of standardization, see Kuhl & Giardina 1982; Rohlf & Archie 1984). The correlation coefficients between these 77 standardized EF descriptors were mostly insignificant (Fig. 8a), but some significant correlations were

175

H. Tatsuta et al.

found between neighboring positions (Fig. 8b). Thus, the independency of Fourier coefficients was not guaranteed in this case. As reported in Rohlf and Archie (1984), the degree of contributions of Fourier harmonics usually varies and an individual Fourier coefficient is almost meaningless as a morphological descriptor. Therefore, shape variation is usually described as a set of Fourier coefficients. When this approach is taken in studies for biologically meaningful shape variation, one should keep in mind that some Fourier coefficients might have non-trivial correlations. Here, using the same example data of Fukudome and Sakamaki (2011), we examine whether the frequency distribution of EF coefficients conforms to normality. The Shapiro–Wilk W test (Shapiro & Wilk 1965) showed that 19 (24.6%) of the 77 coefficients deviated significantly (P < 0.05) from the normal distribution. For these coefficients, there was no particular trend with respect to skewness and kurtosis. In the case of another example dataset of Japanese radish (Iwata & Ukai 2002), 23 (29.9%) of 77 standardized Fourier coefficients deviated significantly from the normal distribution, without any particular tendency for skewness or kurtosis. Thus, in reality, our examples showed that 20–30% of EF coefficients deviated significantly from the normality. Consequently, when EF coefficients are used as variables in conventional statistical analyses that are susceptible to deviation from the normality, it is recommended to test whether or not the coefficients sufficiently satisfy the normality. If the assumption of normality is violated significantly, it is better to choose any resampling procedure (e.g. bootstrap, permutation test, or cross-validation) or kernel partial least-square regression (KPLS; Iwata et al. 2015) for assessing the statistical significance. In some entomological studies, permutation procedures such as non-parametric MANOVA and cross-validation were used for the identification of groups (Márquez & Knowles 2007; Yang et al. 2015). Variables used in these resampling statistics do not need to conform to normality, independency or homogeneity of variance. For exploring relationships between explanatory variables and Fourier coefficients, KPLS is recommended to be carried out; however, as far as we know, there is no such example in the field of entomology. It is common to use PC analysis for summarizing Fourier coefficients (EF-PCA). By conducting EF-PCA, numerous coefficient data are reduced into a few principal components, and then the shape variation can be reconstructed by inverse Fourier transformation along these principal component axes, which enables us to interpret the biological meaning of each axis easily (Rohlf & Archie 1984). In PC analysis, no preliminary

176

check of the statistical property of EF coefficients is necessary because this method implements only the rotation of the raw data (Zelditch et al. 2004). After EF-PCA, the subset of principal component scores may be used for subsequent statistical analyses (e.g. ANOVA, MANCOVA, discriminant analysis and cluster analysis). However, before conducting these analyses in a conventional manner (i.e. under the assumption of normality), it is preferable to check whether the distribution of PC scores satisfies the assumption(s). If the resampling procedures mentioned above are chosen to analyze PC scores, there is no need to check the assumption(s).

USE OF GEOMETRIC MORPHOMETRICS APPROACHES IN ENTOMOLOGY: A 15-YEAR TREND Morphology of various body parts of insects has been shown to be associated with function. As insects, in general, have firm exoskeletons and their morphology does not change much after developing into the adult stage, shape quantification based on geometric morphometrics has been an effective approach to describe shape variation. To understand how geometric morphometrics has been utilized in entomological studies, a bibliographic survey was carried out using Google Scholar, a search engine for academic publications. The search was undertaken for articles that used geometric morphometrics to describe the shape of insects and published between 2000 and 2015. The key words used for the search were “geometric morphometrics” and “insect” or “entomology”. The papers were then screened and broken down by the order of the target insects, body parts analyzed and the geometric morphometric methods used. In this process, the studies in which geometric morphometric methods were not used or those written in languages other than English were excluded. Also excluded were a small number of papers that were difficult to obtain. When there were multiple target body parts and/or multiple approaches used in a single article, the individual sets of body part and geometric morphometric approaches were counted as separate cases. The bibliographic survey resulted in 472 articles that used geometric morphometric methods to describe the shape of insects. The number of papers has continued to increase linearly approximately 20-fold between 2000 and 2015 (Fig. 9). Those studies of geometric morphometrics were biased in number to a few orders of insects: Diptera being the top, followed by Coleoptera and Hymenoptera (Table 1). Within Diptera, studies on the two most important groups, Drosophila flies



Figure 9 Yearly publication of academic articles regarding applied geometric morphometric approaches to analyses of insect shapes between 2000 and 2015, based on a search using Google Scholar.

and mosquitos, accounted for more than half of the studies on this order. The body parts studied by morphometric analyses differed among the three wellstudied orders; wings were major target parts in Diptera and Hymenoptera, whereas wings, genitalia, head, thorax, whole body shape and elytra were evenly studied in Coleoptera (Table 1). The wing shape was also rather well subjected to geometric morphometric analyses in other insect orders such as Lepidoptera, Hemiptera, Odonata and Orthoptera (Table 1). Four major geometric morphometric approaches were used in the papers: landmark-based approach, semi-landmark-

based approach, landmark- and semi-landmark-based approach, and EF analysis (Table 1). The wings of insects, the body part most subjected to morphometric analyses, are greatly deformed in 3-D on flapping in flight to generate the necessary lift for flying. Unlike vertebrate wings, insect wings have flight muscles concentrated at their bases, so that the wings are passively deformed by flapping. Therefore, the wing structure and its material properties are major factors determining the flight characteristics of insects (Wooton 1992). In the wings of insects, chitin strips called wing veins are the main support structure. The arrangement and branching patterns of these wing veins are largely different among various insect taxonomic groups and are often used for classification, especially at the order and family levels. Some studies have shown that such differences in the wing vein pattern lead to functional differentiation. Furthermore, it has been confirmed that the wing vein pattern affects wing flexibility in a wide range of insect groups (Ennos 1988; Combes & Daniel 2003), and probably flight capability in consequence. A recent study by Ray et al. (2016) suggests that wing shape deformation significantly affects aerial agility in D. melanogaster. In many insect taxa, the patterns of wing vein location can be readily characterized by placing homologous landmarks on the intersection of the wing veins (Miguel et al. 2010; Takahashi et al. 2011b; Habel et al. 2012); an example (D. melanogaster wing) shown in Figure 5 (a). Therefore, the homologous landmark-based approach was applied in more than 95% of the geometric morphometric studies on insects (Table 1).

Table 1 Published works that applied different methods of geometric morphometrics to shape analyses in different orders of insects Body part:

Wing

Genitalia

Head

Thorax

Body shape

Elytra

Mandible

Others

Total

8 16 3 19 2 0 4 7

1 14 3 2 14 0 1 1

1 20 5 0 0 0 1 1

0 10 2 0 4 0 1 6

0 16 0 0 0 0 0 0

0 4 3 1 0 0 2 1

8 13 4 4 3 0 2 4

163 107 99 63 55 25 14 21

26 12 15 4 2 59

24 5 1 3 3 36

14 8 3 3 0 28

19 1 0 0 3 23

7 7 2 0 0 16

5 2 3 1 0 11

23 9 5 1 0 38

439 55 29 12 12 547

Insect order Diptera 145 Coleoptera 14 Hymenoptera 79 Lepidoptera 37 Hemiptera 32 Odonata 25 Orthoptera 3 Others 1 Geometric morphometric approach Landmark 321 Landmark & semi-landmark 11 Elliptical Fourier analysis 0 Semi-landmark 0 Others 4 Total 336

Each studied body part was counted separately. Results are based on a Google Scholar search for academic articles published between 2000 and 2015 using the key words “geometric morphometrics” and “insect” or “entomology”. The articles were excluded if geometric morphometric methods were not used, they were written in languages other than English or they were difficult to obtain.


177

H. Tatsuta et al.

When the landmark-based approach is applied to insect wings, a tedious and error-prone process is manually placing of several landmarks to cover the whole wing region. To date, several attempts have been made to develop semi-automated morphometric toolkits to deal with this problem. Dobens and Dobens (2013) developed a set of macros designed to perform a semiautomated morphometric analysis of Drosophila wing images. Although these macros are not designed to obtain the coordinates of homologous landmarks, they do allow the user to count trichomes and measure the wing area semi-automatically. Houle et al. (2003) developed a semi-automated image analysis system, called “Wingmachine”, to obtain the x and y coordinates of homologous landmarks on Drosophila wings. In this system, an a priori B-spline model was fitted to each of the wing images using the pixel brightness of the reversed and filtered images (Lu & Milios 1994; Houle et al. 2003). So far, any fully automated morphometric tool applicable to general insect wings has not been developed, and the effort to place landmarks on wings manually is not avoidable for many insect taxa. Genital morphology of insects, the second wellstudied body part in morphometrics, has been used in taxonomy related, for example, to the morphological and biological species concepts (Shapiro & Porter 1989). Male genitalia have been well studied as important organs providing diagnoses at the species level in many insects, whereas female genitalia, in general, are harder to study because they are often soft and involuted (Eberhard 1985). This was confirmed by the search of published works: among a total of 59 geometric morphometric analyses of insect genitalia published between 2000 and 2015 (Table 1), studies of male genitalia accounted for approximately 85%. Geometric morphometric analyses on insect genital structures were not only used for taxonomy and/or morphological description, but also for investigating reproductive isolation between species. The diverse genital morphology among insects has the potential to contribute to reproductive isolation between species through two different mechanisms: (i) mechanical lockand-key process; and (ii) sensory lock-and-key process (Eberhard 1992; Masly 2012). In the mechanical or structural lock-and-key process, mechanical incompatibility of genitalia, which is due to interspecific differences in genital morphology, prevents or reduces mating success between allospecific pairs. In the sensory lock-and-key process formalized by Eberhard (1992), interspecific differences in genital morphology are perceived by one or both sexes and result in premature termination of mating or post-copulatory

178

reproductive fitness problems. Although many attempts have been made to test this hypothesis, the majority of such studies have failed to find convincing evidence supporting it (Masly 2012), reflecting that the relationship between genital morphology and its function remains unclear. In contrast to wings, the approach based on homologous landmarks was less (only 44%) applied in the 59 studies on genitalia, but outline-based approaches using semi-landmarks or EF analysis instead accounted for approximately 32% of the studies (Table 1). This equivalent use of landmark-based and outline-based approaches seems to reflect variation in the genital morphology. As generally suggested for geometric morphometry, the approach based on homologous landmarks is recommended for quantifying the shape of genitalia as well, if sufficient landmarks are available. Such landmarks have often been placed on the tips of protrusions and/or at junctions of multiple segments. However, if such landmarks are insufficient, the best approach to be used depends on the structural features of the genitalia. When only few homologous landmarks are available for sufficiently covering the entire genital morphology, some semi-landmarks can be added; such combined analyses accounted for approximately 20% of the 59 studies on genitalia (Table 1). However, when distinct structures are not available for placing landmarks, outline-based approaches could be applied.

SOFTWARE AND USEFUL TOOLKITS For many entomologists, the intricate calculations using linear algebra that are necessary to carry out a geometric morphometric analysis is seen as forbidding. Fortunately, many useful computer routines and userfriendly software have now been developed by various researchers. One way of searching the available resources is to visit the Morphometrics site at SUNY Stony Brook (http://life.bio.sunysb.edu/morph/). The graphical user interface (GUI)-based software provided in the website, most of which have been developed by F. James Rohlf, is particularly useful for data acquisition and conducting landmark-based analyses. Recently, the R statistical software, with a variety of R source codes and packages, has also been developed and rapidly spread. Here is a brief introduction to several representative routines using R software and its R packages.

Landmark-based approaches With respect to landmark-based approaches, coordinates of landmarks need to be first acquired. This



entails the use of digitizer software once the images of interest have been obtained. Coordinates of specified points in an image can be obtained by any general image processing software such as ImageJ (http:// imagej.nih.gov/ij/) or any digitizer software designed specifically for acquiring landmark coordinates for geometric morphometric analysis. One of the most frequently used digitizer software programs for acquiring landmark coordinates is tpsDIG2, developed as one of the tps software routines by Rohlf (2015). An advantage of using tpsDIG2 over ImageJ is that it can export the coordinate data in a format that can be directly processed by tps software for further geometric morphometric analyses. Landmark coordinates can also be acquired by using the digitizer function in the geomorph package of a specific statistical software R routine (Adams & Otarola-Castillo 2013). As in the case of tpsDIG2, this also allows for a direct transfer of the coordinate data for further analyses in R. One criterion for determining which software or R package is to be used for landmark coordinate acquisition is that file format transformation is not necessary to transfer coordinate data from software to software. After acquiring landmark coordinates, the landmark configurations may need to be superimposed. Software such as tpsSuper (Rohlf 2016), MorphoJ (Klingenberg 2011) or R packages such as geomorph (Adams & Otarola-Castillo 2013) and shapes (Dryden & Mardia 2016) can perform GLS. Both GLS and resistant-fit superimposition such as RFTRA can be performed with CoordGen and SuperPoser in the IMP suite (Integrated Morphometrics Package, Sheets 2001). At the moment, no R package is available for resistant-fit superimposition, but Claude (2008) provided several helpful R functions for resistant-fit superimposition for general users. Once Procrustes coordinates are obtained, basic geometric morphometric analyses such as thin-plate spline analysis and relative warp analysis can be performed using software routines such as tpsSpline and tpsRelw32 (http://life.bio.sunysb.edu/morph/), MorphoJ, or the IMP suite; other useful R packages are geomorph or shapes. In addition, some software or R packages cover further advanced analyses: MorphoJ specializes in phylogenetics, quantitative genetics and analyses of modularity in shape data; Evomorph, an R package developed by Cabrera and Giri (https://cran.rproject.org/web/packages/Evomorph/index.html), specializes in the simulation of evolutionary process using geometric morphometric data; and TNT developed by Goloboff et al. (2008) calculates the most parsimonious cladogram from landmark coordinate data (Catalano et al. 2010).


Semi-landmark-based approaches To evaluate outline shape by semi-landmark-based methods, semi-landmark coordinates need to be acquired. There are several ways to delimit segments of the outline under analysis and different software routines provide different options. As in the case of landmark coordinate acquisition, tpsDIG2 and geomorph can be used for digitizing semi-landmark coordinates and they delimit the outline by increments along the length of the curve of interest. Two programs, MakeFan and SemiLand, in the IMP suite are also designed for the purpose of acquiring semi-landmark coordinates. They offer two options to delimit the outline: MakeFan draws fans at equal angular intervals on an image, whereas SemiLand slides the semi-landmarks along the curve to minimize the Procrustes distance between the subject (target) and a reference (base). Superimposition of the semi-landmark configurations can be carried out in these programs or geomorph. Once the coordinates of the superimposed configurations are obtained, multivariate analysis such as PCA or discriminant functions can be applied to evaluate the shape variation and the separation of groups.

Elliptical Fourier analysis and related analyses for outlines For the EF method, integrated packages with a digitizer and software or functions for EF analysis are available. In Momocs (modern morphometrics) (https://cran.rproject.org/web/packages/Momocs/index.html), there is an R package for EF analysis, which includes functions for performing chain coding of contour information and calculating EF descriptors, and series of analyses including multivariate ones such as PCA (Bonhomme et al. 2014). Another integrated software package for EF analysis is SHAPE (http://lbm.ab.a.u-tokyo.ac.jp/ ~iwata/shape/, Iwata & Ukai 2002), where all the software programs required for EF analysis are included as in Momocs. Both Momocs and SHAPE are valuable tools to perform EF analyses; SHAPE is more interactive and useful, especially when outline acquisition needs to be checked individually for each sample. However, only bitmap files can be handled in SHAPE and only a Windows OS version is available. For UNIX or MacOS users, it would be feasible to access the page of Shape on R developed by Hiroyoshi Iwata (http://lbm. ab.a.u-tokyo.ac.jp/~iwata/software/shape_r/), where the user can take the R source code together with instructions (unfortunately, only in Japanese!) and the Perl script that can convert the data format of SHAPE into an appropriate format. Finally, another package to be

179

H. Tatsuta et al.

mentioned is MLmetrics, which is available at a nominal cost from [email protected]. The on-screen digitizing routine has now been ported over to MS Windows obviating the need for a separate digitizer. The set of DOS routines is currently in the process of being converted into an MS Windows version. This package allows the user to compute EF function harmonics, predicted points, residuals and ellipse orientations, providing a plethora of other useful functions as well as facilitating the transfer of computed data to other statistical or graphical-plotting packages.

CONCLUDING REMARKS: COMMENTS FOR 3-D AND FINE-SCALE FEATURES Geometric morphometrics technology has gained considerable momentum and is making rapid progress. In entomology, 2-D data have been the primary sources for analyses, particularly on flat organs such as the wings and the posterior lobe of epandrium (a part of male genitalia) in Drosophila (Zeng et al. 2000; Takahashi et al. 2011a), and even on cubic organs such as mandibles of stag beetles (Tatsuta et al. 2004). Obviously, in the latter case, the comparison of shape features should ideally be analyzed in 3-D. The use of 2-D projections of 3-D-structured body parts may lead to possible distortion as well as a substantial loss of morphological information. Recent advances in high-precision X-ray Micro Computed Tomography (Micro-CT) have enabled us to compare 3-D external and internal structures of body parts and organs, such as brain or muscles, in a reasonable fashion (Wipfler et al. 2016). This has now been extensively utilized in the comparison of skulls, bones and teeth in paleobiology, anthropology and ornithology (e.g. Araujo et al. 2017; Bribiesca-Contreras & Sellers 2017; Dou et al. 2017). This approach also displays its usefulness and further progress in the field of entomology (Friedrich et al. 2014). The use of nondestructive inspection of materials was implemented, for instance, for the identification of different social castes in nests (Aguilera-Olivares et al. 2017) and for tracking the movement of insects inside nests (Johnson et al. 2004; Jennings & Austin 2011). As verified in various organisms, it is now feasible to inspect the ultrastructure of internal tissues or muscles (David et al. 2016; Fabian et al. 2016; Smith et al. 2016; Swart et al. 2016) and to visualize the process of metamorphosis (Lowe et al. 2013; Martin-Vega et al. 2016; Hall et al. 2017). In studies of sexual selection, it has almost been impossible to directly observe reproductive organs and tract in vivo during the insemination process, but now time-series changes of concurrent

180

organ rearrangements during the process within the abdomen can be observed using Micro-CT (Mattei et al. 2015). Furthermore, this tool also enables direct observation of interactions between the male intromittent and the female copulatory organs or the process of mating-plug formation (Uhl et al. 2014; Dougherty et al. 2015). As revealed by these lines of evidence, Micro-CT has great potential for understanding of proximate factors that can influence individual fitness during mating. Although 3-D quantitative comparisons have not extensively been applied in entomology yet, such studies on reproductive organs look promising in the future. The long-term controversial debate of cryptic female choice and sexual conflict is beginning to unravel (Orr & Brennan 2015; Firman et al. 2017), with the extensive focus on potential target organs of sexual selection in conjunction with the aid of useful tools in geometric morphometrics.

ACKNOWLEDGMENTS We cordially thank Pete Lestrel for critical reading of the earlier versions of the manuscript and providing thoughtful comments. Thanks are also due to Hiroyoshi Iwata, Chris Klingenberg, and Norman MacLeod for enthusiastic discussion and comments on the contents of the paper. This study was supported in part by Grants-in-Aid for Scientific Research (Nos. JP17H05014 to K.H.T., 24657061 to Y.S., and JP17H03722 and 25304014 to H.T.) from the Japan Society for the Promotion of Science.

REFERENCES Adams DC, Otarola-Castillo E (2013) geomorph: an R package for the collection and analysis of geometric morphometric shape data. Methods in Ecology and Evolution 4, 393–399. Aguilera-Olivares D, Palma-Onetto V, Flores-Prado L, Zapata V, Niemeyer HM (2017) X-ray computed tomography reveals that intraspecific competition promotes soldier differentiation in a one-piece nesting termite. Entomologia Experimentalis et Applicata 163, 26–34. Araujo R, Fernandez V, Polcyn MJ et al. (2017) Aspects of gorgonopsian paleobiology and evolution: insights from the basicranium, occiput, osseous labyrinth, vasculature, and neuroanatomy. PeerJ 5, Article ID e3119. https://doi. org/10.7717/peerj.3119 Atchley WR, Anderson D (1978) Ratios and the statistical analysis of biological data. Systematic Zoology 27, 71–78. Atchley WR, Gaskins CT, Anderson D (1976) Statistical properties of ratios. I. Empirical results. Systematic Zoology 25, 137–148.



Beerends RJ, Morsche HG, Van den Berg JC, Van de Vrie EM (2003) Fourier and Laplace Transforms. Cambridge University Press, Cambridge. Berner D (2011) Size correction in biology: how reliable are approaches based on (common) principal component analysis? Oecologia 166, 961–971. Bolnick DI, Lau OL (2008) Predictable patterns of disruptive selection in stickleback in postglacial lakes. The American Naturalist 172, 1–11. Bonhomme V, Picq S, Gaucherel C, Claude J (2014) Momocs: outline analysis using R. Journal of Statistical Software 56, 1–24. Bookstein FL (1991) Morphometric Tools for Landmark Data: Geometry and Biology. Cambridge University Press, Cambridge. Bookstein FL (1996) Combining the tools of geometric morphometrics. In: Marcus LF, Corti M, Loy A, Naylor GJP, Slice D (eds) Advances in Morphometrics, pp 131–151. Plenum Press, New York, NY. Bribiesca-Contreras F, Sellers WI (2017) Three-dimensional visualisation of the internal anatomy of the sparrowhawk (Accipiter nisus) forelimb using contrast-enhanced microcomputed tomography. PeerJ 5, Article ID e3039. https:// doi.org/10.7717/peerj.3039 Brower DL (1986) Engrailed gene expression in Drosophila imaginal discs. EMBO Journal 5, 2649–2656. Burnaby TP (1966) Growth-invariant discriminant functions and generalized distances. Biometrics 22, 96–110. Catalano SA, Goloboff PA, Giannini NP (2010) Phylogenetic morphometrics (I): the use of landmark data in a phylogenetic framework. Cladistics 26, 539–549. Claude J (2008) Morphometrics with R. Springer Science +Business, New York, NY. Cohen B, Wimmer EA, Choen SM (1991) Early development of lge and wing primordia in the Drosophila embryo. Mechanisms of Development 33, 229–240. Combes SA, Daniel TL (2003) Flexural stiffness in insect wings. I. Scaling and the influence of wing venation. Journal of Experimental Biology 206, 2979–2987. David S, Funken J, Potthast W, Blanke A (2016) Musculoskeletal modelling of the dragonfly mandible system as an aid to understanding the role of single muscles in an evolutionary context. Journal of Experimental Biology 219, 1041–1049. de Boor C (1978) A Practical Guide to Splines. Springer, New York, NY. Dobens AC, Dobens L (2013) FijiWings: an open source toolkit for semiautomated morphometric analysis of insect wings. G3:Genes,Genomes, Genetics 3, 1443–1449. Dou L, Li DJ, Xu TT, Tang Y, Yang D (2017) Root anatomy and canal morphology of mandibular first premolars in a Chinese population. Scientific Reports 7, Article ID 750. https://doi.org/10.1038/s41598-017-00871-9 Dougherty LR, Rahman IA, Burdfield-Steel ER et al. (2015) Experimental reduction of intromittent organ length reduces male reproductive success in a bug. Proceedings


of the Royal Society B – Biological Sciences 282, Article ID 20150724. https://doi.org/10.1098/rspb.2015.0724 Dryden I, Mardia K (1998) Statistical Shape Analysis. John Wiley & Sons, Inc., New York, NY. Dryden I, Mardia K (2016) Statistical Shape Analysis: With Applications in R, 2nd edn. John Wiley & Sons, Inc., New York, NY. Dworkin I, Gibson G (2006) Epidermal growth factor receptor and transforming growth factor-beta signaling contributes to variation for wing shape in Drosophila melanogaster. Genetics 173, 1417–1431. Eberhard WG (1985) Sexual Selection and Animal Genitalia. Harvard University Press, Cambridge, MA. Eberhard WG (1992) Species isolation, genital mechanics, and the evolution of species-specific genitalia in three species of Macrodactylus beetles. Evolution 46, 1774–1783. Eberle J, Walbaum W, Warnock RCM, Fabrizi S, Ahrens D (2015) Asymmetry in genitalia does not increase the rate of their evolution. Molecular Phylogenetics and Evolution 93, 180–187. Engel H (1986) A least-squares method for estimation of Bezier curves and surfaces and its applicability to multivariate analysis. Mathematical Biosciences 79, 155–170. Ennos R (1988) The importance of torsion in the design of insect wings. Journal of Experimental Biology 140, 137–160. Evans DG, Schweitzer PN, Hanna MS (1985) Parametric cubic splines and geological shape descriptions. Mathematical Geology 17, 611–624. Fabian B, Schneeberg K, Beutel RG (2016) Comparative thoracic anatomy of the wild type and wingless (wg1cn1) mutant of Drosophila melanogaster (Diptera). Arthropod Structure & Development 45, 611–636. Ferson S, Rohlf FJ, Koehn RK (1985) Measuring shape variation of two-dimensional outlines. Systematic Zoology 34, 59–68. Firman RC, Gasparini C, Manier MK, Pizzari T (2017) Postmating female control: 20 years of cryptic female choice. Trends in Ecology and Evolution 32, 368–382. Friedrich F, Matsumura Y, Pohl H, Bai M, Hörnschemeyer T, Beutel RG (2014) Insect morphology in the age of phylogenomics: innovative techniques and its future role in systematics. Entomological Science 17, 1–24. Fukudome H, Sakamaki Y (2011) Exploring new diagnostic characters using elliptic Fourier descriptors. Japanese Journal of Entomology (New Series) 14, 232–241. (In Japanese with English summary.) Garcia-Bellido A, Santamaria P (1972) Developmental analysis of the wing disc in the mutant engrailed of Drosophila melanogaster. Genetics 72, 87–104. Goloboff PA, Farris JS, Nixon KC (2008) TNT, a free program for phylogenetic analysis. Cladistics 24, 774–786. Goodall CR (1991) Procrustes methods in the statistical analysis of shape. Journal of the Royal Statistical Society, Series B 53, 285–339. Granier S, Magniez-Jannin F, Rasplus JY, Alibert P (2005) When morphometry meets genetics: inferring the phylogeography of Carabus solieri using Fourier analyses of

181

H. Tatsuta et al.

pronotum and male genitalia. Journal of Evolutionary Biology 18, 269–280. Gurdon JB, Bourillot PY (2001) Morphogen gradient interpretation. Nature 413, 797–803. Habel JC, Engler JO, Rodder D, Schmitt T (2012) Contrasting genetic and morphologic responses on recent population decline in two burnet moths (Lepidoptera, Zygaenidae). Conservation Genetics 13, 1293–1304. Hall MJR, Simonsen TJ, Martin-Vega D (2017) The ‘dance’ of life: visualizing metamorphosis during pupation in the blow fly Calliphora vicina by X-ray video imaging and micro-computed tomography. Royal Society Open Science 4, Article ID 160699. https://doi.org/10.1098/rsos. 16069 Held LI Jr (2002) Imaginal Discs. Cambridge University Press, Cambridge. Hidalgo A (1994) Three distinct roles for the engrailed gene in Drosophila wing development. Current Biology 4, 1087–1098. Houle D, Mezey J, Galpern P, Carter A (2003) Automated measurement of Drosophila wings. BMC Evolutionary Biology 3, Article ID 25. https://doi.org/10.1186/14712148-3-25 Humphries JM, Bookstein FL, Chernoff B, Smith GR, Elder RL, Poss SG (1981) Multivariate discrimination by shape in relation to size. Systematic Zoology 30, 14–29. Iwata H, Ukai Y (2002) SHAPE: a computer program package for quantitative evaluation of biological shapes based on elliptic Fourier descriptors. Journal of Heredity 93, 384–385. Iwata H, Ebana K, Uga Y, Hayashi T (2015) Genomic prediction of biological shape: elliptic Fourier analysis and kernel partial least squares (PLS) regression applied to grain shape prediction in rice (Oryza sativa L.). PLoS ONE 10, Article ID e0120610. https://doi.org/10.1371/journal.pone. 0120610 Jacobshagen B (1997) Craniofacial variability in the Hominoidea. In: Lestrel PE (ed.) Fourier Descriptors and their Application in Biology, pp 227–249. Cambridge University Press, Cambridge. Jennings JT, Austin AD (2011) Novel use of a microcomputed tomography scanner to trace larvae of wood boring insects. Australian Journal of Entomology 50, 160–163. Johnson SN, Read DB, Gregory PJ (2004) Tracking larval insect movement within soil using high resolution X-ray microtomography. Ecological Entomology 29, 117–122. Kendall D (1977) The diffusion of shape. Advances in Applied Probability 9, 89–96. Kendall DG (1984) Shape manifolds, Procruestean metrics and complex projective spaces. Bulletin of the London Mathematical Society 16, 81–121. Kendall DG, Barden D, Carne TK, Le H (1999) Shape and Shape Theory. John Wiley & Sons, Inc., New York, NY. Klingenberg CP (2008) Morphological integration and developmental modularity. Annual Review of Ecology, Evolution and Systematics 39, 115–132.

182

Klingenberg CP (2010) Evolution and development of shape: integrating quantitative approaches. Nature Reviews Genetics 11, 623–635. Klingenberg CP (2011) MorphoJ: an integrated software package for geometric morphometrics. Molecular Ecology Resources 11, 353–357. Klingenberg CP (2016) Size, shape and form: concepts of allometry in geometric morphometrics. Development Genes and Evolution 226, 113–137. Kuhl FP, Giardina CR (1982) Elliptic Fourier features of a closed contour. Computer Graphics and Image Processing 18, 236–258. Lande R, Arnold SJ (1983) The measurement of selection on correlated characters. Evolution 37, 1210–1226. Langlade NB, Feng X, Dransfield T et al. (2005) Evolution through genetically controlled allometry space. Proceedings of the National Academy of Sciences USA 102, 10221–10226. Laparie M, Vernon P, Cozic Y, Frenot Y, Renault D, Debat V (2016) Wing morphology of the active flyer Calliphora vicina (Diptera: Calliphoridae) during its invasion of a sub-Antarctic archipelago where insect flightlessness is the rule. Biological Journal of the Linnean Society 119, 179–183. Lele S (1993) Euclidian distance matrix analysis (EDMA): estimation of mean form and mean form difference. Mathematical Geology 25, 573–602. Lestrel PE (1997) Fourier Descriptors and their Applications in Biology. Cambridge University Press, Cambridge. Lestrel PE (2000) Morphometrics for the Life Sciences. World Scientific, Singapore. Lestrel PE, Wolfe CA (2017) Shape differences in the PlioPleistocene hominin cranial vault: elliptical Fourier descriptors. In: Lestrel P (ed.) Biological Shape Analysis. Abstracts of the 5th International Symposium on Biological Shape Analysis, pp 23–26. University of Tokyo, Tokyo. Libungan LA, Pálsson S (2015) ShapeR: an R package to study otolith shape variation among fish populations. PLoS ONE 10, Article ID e0121102. https://doi.org/10. 1371/journal.pone.0121102 Lohman GP (1983) Eigenshape analysis of microfossils: a general morphometric procedure for describing changes in shape. Mathematical Geology 15, 659–672. Lohman GP, Schweitzer PN (1990) On eigenshape analysis. In: Rohlf FJ, Bookstein FL (eds) Proceedings of the Michigan Morphometrics Workshop, Special Publications No. 2, pp 147–166. University of Michigan Museum of Zoology, Ann Arbor, MI. Lowe T, Garwood RJ, Simonsen TJ, Bradley RS, Withers PJ (2013) Metamorphosis revealed: time-lapse threedimensional imaging inside a living chrysalis. Journal of the Royal Society Interface 10, Article ID 20130304. https://doi.org/10.1098/rsif.2013.0304 Lowrence PA, Morata G (1976) Compartments in the wing of Drosophila: a study of the engrailed gene. Developmental Biology 50, 321–337.



Loy A (2007) Morphometrics and theriology. Homage to Marco Corti. Hystrix the Italian Journal of Mammalogy 18, 115–136. Lu F, Milios EE (1994) Optimal spline fitting to planar shape. Signal Processing 37, 129–140. Lu KH (1965) Harmonic analysis of the human face. Biometrics 21, 491–505. MacLeod N (1999) Generalizing and extending the eigenshape method of shape space visualization and analysis. Paleobiology 25, 107–138. Marcus LF, Bello E, García-Valdecasas A, (eds) (1993) Contributions to Morphomerics, Monografias del Museo Nacional de Ciencias Naturales 8. Consejo Superior de Investigaciones Científicas, Madrid. Marcus LF, Corti M, Loy A, GJP N, Slice DE (1996) Advances in Morphometrics. Plenum Press, New York, NY. Márquez EJ, Knowles LL (2007) Correlated evolution of multivariate traits: detecting co-divergence across multiple dimensions. Journal of Evolutionary Biology 20, 2334–2348. Martin-Vega D, Simonsen TJ, Hall MJR (2016) Looking into the puparium: micro-CT visualization of the internal morphological changes during metamorphosis of the blow fly, Calliphora vicina, with the first quantitative analysis of organ development in cyclorrhaphous dipterans. Journal of Morphology 278, 629–651. Masly JP (2012) 170 years of “lock-and-key”: genital morphology and reproductive isolation. International Journal of Evolutionary Biology 2012, Article ID 247352. https://doi.org/10.1155/2012/247352 Mattei AL, Riccio ML, Avila FW, Wolfner MF (2015) Integrated 3D view of postmating responses by the Drosophila melanogaster female reproductive tract, obtained by micro-computed tomography scanning. Proceedings of the National Academy of Sciences USA 112, 8475–8480. Medal R, Botto-Mahan C, Kalin-Arroyo M (2003) Pollinatormediated selection on the nectar guide phenotype in the Andean monkey flower, Mimulus luteus. Ecology 84, 1721–1732. Miguel I, Baylac M, Iriondo M, Manzano C, Garnery L, Estonba A (2010) Both geometric morphometric and microsatellite data consistently support the differentiation of the Apis mellifera M evolutionary branch. Apidologie 42, 150–161. Minaka N (1999) [Morphometrics]. In: Tanabe K, Mori K (eds) [Morphology and Analysis of Paleoorganism], pp 60–99. Asakura Publisher, Tokyo. (In Japanese.) Mood AM, Graybill FA, Boes DC (1974) Introduction to the Theory of Statistics, 3rd edn. McGraw-Hill, Singapore. Mosimann JE, James FC (1979) New statistical methods for allometry with application to Florida red-winged blackbirds. Evolution 33, 444–459. Orr TJ, Brennan PLR (2015) Sperm storage: distinguishing selective processes and evaluating criteria. Trends in Ecology and Evolution 30, 261–272. Palsson A, Gibson G (2004) Association between nucleotide variation in Egfr and wing shape in Drosophila melanogaster. Genetics 167, 1187–1198.


Penrose LS (1954) Distance, size and shape. Annals of Eugenics 18, 337–343. Persoon E, Fu K-S (1986) Shape discrimination using Fourier descriptors. IEEE Transactions on Pattern Analysis and Machine Intelligence 8, 388–397. Ray RP, Nakata T, Henningsson P et al. (2016) Enhanced flight performance by genetic manipulation of wing shape in Drosophila. Nature Communications 7, Article ID 10851. https://doi.org/10.1038/ncomms10851 Rohlf FJ (1990a) Rotational fit (Procrustes) methods. In: Rohlf FJ, Bookstein FL (eds) Proceedings of the Michigan Morphometrics Workshop, Special Publication No. 2, pp 227–236. The University of Michigan Museum of Zoology, Ann Arbor, MI. Rohlf FJ (1990b) Fitting curves to outlines. In: Rohlf FJ, Bookstein FL (eds) Proceedings of the Michigan Morphometrics Workshop, Special Publication No. 2, pp 167–177. The University of Michigan Museum of Zoology, Ann Arbor, MI. Rohlf FJ (1996) Morphometric spaces, shape components and the effects of linear transformations. In: Marcus LF, Corti M, Loy A, GJP N, Slice DE (eds) Advances in Morphometrics, pp 117–129. Plenum Press, New York, NY. Rohlf FJ (2005) Geometric morphometrics simplified. Trends in Ecology and Evolution 20, 13–14. Rohlf FJ (2015) The tps series of software. Hystrix the Italian Journal of Mammalogy 26, 9–12. Rohlf FJ (2016) tpsSuper. Morphometrics at SUNY Stony Brook. Available from URL: http://life.bio.sunysb.edu/ morph/index.html Rohlf FJ, Archie JW (1984) A comparison of Fourier methods for the description of wing shape in mosquitoes (Diptera: Culicidae). Systematic Zoology 33, 302–317. Rohlf FJ, Bookstein FL (1987) A comment on shearing as a method for “size correction”. Systematic Zoology 36, 356–367. Rohlf FJ, Marcus LF (1993) A revolution in morphometrics. Trends in Ecology and Evolution 8, 129–132. Sanicola M, Sekelsky J, Elson S, Gelbart WM (1995) Drawing a stripe in Drosophila imaginal disks: negative regulation of Decapentaplegic and patched expression by engrailed. Genetics 139, 745–756. Shapiro AM, Porter AH (1989) The lock-and-key hypothesis: evolutionary and biosystematic interpretation of insect genitalia. Annual Review of Entomology 34, 231–245. Shapiro SS, Wilk MB (1965) An analysis of variance test for normality (complete samples). Biometrika 52, 591–611. Sheets DH (2001) IMP Suite (Integrated Morphometrics Package). Available from URL: http://www3.canisius. edu/~sheets/IMP%208.htm Siegel AF, Benson RH (1982) A robust comparison of biological shapes. Biometrics 38, 341–350. Smith DB, Bernhardt G, Raine NE et al. (2016) Exploring miniature insect brains using micro-CT scanning techniques. Scientific Reports 6, Article ID 21768. https://doi. org/10.1038/srep21768

183

H. Tatsuta et al.

Suzuki TK (2013) Modularity of a leaf moth-wing pattern and a versatile characteristic of the wing-pattern ground plan. BMC Evolutionary Biology 13, Article ID 158. https://doi.org/10.1186/1471-2148-13-158 Swart P, Wicklein M, Sykes D, Ahmed F, Krapp HG (2016) A quantitative comparison of micro-CT preparations in dipteran flies. Scientific Reports 6, Article ID 39380. https://doi.org/10.1038/srep39380 Tabata T, Kornberg TB (1994) Hedgehog is a signaling protein with a key role in patterning Drosophila imaginal discs. Cell 76, 89–102. Takahashi KH, Daborn PJ, Hoffmann AA, Takano-Shimizu T (2011a) Environmental stress-dependent effects of deletions encompassing Hsp70Ba on canalization and quantitative trait asymmetry in Drosophila melanogaster. PLoS ONE 6, Article ID e17295. https://doi.org/10.1371/ journal.pone.0017295 Takahashi KH, Okada Y, Teramura K, Tsujino M (2011b) Deficiency mapping of the genomic regions associated with effects on developmental stability in Drosophila melanogaster. Evolution 65, 3565–3577. Tatsuta H, Mizota K, Akimoto S (2004) Relationship between size and shape in the sexually dimorphic stag beetle Prosopocoilus inclinatus (Coleoptera: Lucanidae). Biological Journal of the Linnean Society 81, 219–233. Tatsuta H, Iwata H, Goka K (2013) Morphometric studies on the variation of male lucanid beetle mandibles. In: Lestrel P (ed.) Proceedings of the 2nd International

184

Symposium of the Biological Shape Analysis, pp 55–72. World Scientific, Singapore. Thompson DW (1917) On Growth and Form. Cambridge University Press, London. Uhl G, Kunz K, Vocking O, Lipke E (2014) A spider mating plug: origin and constraints of production. Biological Journal of the Linnean Society 113, 345–354. Wipfler B, Pohl H, Yavorskaya MI, Beutel RG (2016) A review of methods for analysing insect structures – the role of morphology in the age of phylogenomics. Current Opinion in Insect Science 18, 60–68. Wooton RJ (1992) Functional morphology of insect wings. Annual Review of Entomology 37, 113–140. Yang HP, Ma CS, Wen H, Zhan QB, Wang XL (2015) A tool for developing an automatic insect identification system based on wing outlines. Scientific Reports 5, Article ID 12786. https://doi.org/10.1038/srep12786 Zahn CT, Roskies RZ (1972) Fourier descriptors for plane closed curves. IEEE Transactions on Computers 21, 269–281. Zecca M, Basler K, Struhl G (1995) Sequential organizing activities of engrailed, hedgehog and decapentaplegic in the Drosophila wing. Development 121, 2265–2278. Zelditch ML, Swiderski DL, Sheets HD et al. (2004) Geometric Morphometrics for Biologists: A Primer. Elsevier Academic Press, New York, NY. Zeng Z-B, Liu J, Stam LF, Kao CH, Mercer JM, Laurie CC (2000) Genetic architecture of a morphological shape difference between two Drosophila species. Genetics 154, 299–310.
