Application of Shape Analyses to Recording Structures of Marine Organisms for Stock Discrimination and Taxonomic Purposes

Morphometric analysis of biogenic recording structures within marine organisms has applications in stock assessment, taxonomics, and ecomorphological studies, with shape variation markedly influenced by both genetics and the surrounding environment. Geometric morphometrics (GM) is an alternative approach to the “traditional” method of collecting linear measurements and applying multivariate statistical methods to these data. Landmark- and outline-based GM methods are suggested to have several advantages over the “traditional” method. Due to the increasing popularity of GM methods in the modern literature, this chapter first compares different morphometric techniques, and then reviews the methods applied to recording structures, with a focus on GM outline-based analyses. It is clear that outline methods have become a popular method of analysis for structures such as otoliths, particularly for the purpose of distinguishing between population components. However, for other structures such as beaks this technique is only in its early stages of application and is more difficult to apply but shows great promise for future studies. The advantages of using a holistic approach, incorporating several techniques including outline analysis for stock identification purposes is discussed.


Utilization of Recording Structures
Unlike terrestrial species, it is difficult to study marine animals in their natural environment (Arkhipkin, 2005), particularly when there are no obvious physical boundaries to delineate stocks or to limit species distribution. Recording structures hold a wealth of information that could elucidate many aspects of a species life history. These are structures within the body that contain periodical growth increments or layers that can be used to estimate the age and growth of an individual (Dodge and Thomson, 1974;Campana, 1999Campana, , 2001Arkhipkin et al., 2018;Hollyman et al., 2018). As well as age and growth, the shape of the structure can be used to distinguish between species and population components. In addition, the chemical composition of a recording structure can provide a chronological record of the environmental conditions experienced by an individual throughout its lifetime and, thus, a record of the chemical composition of the surrounding seawater at the time of elemental incorporation (Smith et al., 1979;Thresher, 1999;Arkhipkin, 2005;Gillanders, 2005;Semmens et al., 2007). Stable isotope analysis can also be applied to recording structures and is often used to determine trophic position within food webs, to investigate diet, and is occasionally used to infer migration (Kalish, 1991;Gillanders, 2005;Dance et al., 2014;Kerr and Campana, 2014;Kato et al., 2016;Rosas-Luis et al., 2017).
Examples of recording structures available in major taxonomic marine groups (mollusks, teleost fish, and elasmobranchs) can be found in Figure 1. Some of hard body parts are used as aging tools if their increments or layers are distinct, are deposited over a regular and determinable time-scale and the structure continues to grow throughout life (i.e., not subject to resorption or reworking) (Arkhipkin, 2005;Kerr and Campana, 2014). Not all hard body parts fulfill these criteria, e.g., if increments may be visible only for part of an individual's life history. For instance, in squid gladii, the first third of the rachis is unreadable (Arkhipkin et al., 2018). Alternatively, the outer surface could be eroded resulting in underestimation of age either as a result of feeding (in structures such as beaks, teeth, and baleen plates) (Klevezal, 1996;Doubleday et al., 2011;Rodríguez-Domínguez et al., 2013) or wear of the outer shell in bivalves and gastropods from contact with fishing gear or unsuccessful predator attacks (Hollyman et al., 2018). In addition, some structures are subject to resorption or in the case of scale, loss, regeneration (Kerr and Campana, 2014). In older teleost fish, especially epipelagic species (possibly due to their elevated metabolic activity), early annuli of fin rays and spines can be lost to vascularization and resorption of dermal bone (Kerr and Campana, 2014;Tzadik et al., 2017) and may be subject to resorption during periods of nutritional stress (Gillanders, 2001). However, evidence suggests that vertebrae grow by accretion with little remodeling and resorption (Kerr and Campana, 2014). Other recording structures amenable to aging include otoliths, statoliths, the shell (internal growth rings may be available for aging if there is damage to the outer surface), operculum, and eye lenses (Kerr and Campana, 2014;Hollyman et al., 2018). This chapter will focus on recording structures used in both aging and morphological studies, as material that is subject to aging is likely to be more available for study.
Variation in the shape of recording structures is markedly influenced by environment and genetics (Campana and Casselman, 1993;Arkhipkin, 2003;Cardinale et al., 2004;Hüssy, 2008;Vignon and Morat, 2010;Vignon, 2012). Because of this, exploring the morphometry (defined in this chapter as the analysis of the external shape of an object, shape variation between objects and the covariation of shape with other variables such as size) of recording structures has applications in both stock discrimination (Campana and Casselman, 1993;Tracey et al., 2006;Libungan and Pálsson, 2015) and systematics (Tuset et al., 2006;Fang et al., 2018b). Quantifying the morphometry of recording structures has advantages over soft body parts, which can be confounded by changes in animal condition and are subject to warping (Campana and Casselman, 1993;van der Vyver et al., 2016). The shape of recording structures was first analyzed by collecting measurements of linear distances and applying multivariate statistical methods to these data (traditional morphometrics) (Rohlf and Marcus, 1993). Some authors also calculated complicated indices based on distances and angles for stock discrimination purposes, but this is no longer a recommended method (Volpedo and Vaz-dos-Santos, 2015). This approach has its limitations, such as the loss of information by simplifying the shape (in case of all but some otoliths and statoliths the actual shape of the structure is much more complex than it's representation by the idealized form of the shape indices) and the risk of selecting dimensions that do not adequately represent the actual shape variation (Adams et al., 2004). Another issue with treatment of indices is a high degree of multicollinearity between the measurements. This means that, since many of measurements correlate with each other, the impact of the individual measurement decreases, as a result even a small modification of the single measurement could lead to significant change in the resulting indices. Different methods can help to reduce the impact of multicollinearity, specifically only measurements with low correlation could be selected or the data could be transformed using principal component analysis (Jolicoeur and Mosimann, 1960).
Any morphometric data is highly subjected to allometric alteration effects. So, if the ontogenetic modification of the structure's shape is not the aim of the study, the impact of the allometry should be avoided. There are several approaches to this issue: first of all, in case the sample size is large enough and it's possible to estimate the age of animals, the subsample of individuals of a same age cohort should be selected; on another hand the data of measurements could be transformed (Mosimann and James, 1979;Lleonart et al., 2000).
Geometric morphometric (GM) methods are an alternative way of analyzing and quantifying shape, which in theory retains more detail about the geometry of the structure than could be obtained from linear measurements (Adams et al., 2004). Rather than just reporting that shape has changed, in GMbased methods one can report that certain structures have moved relative to others (Rohlf and Marcus, 1993). Landmark-based GM methods (Cadrin, 2014) collect two-or three-dimensional Cartesian coordinates of discrete, homologous, and biologically definable landmarks (Zelditch et al., 2004). Outline-based GM methods (Cadrin, 2000) are principally used to quantify the boundary shape of an object so that patterns of shape variation within and among groups can be evaluated (Cadrin and Friedland, 1999;Libungan and Pálsson, 2015 and references within the introduction).
The purpose of this chapter is to briefly summarize the main GM techniques applied to recording structures, to discuss the main applications of these methods for recording structures that fulfill the main criteria of aging tools and to briefly compare the utility of each outline-based technique for each recording structure. Considerations for interpretation of these data will be highlighted.

Landmark-Based Geometric Morphometrics
Landmark-based methods collect two-or three-dimensional Cartesian coordinates of discrete, homologous anatomical loci (Zelditch et al., 2004). There are three principal types of landmarks as defined in Bookstein (1997): • Type 1: Discrete juxtaposition of tissues (e.g., points in space at which structures meet) • Type 2: Maxima of curvature or other local morphogenetic processes (location defined by obvious geometry, e.g., tips of claws or spines) • Type 3: Extremal points or constructed landmarks (points defined by other points, e.g., end points of diameters, centroids, landmarks on the medial axis of the outline).
Direct analysis of Cartesian coordinates as variables would be inappropriate, as the effects of variation in position, orientation, and scale of imaged specimens are still present (Adams et al., 2004). Non-shape variation can be mathematically removed prior to analysis using several different superimposition methods. For a comparison and discussion of the relative advantages and disadvantages of different superimposition methods, refer to chapter 5 of Zelditch et al. (2004). The use of Generalized Procrustes Analysis (GPA) as a superimposition method is widespread, with non-shape variation removed in the following manner: First, landmark configurations are translated so their centroid (the mean x-and y-coordinates of the landmarks for each specimen) is placed at the origin of the coordinate system. Configurations are then scaled (to unit centroid size) and rotated. Rotation is an iterative process with all landmark configurations rotated around their centroid repeatedly to minimize the squared distance between corresponding landmarks. All configurations are rotated to optimal alignment on the first configuration matrix (the reference configuration, this can be any specimen). The mean shape is then calculated for all rotated configurations and all configurations are rotated to optimal alignment of this mean shape (the new reference). At this stage, the mean shape is re-calculated on the new reference. When two consecutive mean shapes are the same, the process discontinues (Zelditch et al., 2004).
A limitation of landmark analysis is that a sufficient number of landmark locations may not be available on the structure in question to capture its shape, or that important shape differences may be located in the regions between landmarks. For example, there may be no obvious definable landmarks on a particularly rounded structure (Adams et al., 2004). In some cases, it may be appropriate to introduce sliding semi-landmarks to provide a richer description of the shape (Adams et al., 2004).

Outline-Based Geometric Morphometrics
Sometimes it may be more appropriate to use an outline-based approach than a landmark-based approach depending on the structure in question. For a simple structure with no obvious landmarks or conversely a complex, irregular structure where landmarks (and semi-landmarks) may not capture the true complexity of the object in question (Haines and Crampton, 2000), outline-based methods may be more suitable. Commonly applied outline methods are Fourier transform and Wavelet transform (WT), which are discussed in the following sections.

Fourier Analysis
Fourier analysis is a curve-fitting approach commonly used because of its well-known mathematical background (for reviews on the topic please refer to Lestrel, 1997;Claude, 2008). For outlines, direct application of Fourier transforms is not possible since they are defined as a function of x-and y-coordinates in two dimensions. There are two possibilities: • Express the outline as a function of one transformed variable, i.e., Polar transform (PT) of equally spaced radii or fast Fourier transform (FFT) of the tangent angle to the outline • Separately decompose the x-and y-coordinates and expresses these as functions of the curvilinear abscissa, i.e., elliptical Fourier analysis (EFA).
PT methods analyze the shape using a series of equally spaced radii originating from a barycenter of a study specimen (Younker and Ehrlich, 1977;Bird et al., 1986;Torres et al., 2000;Tracey et al., 2006;Claude, 2008). This method reduces the number of parameters well. However, because it is required that intervals are equally spaced, which may lead to more information being gathered in certain parts of the outline than others (Figure 2), PT methods are better suited to simple (and more rounded) outlines (Claude, 2008). In addition, as each radius can only intersect the perimeter once, information on particularly convoluted outlines that curve back on themselves can be lost (Figure 2; Rohlf and Archie, 1984;Tracey et al., 2006;Claude, 2008;Neves et al., 2011). As with all Fourier methods discussed here, PT is sensitive to placement of the starting position, which is often a biologically homogenous point. This can occasionally become a problem when no biologically homogenous point is available from which to start the trace (Haines and Crampton, 2000). Moreover, even though expression of the outline in polar coordinates could provide greater class separation distance, this effect is neglected when more complex statistical methods are applied, resulting in poorer classification .
Outline shape can also be studied using FFT, operating on the tangent angle as a function of arc-length connecting the coordinates (Zahn and Roskies, 1972;Doering and Ludwig, 1990;Ponton, 2006). FFT was improved by Haines and Crampton (2000) to ensure that reconstructed outlines are closed (i.e., determined by the dependence of the first harmonic on all others). This newer method also includes two coefficients per harmonic that describe both amplitude (size) and phase angle (angular offset relative to the starting position) so Fourier FIGURE 2 | Theoretical schematic of how data acquisition takes place on a cephalopod statolith using Polar transform (PT) of equally spaced radii. Data collection begins at a biologically homogenous location (0 • , i.e., the rostrum tip). The red line intersects the perimeter at two locations, suggesting this method would not be suitable for this structure. (a) Intervals describing a small amount of information, (b) region where a large amount of the shape is described by only one interval.
descriptors are computationally independent of each other, preventing spurious results (Haines and Crampton, 2000;Lord et al., 2012). Further improvement of the method was offered by Reig-Bolaño et al. (2010a), who expressed the contours of otoliths in combination of the upper half of the original image and the lower part of the reflected one. Regardless of the improvements, FFT estimation of coefficients appears to be more sensitive to pixel "noise" in the image compared to other methods such as PT (Rohlf and Archie, 1984). "Pixel noise" impact could be reduced by applying smoothing technique, where the weight of each contour point is combined of 50% of the weight of that point and 25% of each neighboring point (Haines and Crampton, 2000), or by "lasso contour" technique described by Harbitz and Albert (2015) where only the points that provide a non-concave contour between the succeeding points are analyzed. Summing up, this method may be used as an alternative to PT when a radius drawn from the center intersects the outline at more than one location (Ponton, 2006).
EFA is often considered the most powerful Fourier method to describe recording structures comprehensively (Mérigot et al., 2007;Ferguson et al., 2011;Leguá et al., 2013;Mahé et al., 2016) and is particularly efficient for reducing the number of variables of the original dataset. EFA works on raw x-and y-coordinates to delineate any type of shape with a closed twodimensional contour. EFA captures outline information in a quantifiable manner by decomposing a curve into a sum of several harmonically related ellipses (Kuhl and Giardina, 1982;Tracey et al., 2006). These ellipses can be combined to reconstruct an approximation of the outline shape (Figure 3). Each harmonic is composed of four coefficients, resulting from the projection of each point of the outline on the x and y axes. The higher the number of harmonics, the greater the complexity of the outline description (Kuhl and Giardina, 1982).
This method is advantageous in that it does not require equal intervals along the outline (like PT methods) and can, therefore, accommodate substantially more complex shapes, for example having data points more closely spaced on segments of high curvature (Crampton, 1995;Tracey et al., 2006). This method can also be applied to shapes where the radius intersects the perimeter twice (Crampton, 1995). However, this method has two disadvantages; firstly, the four Fourier coefficients per harmonic are not computationally independent, which may result in spurious correlations that could compromise statistical analysis (Haines and Crampton, 2000). The second issue is the relative down-weighing of all harmonics after the first. This may not be a problem for outlines that vary in fundamental ways but if a structure differs in more detailed aspects of the morphology, discriminatory power may be lost (Haines and Crampton, 2000).
Univariate characters are often collected as supplementary information in outline studies to calculate shape indices ( Table 1; Jemaa et al., 2015). Form factor is a means to estimate the surface area irregularity, with values of 1 when the outline is a perfect circle and < 1 when irregular (Tuset et al., 2003). These can give an additional insight into how shape varies between specimens (Tuset et al., 2003). The most common indices and their calculations can be found in Table 1.
Despite being considered the most powerful Fourier method, EFA was often underused compared to other Fourier methods (Crampton, 1995;Mérigot et al., 2007;Lord et al., 2012) but is rapidly becoming the preferred method of analysis.

Wavelet Transform Analysis
Whilst Fourier transform provides functions in the form of sines and cosines which are non-local and can, therefore, result in poor approximations of sharp edges (because they expand into a large number of high-order harmonics), WT uses approximating functions that are contained in finite domains (refer to Parisi-Baradad et al., 2005 for equations). This makes WT a robust and powerful alternative, well-suited for approximating sharp edges of complex outlines and transitions and establishing the location of these transitions (Graps, 1995;Parisi-Baradad et al., 2005Libungan and Pálsson, 2015). Mathematically speaking, WT is based on expanding a signal into a family of functions obtained as the dilations and translations of a unique function known as a mother wavelet (Mallat, 1991;Parisi-Baradad et al., 2005;Sadighzadeh et al., 2012). Studying the variability of coefficients at a given angle of an outline is not possible with Fourier methods, which do not provide information on localized differences but on the overall shape (Libungan and Pálsson, 2015;Tuset et al., 2019).
The rest of this chapter will review the existing studies which use GM techniques (predominantly outline-based methods) of FIGURE 3 | Theoretical schematic showing the elliptic Fourier approximation of a cephalopod statolith, using three harmonics (method based on Crampton, 1995).

OTOLITHS
Otoliths (or earstones) are paired calcified structures involved in mediating the senses of balance and hearing in teleost fishes (Gauldie and Nelson, 1990). Within the vestibular organs of fish there are three otolith pairs; the sagitta, lapillus, and asteriscus (Figure 4; Popper et al., 2005). Sagittae are routinely used for age, trace element and morphological analyses, and are composed of polycrystalline aragonite and other inorganic salts in a noncollagenous organic matrix (Campana, 1999;Campana and Thorrold, 2001). Lapillae are also formed of calcium carbonate in aragonite form whilst the asteriscus (and abnormal or crystallized Sagittae) are composed of vaterite (Volpedo and Vaz-dos-Santos, 2015). As composition and function differ between otolith pairs, care must be taken when deciding on which otolith pair to use for morphological analysis. Sagittae are often the largest otoliths, usually display the most diversity in shape and size between species (Campana and Casselman, 1993;Popper et al., 2005) and have a disproportionately large amount of literature available on their morphology compared to the lapillus (see review by Assis, 2005) or the asteriscus (see review by Assis, 2003).
Otoliths are known to grow throughout the entire lifetime of the animal (Campana and Thorrold, 2001). Crystallization processes are thought to be under the control of soluble organic molecules within the endolymphatic fluid surrounding the otolith (Campana, 1999;Campana and Thorrold, 2001). The substantial diversity in the shape and size of otoliths across species is striking and a clear indicator of this carefully controlled biomineralization process. Compared to birds and mammals, fishes possess a high diversity in inner ear morphology (Ladich and Schulz-Mirbach, 2016). Yet, the functional significance of these intricate and highly species-specific shapes is still unknown (Campana and Thorrold, 2001;Popper et al., 2005;Ladich and Schulz-Mirbach, 2016). It is speculated that the shape may have adaptive significance related to its physiological function in hearing and balance or adaption to various ecological conditions (Gauldie and Nelson, 1990;Ladich and Schulz-Mirbach, 2016). For example, fish from highly structured habitats (e.g., reefs or rocky bottoms) are generally considered to have larger otoliths than more migratory species in an open ocean environment that spend most of their time swimming at high speed (Reig-Bolaño et al., 2010b). Larger otoliths are also thought to be associated with groups considered specialists in sound production (sciaenids or haemulids) (Cruz and Lombarte, 2004).
Otolith shape is taxonomically discrete but within a population shape also varies. It is currently recognized that otolith shape is under dual regulation with overall shape regulated genetically (Gauldie and Nelson, 1990;Reichenbacher et al., 2009;Vignon and Morat, 2010), and a wide range of exogenous factors altering the rate of otolith growth, which in turn modifies otolith shape (Gauldie and Nelson, 1990;Smith, 1992;Begg and Brown, 2000;Cardinale et al., 2004;Vignon and Morat, 2010). This produces variation among conspecifics that have experienced contrasting life histories (Vignon, 2012;Tuset et al., 2019). In some cases, genetically regulated shape change has only been found to affect shape change in certain areas of the otolith. For example, Vignon and Morat (2010) found that genetically regulated shape change was prevalent in the rostrum and antirostrum regions of otoliths of coral reef snapper (Lutjanus kasmira), by comparing native and introduced fish. Otoliths are easy to locate and remove, can be stored effectively with minimal preparation and are frequently available as historical series for applications in age and growth analysis (Pawson and Jennings, 1996). Thanks to the pronounced differences in their shapes, even in the beginning of the XX century otoliths were routinely used for species identification at least to the genera level (Frost, 1924). Species-specific shape and size combined with their accessibility made otoliths amenable to morphological studies investigating taxonomy (Phillips, 1942;Bostanci et al., 2015), paleontology (Schwarzhans et al., 2017), diet (Frost, 1924;Pierce and Boyle, 1991), stock identity (Bergenius et al., 2006), ecomorphology (Gauldie, 1988;Lombarte and Fortuño, 1992;Volpedo and Fuchs, 2010), and evolution (Gaemers, 1984;Steurbaut and Nolf, 1989;Stransky and MacLellan, 2005).
Over the years, data collected on the characteristics of otolith shape and size of each species has resulted in the development of several otolith catalogs with detailed descriptions and illustrations of otoliths by region (Hecht, 1971;Morrow, 1979;Williams and McEldowney, 1990;Smale et al., 1995;Rivaton and Bourret, 1999;Assis, 2000;García-Godos Naveda, 2001;Leopold et al., 2001;Campana, 2004;Furlani et al., 2007;Svetocheva et al., 2007;Tuset et al., 2008;Sadighzadeh et al., 2012;Rossi-Wongtschowski et al., 2014;Volpedo et al., 2018). At the larger scale this led to development of the AFORO system which contains a database of than images of otoliths of more than 2300 species inhabiting the waters of the World Ocean (Lombarte et al., 2006a). Moreover, the system includes packages for automatic species identification based on the analysis of otolith outlines (wavelet and EFA). Latest advances in the system development include the proposal to use discrete Fourier transform of the outlines, which would allow classification of otolith images which are not optimally positioned (Marti-Puig et al., 2020).
A summary of otolith outline studies can be found in Table 2. This table provides a summary of 60 otolith outline studies but is not exhaustive given the large amount of literature available for this particular recording structure.
Despite the fact that EFA is considered to be used less frequently than other Fourier techniques in several publications (Crampton, 1995;Mérigot et al., 2007;Lord et al., 2012), this technique appears to be the most popular method of otolith outline analysis, with 61.7% of the publications in Table 2 involving EFA (in some cases more than one method of outlinebased analysis was applied). These publications were published within the last two decades, predominantly from 2010 onwards. This may be due in part to the development of free software allowing Fourier methods, predominantly EFA, to be applied with relative ease. FFT and WT methods were applied in 10 publications each, though FFT methods span three decades and WT methods were only applied to the listed publications within the last decade. Only 5 publications applied PT methods, spanning the last three decades. When the efficiency of EFA, FFT, and WT for species discrimination was investigated in 11 species of snapper, WT provided the best results (Sadighzadeh et al., 2012). Although EFA can represent any outline when an ample number of coefficients are used, in practice, the number of harmonics is limited to a lower value (Tracey et al., 2006), representing the original outline to a lesser degree than other types of analyses (Sadighzadeh et al., 2012). This can be an issue where there is a close similarity between different species or populations (Reig-Bolaño et al., 2010b). Ponton (2006) compared the utility of EFA, FFT, and landmark methods for species discrimination and found that despite the otoliths of four fish species having only a few, sparsely located, homologous landmarks, the landmark method appeared slightly more efficient at discrimination.
Whilst outline techniques are clearly a useful tool for identification of species and stocks, interpretation of results can be complicated with several co-variables resulting in the overall shape variability. Strong variability in shape has been related to sex, otolith position, age, allometry, and year-class (Simoneau et al., 2000;Cardinale et al., 2004).
Sexual dimorphism is often excluded as a co-variable, as different sexes frequently do not show significant differences in otolith shape (Bird et al., 1986;Castonguay et al., 1991;Simoneau et al., 2000;Cardinale et al., 2004;Paul et al., 2013). When the effect of sex is found to be significant, each sex is considered separately (Neves et al., 2011).
Left-and right-side sagittae usually have no consistent morphological differences in "round" fish (Begg and Brown, 2000;Cardinale et al., 2004). However, in flatfish significant asymmetry in weight, shape indices and Fourier descriptors have been identified between the left and right otoliths (Mille et al., 2015), which is supposed to be a consequence of their unique metamorphosis (Mérigot et al., 2007). Substantial variation in shape throughout ontogeny may also confound results. A study by Hüssy (2008) concluded that otolith shape development consists of an ontogenetic component, with otolith size and crenulation exclusively linked to fish size (see also Irgens et al., 2017), and an environmental component, with food consumption affecting the size and number of lobes formed. Higher food rations during early ontogeny lead to a higher number of larger lobes and a more rectangular otolith (Hüssy, 2008). If ontogenetic shape variation is not the focus of the analysis, it is suggested that only individuals from the same ontogenetic stage are selected or otolith shape is adjusted for with respect to allometric relationships with fish length. Castonguay et al. (1991) used PT methods to discriminate between mackerel stocks and expressed concern that discrimination may be a result of age-and year-class effects (samples differences) rather than real differences in otolith shape between groups. To mitigate the potential age effect, analysis can be restricted to individuals of a similar age, otherwise the effect of age should be investigated prior to the study of the shape (Bird et al., 1986;Begg and Brown, 2000). Analysis of only mature individuals is recommended for shape analysis, since Campana and Casselman (1993) proposed that otolith shape analysis may be of negligible value for stock discrimination when applied to immature fish, but it is important to consider what is appropriate for the species in question. In most cases, allometric growth changes in shape associated with age and growth are minimized by correcting for length using various standardization techniques or restricting analysis to individuals within a particular length range (Simoneau et al., 2000;Reig-Bolaño et al., 2010a;Harbitz and Albert, 2015).
Interannual differences in shape may confound the interpretation of results because variations in the environment or age structure of a population may vary between years (Vignon, 2012;Mahé et al., 2016). But frequently when multiple years are included there is temporal stability in the otolith shape (Campana and Casselman, 1993;Begg and Brown, 2000;Jónsdóttir et al., 2006;Denechaud et al., 2020). The use of multiple cohorts collected over several years was considered good practice by Bergenius et al. (2006) to derive a time-averaged assessment of the spatial structure of a stock. It has been recommended that discriminant functions should be recalculated for each of the major age classes to minimize the year class effect, which results from the variable growth conditions that each year-class is subjected to Begg and Brown (2000) and Cardinale et al. (2004), and to improve classification power (Paul et al., 2013).
Studies involving otolith shape have predominantly focused on commercial marine species (Simoneau et al., 2000). As the format of this review implies some limitations we will not review in details all studies of the otolith outlines. Instead we will focus on one iconic example of application of outline-based GM methods to stock units discrimination.
The commercially important cod (Gadus morhua) has an extensive geographical range on the continental shelves and banks of the North Atlantic Ocean. Different populations or stocks of cod, therefore, inhabit very different environments along their geographical range, having a marked effect on life history characteristics such as growth (Petursdottir et al., 2006). In addition, there is a diversity of migratory behavior exhibited, ranging from site-attached residents to wide-ranging dispersers (Galley et al., 2006), which mix to varying degrees. This makes stock assessment a very difficult task. But studies have demonstrated that otolith shape can be used to distinguish between some spawning groups of cod, suggesting that this technique may be a useful tool for stock assessment of this species.
One of the oldest studies on cod was undertaken in the Northwest Atlantic (Campana and Casselman, 1993). Though otolith shape was found to be attributable to several factors such as age, sex, year-class, or stock, shape was found to be a reliable indicator of stock identity. Results were largely consistent with other stock identification techniques including tagging, meristics, morphometrics, and parasite loads. In the Baltic Sea, two cod stocks have been identified (using a combination of meristics, morphometrics, genetics, and tagging experiments) which differ by spawning location and spawning period (Paul et al., 2013). Migrations of juveniles and adults result in the mixing of the two stocks, but none of the previous methods were able to assign individuals from mixed populations to their respective stock. In Paul et al. (2013), EFA was able to allocate otoliths to their correct stock with >90% accuracy, which suggests that EFA techniques can improve our understanding of mixing of stocks in the Baltic Sea. In Icelandic waters, cod have historically been assumed to belong to a single stock, but recent studies have shown that the stock may be composed of multiple units of local populations that have limited interactions (Jónsdóttir et al., 2006;Petursdottir et al., 2006). In Petursdottir et al. (2006), individuals sampled close to the coast differed to those samples further out on the bank and continental shelf in both growth and otolith shape. It was concluded that cod spawning in the coastal area may require special protection. Using PT methods, Jónsdóttir et al. (2006) identified two major groups of cod north and south of Iceland and some evidence of depth separation.
In the northern North Sea and west coast of Scotland, results of FFT outline analysis on spawning groups were consistent with genetic evidence of low gene flow, and tag-recapture evidence of resident populations and high spawning fidelity (Galley et al., 2006). The feasibility of outline analysis for stock identification has also been evaluated (and confirmed) for known-age Atlantic cod from the Faroe Islands, with total classification success between 79 and 85% between stocks for different age classes (Cardinale et al., 2004). Atlantic cod from coastal and offshore regions off northern Norway are routinely allocated to Norwegian coastal cod or Northeast Arctic cod based on internal morphological features of their otoliths (Stransky et al., 2008a).
EFA was found to be a cost-effective method of stock allocation that: (1) was quantitative and not subject to otolith reader-specific experience; and (2) did not require the otolith to be sectioned or broken, making the process less time-consuming and destructive.
Apart from the issues not specific to any of the recording structures, such as allometric alteration of the structure outline during ontogeny (Lleonart et al., 2000), otoliths introduce several challenges in the application of the GM analysis. Probably the most pronounced issue is the complexity of the otolith outline of some taxa (such as Scombridae, Dactylopteridae, Zeidae, Pleuronectidae). Most of the outline-based GM methods are unable to accurately reflect all the features of the otoliths' shapes (Harbitz and Albert, 2015). These issues could be overcome by applying WT method, or by transforming the image (Reig-Bolaño et al., 2010a;Sadighzadeh et al., 2012). However, in some cases minor details of the otoliths' outlines doesn't reflect population or stock specific traits but represent the result of individual variability (Harbitz and Albert, 2015). Thus, pecularities of the studied species biology should be considered in selection of the approach to the analysis of otoliths' outline.
A similar issue could occur when the stock is composed of several cohorts which experienced different environmental conditions during ontogeny (Mahé et al., 2016). There are several practices allowing to avoid these issues. First of all, fish age could be estimated prior to the shape analysis, this would allow comparison of the fish of the same age (in addition to clearer identification of the cohorts it would help to avoid the allometry related alteration of the outline). Alternatively, multiple cohorts collected during several years could be used to perform the analysis of the averaged stock (Bergenius et al., 2006). However, it should be mentioned that this approach is more time consuming and more demanding to the quality of the raw data.
Application of GM methods to otoliths has some advantages over other recoding structures as well. The most importantly different outline analysis approaches were applied to this structure much more than to the others. Thus, a lot of issues which in other structures remain a matter for future research were already addressed in otoliths. Specifically, a standard approach to image collection was developed (Lombarte et al., 2006a), different approaches were tested, and their benefits and limitations were evaluated (Sadighzadeh et al., 2012). Although, possibly the most important benefit over other recording structures is the development of AFORO database. The system provides complex solution to the otolith shape analysis for the species identification (Parisi-Baradad et al., 2010) and a huge database of images which could be used as a reference material.
Regarding the aims for the future research, the most promising fields appears to be the application of three-dimensional image acquisition and related analytical methods (Sadighzadeh et al., 2012). Another promising field is the development of software space integrating both landmark-based and different outlinebased shape description tools.

STATOLITHS
Analogous to fish otoliths, statoliths are paired calcareous concretions found within the statocysts which are responsible for the detection of linear and angular acceleration in aquatic invertebrates such as cephalopods, gastropods, bivalves, and cnidarians (Clarke, 1978;Arkhipkin and Bizikov, 2000;Arkhipkin, 2005). Statoliths of gastropods are composed of alternating layers of calcite and aragonite (Hollyman, 2017), while ones of cephalopods have a more complex structure where most of the statoliths body is composed of layers of calcite and aragonite, whilst the wing (Figure 5) is composed of vaterite (Arkhipkin and Shcherbich, 2012). Unlike statoliths in other animals, cnidarian statoliths are composed of bassanite (Tiemann et al., 2006). For jellyfish and squid, statolith morphometrics has shown promise in discriminating between taxa or stocks, but for other taxa such as gastropods, the statolith morphometry is likely to be too uniform and circular for any distinction (Hollyman, 2017).
Taxonomic identification of jellyfish using soft tissue characters is notoriously problematic as bell and tentacle tissue can be easily damaged (Mooney and Kingsford, 2016). But cubozoan jellyfishes possess one hard structure: a statolith.  (2016) used EFA to successfully discriminate among 12 separate species to family and in many cases species level. This suggests that statoliths may be a robust taxonomic tool when other diagnostic tools are damaged. This paper compared the EFA method (applied to proximal, oral, and lateral statolith faces) with traditional length to width ratios. EFA outperformed the length to width ratio, which did not account for statolith curvature, and the best discrimination was achieved when all faces were combined during EFA analysis. Subsequently, Mooney and Kingsford (2017) used EFA to discriminate among geographically distinct sampling populations of three species of cubozoan jellyfish. Successful discrimination was dependent on the species in question as well as the statolith orientation when the image was captured (proximal, oral, or lateral face).

Mooney and Kingsford
Compared to longer-lived and slower growing teleost fish, identifying squid stocks is challenging, as many have a short lifespan (and, thus, a rapid population turnover), are very mobile, have multiple cohorts, and have a high interannual variability in population size and biomass (Green et al., 2015). In addition, sympatric squid species that are difficult to distinguish in the field are common, making alternative means of identification necessary. The earliest study to analyze statolith shape variability using outline analysis used PT methods to investigate the intraspecific variation of three cephalopod species and the shape variation with growth (Lombarte et al., 1997). Since this initial study, there was a long period of time when outline analysis was not utilized for analysis of squid statolith shape, until 2015.
In recent years, it appears that outline analysis is gaining traction in the study of cephalopod stock structure and taxonomy. In a more holistic approach to stock identification, Green et al. (2015) combined FFT and statolith elemental composition to investigate stock separation in Nototodarus gouldi, which is currently managed as a single population. Statolith shape is considered a function of genetics and the summation of the environment an individual has experienced throughout ontogeny and elemental composition can be quantified at specific phases of an individual's ontogeny. It is considered beneficial to use these techniques' concurrently so stock structure can be investigated on two different time scales. Statolith shape provided evidence that adults caught in the two locations belonged to separate stocks, whilst elemental composition suggested hatching occurred throughout their distribution. Evidence of asymmetry in ontogenetic movement suggested that one location may need to be managed as the source stock (Green et al., 2015). This asymmetry would not have been evident had either method been used in isolation.
Outline analysis has also shown some promise for species identification in cephalopods. Statolith shape was compared in three ommastrephid squid species from the Pacific and the Southwest Atlantic Oceans using both traditional methods and EFA . The study concluded that use of the traditional method is effective when species are distantly related phylogenetically, and EFA is more appropriate at the genus level. Traditional methods are not likely to be more effective when individuals are in different phylogenetic groups per se, it is probably more appropriate to conclude that EFA found much higher within-species variation than traditional methods did. A similar study compared three loliginid squid with highly overlapping habitat ranges and similar body morphologies (Jin et al., 2017). Unlike ommastrephids, loliginid squid are a neritic family which undergo much less extensive diel-vertical and horizontal migrations (Jin et al., 2017) which may have an impact on shape variability. Statoliths and beaks were used for species identification using EFA techniques. It was concluded that a combination of statolith and beak was best for identification. Fang et al. (2018a) is the only study to date to use WT methods for statolith shape analysis. The results of the study supported previous conclusions that the Humboldt squid Dosidicus gigas has two stocks separated by hemispheres. As WT methods have the advantage of establishing locations of variance along the statolith outline, this study identified the rostrum and wing as the most variable regions among purported stocks. It is speculated that localized shape variations are partially induced by variable movement patterns of statoliths within the statocysts, which may be controlled by current characteristics or foraging behavior (Arkhipkin and Bizikov, 2000). Statoliths of pelagic and nearbottom squid may be distinguished by a number of differences in their shape (namely in the rostrum and wing regions), as the type of movement and, therefore, sensitivity to different forces of acceleration differ between the groups. In D. gigas, one stock was characterized by a wider, larger rostrum and a short narrow wing, characteristic of a near-bottom lifestyle which requires a higher sensitivity to low angular acceleration for "plane-like" movements near the seafloor, occasional settlement and for negotiating obstacles (Arkhipkin and Bizikov, 2000;Arkhipkin, 2003). Regardless of the origins of shape variation, it is clear that statoliths from different population components are often able to be distinguished using outline methods, but that this method is underused for this recording structure when compared to otoliths.
Thus, the potential of application of outline-based GM methods to study the shapes of statoliths highly depends on the taxa of interest. In the case of simple-shaped statoliths of cubozoan jellyfishes, EFA provides the best possible result, but in application to statoliths of cephalopods which tend to have a much more complex outline, this method has some limitations (Green et al., 2015;Jin et al., 2017). First of all, this method is limited in identifying the locations of variance (it could be done, by the "manual" comparison of the outlines representing the most deviant groups). Additionally, EFA is unable to reflect the minor details of the outline. However, the information value of these elements is questionable, as many of them represent individual differences or a result of statolith damage during extraction or processing.
The last is an issue of application of every outline-based GM method, especially when it comes to a statolith wing (Figure 5). Statolith wing is a very fragile structure and is often damaged during the capture of animals and extraction of statolith from a statocyst. Later, the analysis could be biased by samples with different degrees of damage on the wing. The most reliable approach to solving this issue is the removal of the wing area from the image as it was performed by Green et al. (2015).
Despite these limitations, outline-based GM analysis studies appear a powerful and time-/cost-efficient method of statolith discrimination in cephalopods and cubozoan jellyfishes. Promising fields for the future application of GM methods to statoliths research include the steps already performed in the studies of otoliths. Specifically, the classification power of different descriptive methods needs to be tested, their limitations and benefits should be evaluated. The worldwide database of cephalopod statoliths is a tool of great interest, especially for the diversity hotspots or the regions where they are intensively exploited.

BEAKS
Cephalopod mandibles are called beaks for their resemblance to the beaks of birds. However, unlike the latter, in cephalopods this term is applied to a single mandible (Clarke, 1962a). Being one of the hardest parts of the cephalopod body, beaks serve predominantly for biting and drilling prey. The beak's rostrum tip represents one of the hardest organic structures, while lateral wall and wing areas, which are connected to the muscular buccal mass, are very soft and flexible (Miserez et al., 2008). These traits are determined by the chemical composition of the beaks. They are composed of a protein-chitinous complex, with the proportions of water, protein, and chitin differing greatly from the rostral tip to the soft parts of wings or lateral walls (Miserez et al., 2008).
Due to their function, the application of beak outlines in research holds certain challenges when compared to other recording structures. In particular, the "working" surface of the beak is subject to erosion, whilst soft parts of the beak undergo significant deformations if it is preserved dry or in alcohol even for a short time (Figure 6). In addition, beaks are relatively large structures that have surfaces covered by a number of protrusions, ridges, folds, and depressions. Finally, beaks also have uneven pigmentation, which provides an additional layer of information that is often taken into consideration in taxonomy (Figure 6). At present, the complexity of the beak's shape and its susceptibility to both erosion and deformation limit the use of complex outline analysis methods. To define the preconditions and prospects of the application of the GM techniques in beak shape analysis, it seems worthy to mention several major stages of beak use in cephalopod studies.
The first descriptions of the beaks as part of the species description were rather simple and far from the modern standards of morphological descriptions (Kent, 1874;Verrill, 1881). Beak shape description was further developed in the second half of the twentieth century (Clarke, 1962a,b). During this period, it was demonstrated that beaks could serve as a powerful tool for identification of cephalopods in the stomachs of predators, not only at the genus or species level but at the population or stock level as well (Clarke, 1962b). Twenty years later, beak shape descriptions were recognized as an essential tool for the description of new cephalopod taxa (Roper and Voss, 1983). Verbal descriptions and measurements were included in the list of mandatory procedures for species description (Robinson and Hartwick, 1983;Clarke, 1986).
Nowadays, quantifying beak shape is routinely used as a tool for cephalopod identification, particular in the stomach contents of predators (Willassen, 1986;Smale et al., 1993;Lu and Ickeringill, 2002;Xavier et al., 2007;Xavier and Cherel, 2009). In the description of a newly discovered octopus species, beak shape descriptions were used as an argument to justify the assignment of a species to a specific genus (Allcock et al., 2007). The potential for using beak morphometry for stock and population discrimination has been examined in several studies. Analysis of measurement indices showed that beaks of Loligo forbesii are significantly different between squids caught in the Azorean waters and those caught in the European continental shelf area (Pierce et al., 1994). Differences in beak morphometry among regions of their species range were also shown for the ommastrephid squids Illex coindetii and Illex illecebrosus (Martínez et al., 2002). Finally, differences in beak morphometry were used as a supplementary tool in the separation of the Atlantic and Pacific populations of Doryteuthis gahi (Vega et al., 2002).
With the proliferation of beak-related morphological studies, several limitations in the traditional approaches to shape description were identified. The study of beak shape of octopuses in the Southern Ocean showed that whilst the resolution of measurements and verbal descriptions allowed for discrimination between genera, it was insufficient for reliable species determination or for phylogenetic purposes (Ogden et al., 1998). As well as traditional methods not being able to reliably discriminate between lower taxonomic levels, they have also been shown to produce conflicting results. Research within the genus Sepia indicated that the dorsal and ventral beaks followed different interspecific patterns in shape and that the dorsal beak represented a more reliable source of information (Neige, 2006). However, for the genus Idiosepius, beak shape analysis did not reveal any significant interspecific differences in the appearance of the dorsal beaks among the representatives of the genus (von Byern and Klepal, 2010). Neither verbal descriptions nor traditional linear measurements can apparently convey all the structural features of beaks for a particular species or population. Luckily, the solution to this issue could be found in the field of geometric morphometrics.
As mentioned above, the application of GM-based methods to the beaks of cephalopods is hampered by their susceptibility to erosion and deformation, the complexity of their shape and their large sizes. So, the first issue to solve was the way to represent a three-dimensional structure with a complex shape in a two-dimensional space. The first attempts to display beaks as three-dimensional objects used specific cameras equipped with a system of mirrors (Crespi-Abril et al., 2010). However, the use of such a camera only obtained images of beaks from the lateral and frontal sides. Other authors studied "simplified" beak shapes. In this case, only the lateral profile of the beak was subject to analysis (Fang et al., , 2018bJin et al., 2017). Analyzing only the lateral beak profile led to a loss of information, as some distinguishing traits are visible only when examining beaks in the dorso-ventral plane (Smale et al., 1993).
On a par with applying different approaches to capture images, different methods of outline analysis itself were tested. At present, the most commonly used method is the landmark analysis (Crespi-Abril et al., 2010;Fang et al., , 2018b. However, as this approach is out of the scope of the present review, this will only be briefly mentioned. In general, landmarkbased GM results are rather inconsistent. The study of Illex argentinus beaks revealed no differences between sex, between immature and mature squids, or between spatial groupings (Crespi-Abril et al., 2010). However, the same method revealed significant differences between western and eastern stocks of Ommastrephes bartramii . The study of the beak outlines of four abundant species of octopuses from the East China Sea demonstrated twofold superiority of the landmark-based analysis in the accuracy of animal identification in comparison with traditional measurements (Fang et al., 2018b). More importantly, this study demonstrated the advantage of machine learning methods in the shape analysis data, in some cases it led to an increase in the accuracy of a species determination from 31.4 to 77.3%. The most recently, landmark-based outline analysis of cephalopod beaks was applied to test the differences between Eledone cirrhosa experiencing developmental malformations (bearing additional buccal pieces at the location of the first sucker of ventral arms) and ones with the normal structure of the raptorial apparatus (Fernández-Álvarez et al., 2020). Analysis revealed significant differences between these individuals, although authors assume that identified differences don't affect octopus's trophic position (which resulted in occupation of the same isotopic space) or survival.
To date, EFA has been applied to the study of the beak outlines only once. The study of the beak shape patterns of three loliginid squids from the South China Sea demonstrated comparable levels of interspecific discrimination (Jin et al., 2017). However, this study also showed that a combination of two recording structures (statoliths and beaks) is more reliable for species identification. Thus, geometric outline analysis of cephalopod beaks is a promising method for species and stock determination, but some problems obstructing the application of these methods need to be resolved. First of all, it seems essential to develop a method that would allow images of beaks from the lateral, frontal and dorsal sides to be collected. The method used by Crespi-Abril et al. (2010) seems to be limited in efficiency due to the complexity of the positioning of the camera and mirrors and the method does not collect images of all three sides of the beak. Methods used by other authors (Fang et al., , 2018bJin et al., 2017) can be modified, so images of the dorsal surface of the beak are obtained first, and then beaks are processed to obtain images of the lateral side. Regardless, following both approaches stated above, every side of the beak is still analyzed as a separate structure. The alternative option is the application of scanning devices or 3D modeling software, but this approach has some limitations as well. First of all, it's much more time consuming than the processing of standard two-dimensional images. Moreover, careless application of this approach could lead to a significant loss in model accuracy, which is hardly traceable by an inexperienced observer. However, none of these approaches avoid the issues caused by deformation of the soft parts of the beak. A possible solution is inspired by Robinson and Hartwick (1983), who found that measurements of pigmented elements of the beak are a more reliable source of information than measurements that included both pigmented and non-pigmented parts of beaks. Similarly, analysis of the pigmented part of the beak only would exclude any inaccuracy caused by the deformation of part of the beak. Moreover, the combination of the outline analysis with the coloration analysis would lead to a further increase in the accuracy of the analysis . The most promising way to increase the efficiency of the analysis itself is the application of machine learning methods, as shown by Fang et al. (2018b).

SHELLS
Shells are complex structures that have a similar function to vertebrae in mammals. The main function of the shell, which is characteristic of all classes of the phylum Mollusca, is the support of the muscles and other soft tissues (Bizikov, 2008). Another equally important function is the protection of the animal, which is more characteristic of mollusks with external shells, such as bivalves and gastropods. The shells of most mollusks are calcified structures formed by epithelial cells of the mantle (Marin and Luquet, 2005). However, in the class Cephalopoda, chitinous shells are more common (Bizikov, 2008).
Since mollusk shells are extremely diverse in terms of composition and their location in relation to soft tissues, the limitations on their use for GM analysis varies greatly between taxa (Figure 7). For example, cephalopod vestigial shells have similar challenges for GM-based studies as beaks in that if they are stored dry, or in alcohol, they undergo significant deformation, making analysis impossible. However, gastropods and bivalves do not have issues with deformation but the outer shells of some specimens can be so densely covered with fouling that the GM analysis becomes impossible.
Largely due to the specifics of the taxa themselves, the shells of different groups of mollusks have been studied extremely unevenly. For example, in studies of the modern cephalopod shells, GM-based analysis has only been applied once (López-Galán and del Carmen Alejo-Plata, 2015). In this research, traditional morphometrics and landmark-based GM analysis were applied to test gender differences in Octopus hubbsorum stylets. Landmark-based methods revealed differences between males and females, but the accuracy of sex determination was as low, as was the case with traditional morphometrics. It should be noted that the sample size in this study was rather small (160 individuals) and that stylets (the internalized shell of octopus) are subject to pronounced variations. Nevertheless, as these results are generally compatible with the attempts to find the genderspecific differences in other recording structures of cephalopods (Lombarte et al., 2006b;Crespi-Abril et al., 2010), this approach is probably not accurate enough, as concluded by the authors.
GM-based methods are more prevalent in paleontological studies of extinct cephalopods, such as ammonites. The main application in this field is the determination of interspecific differences (Reyment and Kennedy, 1998;Reyment, 2003;Gerber et al., 2007), but data can also be used to determine phylogenetic relationships (Reyment, 2003) or to describe ontogenetic variability (Gerber et al., 2007). In all the cases mentioned, the approach to the analysis of the shell appearance was rather conservative. Reyment and Kennedy (1998) used landmarkbased relative warp analysis to test for differences in shell shape of three ammonite species of the genus Neogastroplites. The authors indicated that use of counts of the shell ribs as supplementary information increases identification accuracy. Later, the same methods were applied to test the hypothesis that the Nigerian fauna of ammonites was highly connected and at some point even derived from Northern African fauna (Reyment, 2003). This hypothesis was confirmed with the study showing that despite high intraspecific plasticity, there is clear morphologic evidence of connectivity between these groups of ammonites. Finally, Gerber et al. (2007) addressed intraspecific differences in juvenile and adult ammonites of the family Hildoceratidae using landmark-based methods. Thus even though geometric morphometric analysis never was used in stock discrimination of cephalopods its application in the paleontological studies shows that this approach could have a great potential in this field.
Except for a couple of studies (Sobrepeña and Demayo, 2014;Doyle et al., 2018;Tuset et al., 2020), the main approach to the GM analysis of gastropods shells is a landmark-based analysis. It is noteworthy that in both cases of using EFA, high classification accuracy of mollusks was achieved. In the study of the shells of Bolinus brandaris bearing highly pronounced spines (Tuset et al., 2020) potential of EFA and WT analysis to accurately describe their shells and provide correct classification of these molluscs from spatially different populations was tested. The study showed that WT analysis provides slightly better classification than EFA, 82.7 vs. 75.3% of correct classification, respectively. The former allows more accurate description of objects with irregular shape (like the spiny shells of Bolinus brandaris) than the latter. Doyle et al. (2018) compared the classification power of landmarkbased analysis and EFA and showed that the latter provides higher accuracy in the classification of cases. Nevertheless, the landmark-based analysis, which is used routinely in studies of the intraspecific variability of gastropod shells also provides a reasonable level of classification accuracy (Becker et al., 2016;Cazenave and Zanatta, 2016). Despite the fact that landmarkbased GM has been used in studies of gastropod mollusks for about 30 years (Johnston et al., 1991), this method truly flourished only during the last decade, when the development of imaging technologies improved the quality and resolution of digital images, allowing for a sufficient number of landmarks to be distinguished in the images. In a number of studies (Moneva et al., 2012b;Sepúlveda and Ibáñez, 2012;Zhao et al., 2019), environmental factors determining such variability were analyzed. Zhao et al. (2019) indicated that temperature was the key environmental factor determining the variability of the Monodonta labio shell. Other studies (Moneva et al., 2012b;Sepúlveda and Ibáñez, 2012) concluded that the shape of the shell is determined by a combination of factors, including temperature, chemical characteristics of the water, and the impact of currents. Among the key factors determining the appearance of the shell, different authors identified physical factors such as temperature (Amini-Yekta et al., 2019) and power of currents (Vuolo et al., 2011), chemical factors such as waste from the transport and oil-refining industries (Márquez et al., 2011;Abdelhady, 2016;Primost et al., 2016), and ecological factors such as parasite infestation (Vasallo et al., 2013. The role of heredity in the formation of a particular shell shape is not uniquely determined (Madeira et al., 2012;Dillon and Jacquemin, 2015). According to a study of the heritability of shell appearance in Physa acuta and Physa carolinae, genetic line affiliation accounts for up to 68.5% of the variability. However, it should be mentioned that these results were obtained in laboratory conditions, where environmental properties are less variable than in the wild. On the other hand, the study of the geographical variation of Cerithidea decollata (Madeira et al., 2012) showed that the genetic variation and mollusk shell shape variation can be independent of one another. The study suggests that genetic characters are determined by the distribution history of the species, whilst the shell structure to a greater extent reflects the influence of local environmental factors. This probably explains the fact that in some studies, GM-based analysis and molecular analysis showed similar results (Cruz et al., 2012;Vaux et al., 2018), whilst in others they did not (Becker et al., 2016).
Finally, when discussing the application of GM-based methods to the gastropod shells, the existence of critical comments on the modern approach should be addressed. Malvé et al. (2018) noted that despite the unquestionable advantages, the use of such methods is laborious and timeconsuming. This obstructs its application to the examination of large collections, where traditional morphometrics may be much easier to use and the accuracy is sufficient to detect intraspecific variability. A comparison of different GM-based methods to describe the variability of Littorina littorea shell shapes (Doyle et al., 2018) showed that the EFA descriptors provided higher accuracy in the classification of mollusks. In this case, higher efficiency of EFA could be explained by the fact that the Littorina littorea shell has a relatively simple sculpture (i.e., it bears no spikes and ridges). Larsson et al. (2020) noted that traditional morphometric approaches are unable to convey a significant part of the information actually "recorded" in the shape of shells. Finally, research on the impact of environmental factors on the color patternation of Littorina saxatilis is worthy of mention (Sokolova and Berger, 2000). The study found that color can not only be used as a source of supplementary information in the analysis of shell shape, but can also serve as a reliable tool for intraspecific group identification.
Similarly to gastropod devoted studies, the main fields of application of GM-based studies of bivalve mollusks are: identification of intraspecific groups (Costa et al., 2008(Costa et al., , 2010Funk and Reckendorfer, 2008;Márquez et al., 2010;Rufino et al., 2013;Boretto et al., 2014), assessment of the impact of environmental factors (Innes and Bates, 1999;Krapivka et al., 2007;Funk and Reckendorfer, 2008;Márquez et al., 2010;Boretto et al., 2014;Telesca et al., 2018) and investigating ontogenetic and evolutionary changes (Serb et al., 2011;Boretto et al., 2014;Alvarez and Pérez, 2016). Funk and Reckendorfer (2008) addressed the morphological variability of Pisidium subtruncatum shells. EFA indicated that shells are highly variable in size and shape and that high intrapopulation variability was observed in the sites with higher environmental variability. It was found that the amount of sediment is the key shapedetermining parameter. The connection between the shell outline and hydrological conditions (such as presence and power of currents) was also shown. Costa et al. (2010) took digital images of the left valve in mollusks of the genus Tapes, performed an EFA and analyzed the harmonic coefficients using multivariate classification. Results showed that the contour accurately reflected the general shape of the shell and highlighted that structural traits that are not distinguishable by the naked eye. Correct classification was achieved using partial least square discriminant analysis in over 95% of cases, which showed that EFA can be used as an identification tool. An analysis of stock variability of the striped clam Ameghinomya antigua in the San José Gulf showed that every specific fishing site corresponded to a specific shell shape (Márquez et al., 2010). The use of a landmark-based analysis and EFA allowed for correct classification in up to 90% of specimens. A similar approach was applied to the pod razor shell Ensis silique in the waters of the Eastern Atlantic (Rufino et al., 2013). In the framework of this research, shell shapes were described using a combination EFA of the outer surface of the shell and landmarks located on the inner surface. Both approaches demonstrated that molluscan shells from Central Portugal differed from those from neighboring regions, but were similar to the shells from Irish waters. Another study on the shell shapes of Tapes decussatus and T. philippinarum was devoted to analyzing factors affecting shell outline variability (Costa et al., 2008). The images of the left shell valves were studied using EFA and partial least square discriminant analysis. The study revealed that morphological distances between two species were higher than between populations of the same species. Finally, the temporal variation of the striped clam shell was analyzed by Boretto et al. (2014). Using landmark-based methods and EFA they investigated the changes of the striped clam shell shape from Pleistocene to the present day. Both approaches clearly showed significant differences between mid-Pleistocene and modern shells. Authors suggested that these differences are likely a result of habitat change.
Environmental impact on shell shape variability in bivalve molluscs was studied on the mussel Mytilus chilensis (Krapivka et al., 2007). Images of the right shell valves were digitized using EFA and analyzed by canonical variate analysis. The study showed that, though there was an absence of genetic differences, mussel populations had distinct differences in shell shape, most likely determined by environmental factors such as temperature, salinity, and wave exposure. This issue was further studied using two species of mussels (Mytilus edulis and M. trossulus) from the Eastern Atlantic Ocean (Telesca et al., 2018). Images of mussel shells were expressed using EFA and analyzed by principal component analysis. The main principal components were used to model the environmental factors influencing shell shape variability. The study showed that harsh living conditions (such as low temperatures and poor food availability) had a similar effect on the exterior appearance of shells in both species at all geographic sites examined, and that the shape of mussel shells can serve as an identifier of adverse conditions for mollusks. The influence of environmental factors on mussel shell shape was also investigated by Innes and Bates (1999). Unlike Telesca et al. (2018), these researchers attempted to find out whether the appearance of shells is determined by species affiliation or ontogenetic changes. Shell images were studied using EFA and canonical variate analysis. Mussel shell shape was found to differ between species, sampling site and size group. The greatest differences were found between spatial groupings, which lead to the suggestion that mussel shell shape is more determined by environmental influences than by genetics or ontogenetic stage.
Molluscan shell outlines of long-distance swimming scallops were examined using landmark-based GM methods (Serb et al., 2011). Shells of such long-distance swimmers were found to have very similar and non-variable shape patterns. Unlike long-distance swimmers, byssal attached, freeliving epifaunal and even relatively active Pseudamussium septemradiatus scallops had high variability in shell shape, both at the inter-and intraspecific levels. Shells of two bivalve species, Retrotapes exalbidus and R. antarcticus (Alvarez and Pérez, 2016) were examined to check if the morphs of the exterior shells presented in both species are environmentally determined, or whether shells were undergoing ontogenetic modifications. EFA showed that the globoid shape of the shells does not reflect any species or population features but is a typical characteristic for gerontic representatives in both species.
Summarizing, application of GM-based methods in studies of molluscan shells is as uneven (in terms of methods used and the results achieved) as the shells themselves. The most commonly used method is landmark-based analysis. However, methods such as EFA represent a promising tool for studies of gastropod and bivalve shells. In the framework of the cephalopod shell studies, application of GM-based methods is still quite limited and comes down mainly to the characterization of interspecific differences. The positive results achieved in paleontological studies highlight the potential of using such an approach, at least in the field of phylogenetic studies. For research on gastropod shell shape, GM-based methods represent a powerful tool for separation of intraspecific groups and stock units and for the determination of factors responsible for a specific shape. In this regard, studies on the effect of pollution are of particular interest, since, according to their results, the shape of shells can serve as an indicator of contamination of individual stock units in concentrations that are dangerous to humans. Caution should be taken as the use of GM-based methods has some limitations, both in comparison with classical morphology (Sokolova and Berger, 2000;Malvé et al., 2018) and with more advanced approaches (Larsson et al., 2020). In the studies of bivalve shells, geometric morphometrics are a routine tool for the shape analysis. Due to the fact that, unlike other recording structures, shells of the taxon usually have rounded outlines, EFA appears to be the most common and convenient method of analysis (Innes and Bates, 1999;Krapivka et al., 2007;Funk and Reckendorfer, 2008;Márquez et al., 2010;Rufino et al., 2013;Boretto et al., 2014;Alvarez and Pérez, 2016;Telesca et al., 2018). In cases where the use of EFA is limited, application of landmarks can provide accurate discrimination (Márquez et al., 2010;Rufino et al., 2013;Boretto et al., 2014).

CONCLUSION AND FUTURE DIRECTIONS
A holistic approach toward stock identification is encouraged to improve the accuracy of stock discrimination and provide greater insight into stock structure and the mechanisms that underlie it (Begg and Waldman, 1999;Vignon and Morat, 2010;Ferguson et al., 2011;Baldwin et al., 2012). A holistic approach uses a broad spectrum of complimentary techniques to provide information on stock structure at different temporal scales. From this review, it is clear that geometric outline analysis is a powerful tool which has the potential to complement other population markers and improve stock estimates. Other population markers may include genetics, body morphometrics and meristics, artificial tags, trace elements, stable isotopes, parasite markers, and life history parameters (Pawson and Jennings, 1996;Baldwin et al., 2012;Brophy et al., 2016). Some investigations have already combined outline analysis and another method into a single study, which appears to be a useful exercise (Ferguson et al., 2011;Green et al., 2015).
The advantage of using outline analysis in population or taxonomic identification is that this method is cost effective and only requires images from which outlines can be extracted (Libungan and Pálsson, 2015), which is less time-consuming than other techniques. Burke et al. (2008a) explained that herring (Clupea harengus) populations are often differentiated with microstructure analysis of the otolith, which requires time consuming polishing and sectioning prior to analysis. Outline analysis provided a fast and reliable alternative, successfully separating the two components without the need for substantial destructive pre-processing steps. It must be noted that careful consideration must be made prior to analysis on which type of outline analysis is required in relation to the general shape of the recording structure in question.

AUTHOR CONTRIBUTIONS
JJ was responsible for developing the general structure of the manuscript, writing sections "Introduction" to "Statoliths" and "Conclusion and Future Directions, " and drawing Figures 2, 3. FL was responsible for developing the general structure of the manuscript, writing sections "Introduction" to "Shells, " and drawing Figures 1, 4-7. Both authors were involved in the development of the reviewed version of the manuscript and contributed to the development of this manuscript.