Unravelling the phylogenetic and ecological drivers of beak shape variability in cephalopods

Cephalopod beaks are essential for prey acquisition and fragmentation during feeding. Thus, it is expected that ecological pressures affect cephalopod beak shape. From a practical perspective, these structures are also used to identify gut contents of marine megafauna, such as toothed whales, sharks, seabirds, and large pelagic fishes. Here, we investigated the relative importance of ecological pressures and phylogenetic relatedness in the evolution of beak shape using a wide range of Mediterranean cephalopod species. Phylogenetic analyses based on complete mitogenomes and nuclear ribosomal genes provided a well-supported phylogeny among the 18 included cephalopods. Geometric morphometric and stable isotope methods were implemented to describe interspecific beak shape and trophic niche variability, respectively. Phylogenetic signal was detected in the shape of both parts of the beak (upper and lower). However, lower beak shape was more distinct among closely related species, in line with the empirical notion that lower beak morphology is more useful as an identification tool in cephalopods. Interestingly, no association between beak shape and trophic niche (stable isotope values) was found. These results suggest that the evolution of cephalopod beak shape as quantified here is mainly driven by phylogenetic relationships, while feeding habits play a minor role.

of beak shape using a wide range of Mediterranean cephalopod species. Phylogenetic analyses based on complete mitogenomes and nuclear ribosomal genes provided a well-supported phylogeny among the 18 included cephalopods. Geometric morphometric and stable isotope methods were implemented to describe interspecific beak shape and trophic niche variability, respectively. Phylogenetic signal was detected in the shape of both parts of the beak (upper and lower). However, lower beak shape was more distinct among closely related species, in line with the empirical notion that lower beak morphology is more useful as an identification tool in cephalopods. Interestingly, no association between beak shape and trophic niche (stable isotope values) was found. These results suggest that the evolution of cephalopod beak shape as quantified here is mainly driven by phylogenetic relationships, while feeding habits play a minor role.
Abstract Cephalopod beaks are essential for prey acquisition and fragmentation during feeding. Thus, it is expected that ecological pressures affect cephalopod beak shape. From a practical perspective, these structures are also used to identify gut contents of marine megafauna, such as toothed whales, sharks, seabirds, and large pelagic fishes. Here, we investigated the relative importance of ecological pressures and phylogenetic relatedness in the evolution

Introduction
Studying the morphology of animal structures and their functionality is important for understanding the evolutionary trends that generate biodiversity (Williams 1972;Karr and James 1975). This assumes a given morphology is fitted for a certain function, and its functionality may act as an intermediate link between morphology and ecology (Arnold 1983;Ricklefs and Miles 1994;Collar and Wainwright 2006). Although it may seem that a certain morphology results in a set of ecological consequences, these biological processes are not directly related to morphology, rather they are related to the organism's capacity to perform a certain task (Losos 1990). Thus, performance, not morphology, expressed though the traits of organisms, is selected by natural selection (Wainwright 2007;Ibáñez et al. 2021a, b) and the outcome of this selection is observed as morphological variation across different organisms (Fernández-Álvarez et al. 2020;Masello et al. 2022).
Biological structures and their morphology can be subject to stronger or weaker selective pressures depending on their importance for fitness, the ultimate performance trait (Arnold 1983;Wainwright 2007). The link between morphology and functionality is generally expected to be tighter for structures that accomplish a vital biological function and have a significant impact on fitness (Collar and Wainwright 2006). For example, head morphology associated with bite performance is a major model to investigate the integration of morphology and ecology through performance, and this has been studied in several taxa, from vertebrates such as mammals (Santana et al. 2010), birds (Herrel et al. 2005), lizards (Kaliontzopoulou et al. 2012), or fishes (Herrel et al. 2002), to invertebrates (Püffel et al. 2021). These studies represent different models that have been explored exhaustively and which revealed the high importance of this structure for catching prey, ensuring reproduction or establishing dominance among conspecifics. All these models include crushing organs with hardened and articulated structures, usually called mandibles, which are basic tools for food acquisition and processing (Turnbull 1970;Boyle and Rodhouse 2005). These structures are commonly highly plastic in shape, and valuable to study given their wide ecological and phylogenetic variation (Renaud and Auffray 2010;Booher et al. 2021).
In marine ecosystems, cephalopods are a group of marine invertebrates of high interest to biologists as they are highly diverse in both form and life strategy (Boyle and Rodhouse 2005). They also have a hardened and articulated crushing structure, known as the beak. Cephalopods are ecologically important in marine food-webs worldwide, as they occupy a wide range of ecological niches Jereb et al. 2014), and constitute a major component of marine biomass (Boyle 1996;Smale 1996). Thus, they act as key species for the bottom-up transfer of energy within marine ecosystems (Boyle and Rodhouse 2005). Cephalopods are fast-moving semelparous predators with short lifespans (generally < 2 years; Boyle and Rodhouse 2005), and high metabolic and growth rates (Rodhouse and Nigmatullin 1996). They prey voraciously on a high diversity of organisms, which varies depending on the cephalopod species, habitat or ontogenetic stage (Rodhouse and Nigmatullin 1996;Boyle and Rodhouse 2005;Navarro et al. 2013;Villanueva et al. 2017). Their beaks consist of two chitinized parts that act together for prey capture and subjection, injection of venom and of digestive fluids, food fragmentation and ingestion (Boyle and Rodhouse 2005). Beaks vary enormously in shape across species, and they can completely change in proportions and size depending on the species and ontogeny (Xavier and Cherel 2009;Franco-Santos and Vidal 2014;Fang et al. 2018). It is thought that their wide variability in shape is a reflection of the wide variability in cephalopod feeding habits (Boyle and Rodhouse 2005). Some authors suggested that beak shape may exhibit specific adaptations for prey consumption (Roscian et al., 2022), at least in certain ontogenetic stages (Franco-Santos and Vidal 2014). Furthermore, due to the interspecific cephalopod beak shape variability and to its resistance to digestive fluids (Boyle and Rodhouse 2005), beaks are used as an identification tool of cephalopod remains in gut contents of fishes, their direct competitors (Packard 1972;Smale 1996) as well as marine megafauna, including species such as toothed whales, sharks and seabirds (Clarke 1986;Xavier and Cherel 2009;Tan et al. 2021). As such, although previous studies have suggested that the beak of cephalopods is under certain ecological selective pressures, it is reasonable to assume that it is also a structure that reflects the phylogenetic relatedness between species (Clarke 1986;Clarke and Maddock 1988;Xavier and Cherel 2009;Tanabe et al. 2015). This suggests there could be critical information stored in beak shape that could improve our understanding of cephalopods, yet beak shape variation has never previously been explored in the context of combining both phylogenetic relatedness and trophic niche based on stable isotopic analyses.
In the present study, we investigated the relative importance of ecological mechanisms and phylogenetic relatedness for determining the morphological diversity of cephalopods beaks. For this, we examined the beak shape variation (geometric morphometric methodology), phylogenetic (mitogenomes and nuclear ribosomal data) and trophic niche (stable isotope values) of 18 cephalopods sampled along the western Mediterranean Sea. Using multivariate analyses that combined phylogeny, morphology and trophic values, we were able to quantify and test the relative influence of ecological factors and phylogenetic relationships on cephalopod beak morphological diversity.

Data collection
In total, 214 individuals of 18 cephalopod species were collected along the North-western Mediterranean Sea (Fig. 1) between the years 2013 and 2019, representing four different orders [Octopoda (octopus), Oegopsida (oceanic squids), Sepiolida (bobtail squids) and Sepiida (cuttlefishes)], including species from every marine habitat occupied by cephalopods in the study area (from 0 to 2000 m depth, Table 1). To avoid potential ontogenetic effects on morphological and trophic niche analyses, only mature individuals were included in the analysis (maturity stages were identified following Cuccu et al. 2013 for octopus and Lipiński and Underhill 1995 for the other groups). No mature individuals of Myopsida from this area were retrieved. For each individual, the dorsal mantle length (DML ± 1 mm) and body mass (± 0.1 g for individuals of less than 50 g and ± 1 g for heavier ones) were measured. Also, the two parts of the beak (Upper Beak-UB, and Lower Beak-LB; see Fig. 2) were extracted for subsequent morphometric and isotopic analyses. A total of 202 UBs and 205 LBs were successfully extracted from the buccal mass and preserved in vials with 70% ethanol prior to further analyses. A subsample of 116 UBs was used for the isotopic analyses (Table 1).

Phylogenetic analyses
Phylogenetic analyses were performed based on a phylogenomic database (GenBank Accession numbers in Table 2). New data were produced through Genome Skimming (Dodsworth 2015), a shallow whole genome sequencing method that allows large regions of the genome, such as the complete mitogenome and the complete nuclear ribosomal cluster, to be obtained. DNA was extracted using the Purelink genomic DNA Mini kit (Invitrogen, MA, US) following the manufacturer instructions. Indexed libraries were prepared using a BGI Library Kit and sequenced 9 Gb/sample in an DNBseq-G400 (Beijing Genomics Institute, Shenzhen, China). The quality of the reads was assessed through FastQC (Andrews 2010). Mitochondrial and nuclear ribosomal DNA were assembled de novo using NOVOPlasty3.8.3 (Dierckxsens et al. 2016) using a reference sequence (either the complete mitogenome or the complete nuclear ribosomal gene cluster of a closely related species), and a fragment of cytochrome c oxidase subunit I (cox1), 12S rRNA or 16S rRNA (for the mitogenomes) or a fragment of 18S or 28S rRNA (for the nuclear markers) as a seed. For mitogenome gene annotations Mitos2 (Bernt et al. 2013) was used, with NCBI Ref-Seq 63 Metazoa database reference and genetic code 5, for invertebrates. Gene annotations were checked and corrected by hand. Because of the presence of duplicate genes in the mitochondrial genome of oceanic squids, NOVOPlasty did not return circularized genomes (see Fernández-Álvarez et al. 2022). For solving this methodological artifact, we used the mitogenome gene orders established for the flying squids Todarodes pacificus and Watasenia scintillans by Yokobori et al. (2004) using long PCRs. Nuclear 18S and 28S rRNA were annotated using RNAmmer (Lagesen et al. 2007).
Individual genes were concatenated and a database was constructed using the 13 protein-coding genes, and the mitochondrial 12S and 16S and the nuclear 18S and 28S ribosomal genes. Protein-coding genes were manually aligned, while the ribosomal genes were aligned with the MAFFT server (https:// mafft. cbrc. jp/ align ment/ server/, Katoh et al. 2009) using the Q-INS-i iterative refinement method. Conserved blocks were obtained from this alignment through GBlocks (Castresana 2000) using the less restrictive parameters from the GBlocks server (http:// molev ol. cmima. csic. es/ castr esana/ Gbloc ks_ server. html). These conserved blocks were used in the analyses. The dataset included all mitochondrial protein and ribosomal RNA genes and the nuclear 18S and 28S genes from all individuals, accounting for 19,584 nucleotides and 17 genes ( Table 2). Nautilus macromphalus was added to the matrix as outgroup.

Morphological analyses
To quantify beak shape, standardized digital photographs of the lateral side of each individual beak (both UB and LB) were obtained using two different camera sets depending on beak size to avoid distortion. Larger beak pictures were obtained with a digital camera (Olympus-Tough TG-5) mounted on an articulated arm; whereas pictures of smaller beaks were taken with a digital camera (Leica DFC450) attached to a binocular stereomicroscope (Leica MDG41). Scale was included in every picture and accounted for the further analytical process. To avoid any beak shape distortion due to dryness during the manipulation, they were humidified with 70% ethanol if any dryness or shape changes were detected. Then, a total of eight landmarks for the UB and ten for the LB were selected to describe beak shape based on the landmark configuration of Fang et al. (2017), which was modified from Neige and Dommergues (2002) (Fig. 2, Table 3). The Cartesian coordinates of each landmark were recorded on each picture using the software tpsDig2 v. 2.31 (Rohlf 2018). To verify that digitization was consistently performed and that digitization error did not influence our inferences, digitization was repeated twice by the same observer, for both the upper and lower beak. Then, we performed an ANOVA analysis with individual and repetition as factors, to investigate the contribution and relative variance explained by both factors on shape variation. Then, analyses of UB and LB landmark configurations were carried out by first applying a Generalized Procrustes Analysis (GPA), implemented using the package geomorph v. 4.0.0 (Adams et al. 2021;Baken et al. 2021) in the software RStudio v. 1.4.1106, working under R v. 4.0.5 (R Core Team 2021). This analysis eliminates effects not related to shape, including orientation, position, and size (Rohlf and Slice 1990), allowing the comparison of homologous structures, and provides shape variables for further analyses (see below).

Stable isotopes analyses
To obtain long-term information of the trophic niche of cephalopods, we analyzed the stable isotope values of C (δ 13 C) and N (δ 15 N) in their beaks Vigo et al. 2022). δ 15 N values are related to trophic level (Post 2002). δ 13 C values are related with primary production: in marine ecosystems the most productive areas are near the shore, therefore high levels of this marker indicate proximity to coastal habitats, while lower levels indicate greater association with more deep oceanic waters (Layman et al. 2011). Complete UBs were analyzed. Specifically, we analyzed the stable isotope values of ten individuals (five males and five females) per species, when available ( where isotopic analyses were performed. Samples were combusted at 1020 °C using a continuous flow isotope-ratio mass spectrometry system via a CON-FLO IV interface (Thermo Fisher Scientific, Bremen, Germany). All isotopic results were reported in the conventional delta (δ) per mil notation (‰), relative to Vienna Pee Dee Belemnite (δ 13 C) and atmospheric N 2 (δ 15 N). This process showed analytical measurement errors of ± 0.1 ‰ and ± 0.2 ‰ for δ 13 C and δ 15 N, respectively. Laboratory standards were previously calibrated with international standards supplied by the International Atomic Energy Agency (IAEA, Vienna).

Phylogenetic analyses
The partition scheme of the matrix included three segments: mitochondrial proteins, mitochondrial ribosomal RNA and nuclear ribosomal RNA genes. The Maximum Likelihood analysis (ML) was performed in IQTREE server (Hoang et al. 2018;Nguyen et al. 2015). The statistical support for each node is indicated after 1000 generations of the Shimodaira and Hasegawa like interpretation of aLRT statistic (SH-aLRT; Anisimova et al. 2011) and 10,000 ultrafast bootstrap iterations. We implemented the ModelFinder tool (Kalyaanamoorthy et al. 2017) in the IQTREE portal to estimate the best fitting model of substitution for each partition following the Bayesian Information Criterion (BIC), and selected these for the downstream analyses. A coalescent phylogenetic inference analysis was performed in BEAST v. 2.6.4 (Bouckaert et al. 2019). The input file was created using BEAUti. Site and clock models were independently set for each partition, based on the initial results of ModelFinder and selected using the extended options of the BEAST Package Standard Substitution Models SSM v. 1.0.1 (Bouckaert and Xie 2017). Clock models were set to relaxed log-normal models (Drummond et al. 2006). The prior of the species tree model was set to Yule model, and the birth rate was estimated by the analysis. Additionally, three fossil calibrations were applied to the analysis: crown Cephalopoda (Kröger and Mapes 2007), crown Coleoidea (Kröger and Mapes 2007) and crown Decapodiformes (Fuchs et al. 2013) with a minimum age of 408, 240 and 75 Ma, respectively. A Markov Chain (Drummond et al. 2002) of 100 million generations was run sampling every 10,000 generations. Chain convergence was examined with Tracer v. 1.7.2 (Rambaut et al. 2018) and ensured that ESS values over 200 were obtained. Finally, the initial 25% tree configuration was discarded as burn-in and the majority consensus tree obtained using TreeAnnotator. According to the IQTREE manual (Hoang et al. 2018;Nguyen et al. 2015), the SH a-LRT and ultrafast bootstrap values are considered accepted when above 80 and 95%, respectively. Posterior probabilities were accepted when values were above 0.95.

Ecomorphological analyses
Firstly, individual-level geometric morphometric and stable isotopic data were used to obtain a general characterization of variation among species in morphology and trophic niche. To this end, ANOVA tests were conducted using the function lm.rrpp of RRPP R package v. 1.0.0 Adams 2018, 2019) to test for interspecific differences in UB and LB shape, and δ 13 C and δ 15 N values, using residual randomization procedures for statistical evaluation. This approach was selected as it is known to perform well with highly dimensional (and univariate) data, such as those obtained through geometric morphometric techniques (Collyer et al. 2015;Collyer and Adams 2018). When ANOVA tests indicated significant differences among species, post-hoc pairwise analyses were performed as implemented in the function pairwise of RRPP Adams 2018, 2019). Afterwards, species means for both geometric morphometric and isotopic data were calculated, which were used for all subsequent analyses.
To examine main patterns of evolutionary shape variation considering phylogenetic relationships among species, a principal component analysis (PCA) was performed on species average values of UB and LB sets of Procrustes coordinates separately. We examined the position of different species in a phylomorphospace plot (Sidlauskas 2008) representing the two first principal component (PC) values of beak shape and the species phylogenetic relationships from BEAST analysis. PCA analyses and plots were implemented through the function gm.prcomp of geomorph R package (Adams et al. 2021;Baken et al. 2021). Moreover, for each PC we produced deformation grids to visualize specific morphological changes with respect to overall mean shape, using the function plotRefToTarget of geomorph R package (Adams et al. 2021;Baken et al. 2021). The degree of phylogenetic signal of each set of Procrustes coordinates, and of each of the isotopic markers, was estimated with the Kmult statistic using the function physignal of geomorph R package. This function evaluates the degree of phylogenetic signal in a dataset relative to that expected under a Brownian motion model of evolution for either multivariate (i.e. shape: Adams 2014) or univariate (i.e. isotopic markers: Blomberg et al. 2003) data.
The multivariate association between shape (UB and LB Procrustes coordinates) and the ecological variables (stable isotopes) taking the phylogenetic relationships among species into account was evaluated through the function phylo.integration of geomorph R package. For this purpose, the function performs a Partial Least Squares (PLS) analysis after projecting the data onto the phylogenetic covariance matrix, and then assesses the significance of phylo-projected PLS vectors using residual randomization procedures. We performed three sets of analyses, to explore the relationship of each of the UB and LB with both stable isotopic values, only δ 13 C and only δ 15 N.

Phylogenic results
The ML analysis performed with IQTREE (Supplementary files 1-2) recovered the monophyly of Octopoda, Oegopsida, Sepiida and Sepiolida. The splitting between the outgroup N. macromphalus and the remaining cephalopods was not statistically supported in our analysis, neither was the relationship between Octopodiformes and Decapodiformes, cephalopods with eight arms or with eight arms and two tentacles, respectively. Within Octopoda, our analysis recovered a clade of Octopus vulgaris and O. salutii sister to Scaeurgus unicirrhus and Pteroctopus tetracirrhus. The genus Bathypolypus clustered with the previous clade with no support (71% SH-aLRT and 55% ultrafast bootstrap), forming a group sister to Eledone. Decapodiformes received strong statistical support, while the relationships among the three decapod orders represented in this work (Oegopsida, Sepiida and Sepiolida) was not supported. Within Oegopsida two sister clades were defined with high support (99.4/100%). The first clade consisted of Histioteuthis reversa and Histioteuthis bonnellii. The second clade consisted of Illex coindetii and Todaropsis eblanae sister to Abralia veranyi (99.5/100%). The order Sepiida was represented by Sepia orbignyana and Sepia elegans clade sister to S. officinalis. Within Sepiolida, Rossia macrosoma and Neorossia caroli formed a clade. This clade was sister to Heteroteuthis dispar, with moderate or no statistical support (89/86%), and sister to Sepietta oweniana.
The BEAST analysis ( Fig. 3; Supplementary file 3) resulted in a different topology regarding the relationships among Oegopsida, Sepiida and Sepiolida, and between Eledone cirrhosa and Bathypolypus sponsalis.
The UB phylomorphospace (Fig. 4A) captured the highest amount of morphological variation, where PC1 explained 77.53% of the total shape variation and PC2 11.75%. The examination of PC1 deformation grids (Fig. 4A) indicated that species were differentiated based on hood proportions and lateral wall size. Negative PC1 values described beaks with expanded hoods, sharper rostra and generally smaller lateral walls (e.g., S. officinalis and H. bonnellii). Positive values described beaks with reduced hoods and expanded lateral walls (e.g., E. cirrhosa). PC2 deformation grids (Fig. 4A) revealed that species varied across this axis based on wing size, hood width and lateral wall and crest proportions. Negative values corresponded to beaks with reduced wings, wider hoods, dorsoventrally compressed lateral walls and slightly reduced crests (e.g., T. eblanae). Positive values represented beaks with more developed wing and crest, while having reduced lateral walls and narrower hoods (e.g., S. elegans).
The LB phylomorphospace associated to the first two principal component axes (Fig. 4B) captured a lower amount of morphological variation, with PC1 explaining 65.74% of shape variation and PC2 12.55% of the total variability. Corresponding PC1 deformation grids (Fig. 4B) showed that LB shape variation was related to lateral wall size proportions, wing width and rostrum-jaw angle distance (RJD). PC1 negative values described beaks with laterally compressed lateral walls, especially in the crest, narrower wings and stretched RJD (e.g., A. veranyi). In contrast, positive values along PC1 described beaks with elongated and dorsoventrally compressed lateral walls, wider wings and shortened RJD, with a residual jaw angle (e.g., B. sponsalis). PC2 (Fig. 4B) differentiated species based on lateral wall and crest length, as well as wing length and width. PC2 negative values described beaks with shorter and narrower wings and antero-posteriorly elongated lateral walls and crests (e.g., S. oweniana), and positive values corresponded to beaks with more enlarged wings in both width and length, and narrower lateral walls and crests (e.g., S. officinalis).
The comparison of phylomorphospace plots for both beak parts yielded different patterns of morphospace occupancy for the examined species (Fig. 4). The distribution of species in the UB phylomorphospace (Fig. 4A) showed a clear differentiation between species of the four orders, which also corresponded to significant phylogenetic signal (K = 1.39, p < 0.001). Octopuses were the most differentiated group, occupying the morphospace area Values on nodes indicate node ages (Ma) and posterior probability from the Bayesian inference analysis, and SH-aLRT support (%) and ultrafast bootstrap support (%) from the Maximum Likelihood analysis, respectively. The asterisk indicates the maximum support value, NA indicate differences in topology between both analyses tree and the dash indicates that the node was not supported. Yellow dots indicate fossil calibration point positions. The four main orders of the study and Nautilida (outgroup) are indicated. Illustrations correspond to each order represented in the tree and are based on Illex sp., Euprymna scolopes, Sepia officinalis, Octopus vulgaris and Nautilus macromphalus, from top to bottom that corresponds to maximum PC1 values and intermediate PC2. Within them, B. sponsalis was the most differentiated species, occupying an intermediate position in the morphospace between Octopoda and the other taxa. The other orders exhibited negative PC1 values, spreading across PC2; the negative area of PC2 was occupied by squids and the positive one by cuttlefishes and bobtail squids. Squids were clearly differentiated from the other groups and showed a relatively scattered distribution. Cuttlefish and bobtail squids ranged from minimum PC1 values to slightly higher than zero, Sepiida were located at Each axis is labeled with the percentage of the total shape variability explained by each principal component. Deformation grids show the morphological changes represented by principal components 1 (PC1) and 2 (PC2). In the deformation grids, colored dots indicate the landmark mean coordinates between all species and arrows indicate deformation from the mean. Shape patterns were exaggerated 1.5 times for PC1 and 2 times for PC2. The name of each species corresponding to the abbreviation are indicated as follows: Aver = Abralia veranyi; Bspo = Bathypolypus sponsalis; Ecirr = Eledone cirrhosa; Hdis = Heteroteuthis dispar; Hbon = Histioteuthis bonnellii; Hrev = Histioteuthis reversa; Icoi = Illex coindetii; Ncar = Neorossia caroli; Osal = Octopus salutii; Ovul = Octopus vulgaris; Ptet = Pteroctopus tetracirrhus; Rmac = Rossia macrosoma; Suni = Scaergus unicirrhus; Sele = Sepia elegans; Soff = Sepia officinalis; Sorb = Sepia orbignyana; Sowe = Sepietta oweniana; Tebla = Todaropsis eblanae lower PC1 values than Sepiolida. Within Sepiolida, Heteroteuthis dispar was located between Octopoda and the other clades, being the most differentiated of its group. Across the LB phylomorphospace (Fig. 4B) differentiation between orders was not as clear as in the UB, also corresponding to lower but significant phylogenetic signal (K = 0.76; p < 0.001). In this morphospace, differences among orders were still readily visible, but the key difference between this pattern of distribution and the previous is the differentiation of species within orders. Here, differences in distribution between close relatives appear to be more prominent and differences between orders seem to be less relevant, resulting in a less structured distribution of the species. Octopuses were distributed in the maximum PC1 values, scattered along the PC1 axis. Squid species were distributed in the lowest PC1 values and near zero PC2 values. This clade was the least scattered, except for H. reversa which was isolated at near maximum PC2 values. Cuttlefishes and Sepiolida had a wide distribution from intermediate negative PC1 values to positive near zero and from minimum to maximum PC2 values, represented by S. oweniana and S. officinalis, respectively. Sepiida occupied a position with higher PC1 and PC2 values than Sepiolida.
In stable isotope space (Fig. 5), species were distributed mainly in one central group, except for O. vulgaris and S. officinalis, which both had higher δ 13 C and δ 15 N values. Within the central group, we observed more variation along the δ 13 C values and high overlap among species. Regarding the δ 15 N values, less overlap was observed and two groups were distinguished. Octopoda species ranged between low and high δ 13 C values and were located at lower δ 15 N values (except for O. vulgaris). Squids shared low to medium δ 13 C and δ 15 N values. Sepiida were highly scattered, from high δ 13 C and δ 15 N values to both maximums, where S. officinalis was located. Bobtail squids were spaced from low to medium δ 13 C values and low-medium δ 15 N values. Isotopic values did not exhibit significant phylogenetic signal (p > 0.05 in both cases). Integration tests between both beaks and the stable isotopic values did not reveal any significant relationship (p > 0.05; Table 4).

Discussion
Here, we obtained a robust phylogeny and described the morphological variation in beak shape among the selected species to test the influence of phylogenetic relatedness and trophic ecology over beak shape diversity. On one side, a high phylogenetic signal has been detected for both beaks indicating a close relationship between shape and phylogenetic relatedness, while on the other no relationship has been found between trophic parameters and shape. These results suggest a strong influence of phylogenetic relatedness over beak shape.
The relationships obtained through the ML analysis were well-supported. In previous studies, the phylogeny among the four studied orders have been explored (Anderson and Lindgren 2021;Lindgren 2010;Lindgren et al. 2012) without reaching a conclusive result regarding the relationship among decapodiform lineages. In this work, a tritomy appeared in the relationship among Oegopsida, Sepiida and Sepiolida. In the same way, the relationships among families within the order Octopoda were not clear as we obtained a polytomy among Eledonidae, Bathypolypodidae and Octopodidae (including O. salutii, O. vulgaris, P. tetracirrhus and S. unicirrhus). Although the number of analyzed individuals was small, the relationships found between Abralia, Todaropsis and Illex was concordant with the results of Fernández-Álvarez et al. (2022), while the position of Histioteuthidae was unresolved in that work. Within Sepiida, our phylogenetic tree indicated a clear differentiation between the clade formed by S. elegans and S. orbignyana respective to S. officinalis with strong statistical support, similar to previous data (Khromov 1987;Sanjuan et al. 1996;Yoshida et al. 2006). The complete mitochondrial phylogeny from the recent study of Sanchez et al. (2021) supported our topology with Rossinae and Heteroteuthinae sister to Sepiolinae, although the analyses based on nuclear data did not.
Geometric morphometrics revealed clear morphological differences between the 18 cephalopod species included in the present study for both UB and LB. Interestingly, these interspecific differences are also accompanied by significant phylogenetic signal, which suggests that beak shape is tightly associated with phylogenetic relationships, at least for the UB. This result likely explains why beak morphology has been so successfully exploited as an identification tool (Clarke 1986;Xavier and Cherel 2009). The UB shape is clearly more influenced by phylogeny than the LB, as reflected by the more structured distribution of the species in morphospace as well as the higher phylogenetic signal. UB shape reflects more of those differences related to distant groups than closer ones, thus resulting in a better taxonomic order differentiation and consequently in a higher phylogenetic signal. The LB still exhibits a significant phylogenetic signal, but the distribution of species in the morphospace is not as ordered with respect to the phylogeny. Besides, shape differences among closely related species are more marked when considering the lower than the upper part of the beak. This capability to distinguish between closely related species is a key aspect that makes the LB much more suitable for species identification than the UB. The literature also supports this idea, as in cephalopod identification beak guides the LB is preferentially used and empirically achieves better results (Clarke 1986;Tan et al. 2021;Xavier and Cherel 2009).
Stable isotope values revealed interspecific ecological differences, reflecting variations in the trophic niche and main habitats of the studied species (e.g. Hernández-García 1992; Jereb et al. 2014;Regueira et al. 2017). Based on the stable isotope values of N and C together, we found a clear differentiation into two groups: a larger group formed by species that inhabit deeper habitats and a smaller one formed by species present on the continental shelf (O. vulgaris and S. officinalis). The larger group appears to be divided along the δ 15 N values. The group with lower δ 15 N values is formed by A. veranyi, I. coindetii, S. orbignyana, all sepiolids and the remaining Octopoda species; while the group with higher values is formed The first group is represented by offshore demersal species belonging to every order that preys mainly on crustaceans-Octopoda (Jereb et al. 2014;Quetglas et al. 2001Quetglas et al. , 2005Quetglas et al. , 2009Regueira et al. 2017), Oegopsida (Hernández-García 1992), Sepiida and Sepiolida (Guerra-Marrero et al. 2020;Roper 2005, 2010;Vafidis et al. 2009)-, thus showing lower δ 15 N levels; while the second one is formed mainly by pelagic species that prey mainly on fishes, followed by crustaceans, and therefore belonging to higher trophic levels (Castro and Guerra 1990;Hernández-García 1992;Roper 2005, 2010;Quetglas et al. 2010). Although the stable isotope values and their interpretation are concordant with the published diet of these species, it is important to note that these trophic markers do not provide an accurate estimation of their diet composition or the physical quality of the prey tissues. For this reason, new studies including the type and morphology of prey present in the diet of these cephalopods could help to a better understanding of the functional relationship between diet and beak shape. Integration tests performed to explore the potential association between isotopic values and UB or LB shape did not identify any significant covariation pattern. Thus, this study could not establish a relationship between beak shape and dietary habits, in this case estimated by isotopic signatures. Hence, our results suggest that beak shape is not tightly related to the animal's capability to capture and process prey or to its feeding performance. A possible explanation for this might be many-to-one mapping, which is the capability of different phenotypes to achieve a similar functional property, a phenomenon previously described for many other systems (Blob et al. 2006;Guderley et al. 2006;Lappin and Husak 2005). In such a case, the remarkable variation in beak shape observed across the species investigated here, would not be reflected in their dietary habits. Similarly, functional decoupling might weaken the link between beak morphology and functionality: in this process, an important, already existing functionality can be performed by a different structure of the organism. Indeed, the beak is not the only structure related to prey capture and feeding. Structures such as the arms or tentacles, equipped with suckers, sometimes modified into capture hooks, and the radula help in hunting, seizing and fractioning prey (Villanueva et al., 2017). The involvement of cephalopod limbs in prey capture, the external digestion and the radular abrasion in food fractionation may result in lower selective pressure on the beak, allowing for non-adaptive morphological variation, at least with respect to feeding ecology. The morphological variation we observe here in beak shape may be related to these factors, resulting in the observed pattern of phenotypic differentiation across species that corresponds to shared evolutionary history rather than dietary habits. This would also explain the high phylogenetic signal of beak shape: as closely related species are more similar to each other than distant ones, beak shape change would correspond to the phylogenetic differentiation of the group. Using novel 3D morphometric analyses, Roscian et al. (2022) found that the shape of the cephalopods can be also associated with particular ecological factors (e.g., their habitat) or mastication behaviours (e.g., fast closing or hard biting). The results of both methodologies suggest that different morphometric approaches might shed light on different evolutionary trends. Thus, future studies combining both 2D and 3D morphometric analyses under well-resolved phylogenies might help to understand how both phylogenetic, functional and ecological constraints drive cephalopod beak shape evolution.

Conclusions
Genome skimming provides robust phylogenetic results in analyses based on large regions of the genome, such as the complete mitogenomes and the nuclear ribosomal genes used in this study. This robust phylogeny showed that cephalopod beak morphological evolution is strongly driven by the phylogeny, while no effect of ecological pressures on beak morphology was detected. The strong phylogenetic signal for both UB and LB supports the use of beaks as an identification tool since beak shape strongly reflects phylogenetic relationships. Our results also explain why LBs have traditionally been more successful for identification purposes, as this structure allows for a more precise identification between closely related species. Our approach proves that well-supported phylogenies can be used for assessing the most important drivers in cephalopod morphological evolution. It is expected that 3D morphometric methods and more accurate trophic estimations could provide better results in studies among additional related taxa.