Sexual Selection and the Evolution of Male Reproductive Traits in Benthic Octopuses

Competition between same-sex organisms, or intra-sexual selection, can occur before and after mating, and include processes such as sperm competition and cryptic female choice. One of the consequences of intra-sexual selection is that male reproductive traits tend to evolve and diverge at high rates. In benthic octopuses, females often mate with more than one male in a single reproductive event, opening the venue for intra-sexual selection at multiple levels. For instance, males transfer spermatophores through hectocotylus, and can remove the spermatophores left by other males. Considering the limited evidence on post-copula competition in benthic octopuses, and the potential to affect the evolution of reproductive traits within octopodids, we put this hypothesis to a test employing a phylogenetic comparative approach. We combined data on hectocotylized arm length (HAL), ligula length (LL), spermatophore length (SL) with a Bayesian molecular phylogeny of 87 species, to analyze how reproductive traits have diverged across lineages and covary with body size (mantle length; ML). First, additionally to ML, we estimated the phylogenetic signal (λ) and mode of evolution (κ) in each reproductive trait. Second, we performed phylogenetic regressions to quantify the association among reproductive traits and their co-variation with ML. This analysis allowed us to estimate the phenotypic change along a branch into the phylogeny, and whether selection may have played a role in the evolution and diversification of specific clades. Estimations of λ were always high (>0.75), indicating concordance between the traits and the topology of the phylogenetic tree. Low values of κ (<1.0) suggested that evolution depends on branch lengths. All reproductive traits exhibiting a positive relation with ML (β > 0.5 in all cases). Overall, evolutionary rate models applied to the SL-ML regression suggested that octopuses of the family Megaleledonidae have evolved larger spermatophores than expected for their size. The regression HAL-ML indicated that HAL was more variable in Megaleledonidae than in the remaining clades, suggesting that the high divergence across species within this group might partially reflect intra-sexual selection. These results support the hypothesis that, at least in some lineages, sexual selection may account for the divergence in reproductive traits of male octopuses.

Competition between same-sex organisms, or intra-sexual selection, can occur before and after mating, and include processes such as sperm competition and cryptic female choice. One of the consequences of intra-sexual selection is that male reproductive traits tend to evolve and diverge at high rates. In benthic octopuses, females often mate with more than one male in a single reproductive event, opening the venue for intrasexual selection at multiple levels. For instance, males transfer spermatophores through hectocotylus, and can remove the spermatophores left by other males. Considering the limited evidence on post-copula competition in benthic octopuses, and the potential to affect the evolution of reproductive traits within octopodids, we put this hypothesis to a test employing a phylogenetic comparative approach. We combined data on hectocotylized arm length (HAL), ligula length (LL), spermatophore length (SL) with a Bayesian molecular phylogeny of 87 species, to analyze how reproductive traits have diverged across lineages and covary with body size (mantle length; ML). First, additionally to ML, we estimated the phylogenetic signal (λ) and mode of evolution (κ) in each reproductive trait. Second, we performed phylogenetic regressions to quantify the association among reproductive traits and their co-variation with ML. This analysis allowed us to estimate the phenotypic change along a branch into the phylogeny, and whether selection may have played a role in the evolution and diversification of specific clades. Estimations of λ were always high (>0.75), indicating concordance between the traits and the topology of the phylogenetic tree. Low values of κ (<1.0) suggested that evolution depends on branch lengths. All reproductive traits exhibiting a positive relation with ML (β > 0.5 in all cases). Overall, evolutionary rate models applied to the SL-ML regression suggested that octopuses of the family Megaleledonidae have evolved larger spermatophores than expected for their size. The regression HAL-ML indicated that HAL was more variable in Megaleledonidae than in the remaining clades, suggesting that the high divergence across species within this group might partially reflect intra-sexual selection. These results support the hypothesis that, at least in some lineages, sexual selection may account for the divergence in reproductive traits of male octopuses.

INTRODUCTION
The evolution of mate choice and mating competition has been a major component of Darwin's theory of sexual selection (Darwin, 1871). Since then, it has been generally accepted that different selective pressures acting on male and female attributes may give rise to sexually dimorphic traits, which are often interpreted as evidence of direct competition for mates within a given sex, differential success to attract potential mates from the opposite sex, or gametic competition (see Basolo, 1990;Birkhead, 1992;Evans and Sherman, 2013). At the gametic level, competition can occur before and after mating through sperm competition or cryptic female choice. Whereas sperm competition involves strategies by the male to either remove, displace or inhibit the sperm of other males, cryptic female choice constitutes femalebiased selection through the use/removal of sperm to fertilize their eggs. Competition at the gametic level has been described in organisms from several phyla, such as insects, molluscs, birds and mammals (Mann, 1984), and is generally enhanced in polyandric species where females can mate with multiple males in a single reproductive episode (Birkhead, 1998;Gomendio, 2002). Because of its relevance to males' reproductive success (Eberhard, 1998;Snook, 2005), multiple responses have evolved to outcompete rival males, including: (i) the production of spermatophores or packages of sperm (Mann, 1984;Nigmatullin et al., 2003) that, when transferred to females, may occupy considerable space within the storage organs preventing other spermatophores from being stored (Thornhill and Alcock, 1983), (ii) the removal of other males' spermatophores during copulation (Cigliano, 1995), or (iii) the production of sperm with inhibitory effects on the rival males' sperm function (Snook, 2005). Because of this evolutionary arms race, the genitalia of several organisms exhibit extreme differences in size and shape across closely related species and are presumed to evolve faster than other traits (Eberhard, 1985;Genevcius et al., 2017).
Polyandry, sexual dimorphism and sexual selection have been described in several lineages of cephalopod molluscs and is widespread in this group (Mann, 1984;Hanlon and Messenger, 1996;Squires et al., 2012). Polyandrous behavior has been observed in female octopuses, potentially increasing post-mating sexual selection and driving the evolution of a myriad of sperm transfer strategies (e.g., male mounting female; Figure 1A; for other examples see Hanlon and Messenger, 1996;Cheng and Caldwell, 2000;Huffard et al., 2008;Gutierrez et al., 2012;Morse and Huffard, 2019). In all cases, male octopuses pack their sperm into spermatophores and transfer them to females by using a modified arm called hectocotylus. This specialized arm employed to deposit the sperm packages into the females' pallial cavity is characterized by two well defined segments, the calamus and the ligula (see Hanlon and Messenger, 1996;Wodinsky, 2008;Marian, 2015), and by a considerable inter-specific morphological variation ( Figure 1B). This variation has been associated with the successful transference of spermatophores during mating (Robson, 1926); however, direct behavioral evidence on their role in removing or breaking down spermatophores from rival males remain speculative (Cigliano, 1995;Hanlon and Messenger, 1996;Norman et al., 2004), providing an important framework for evaluating untested hypotheses on sexual selection (see Voight, 2009). Furthermore, spermatophores also exhibit considerable inter-specific differences in size (e.g., ranging from 7 to 1130 mm in length), even after accounting for size effects (Voight, 2009). This is likely because individuals with larger spermatophores have greater sperm reservoirs and consequently much more sperm to fertilize females' eggs (Voight, 2001). Nonetheless, in spite of the high levels of morphological variation in these traits, the evolution of hectocotyli and spermatophores across benthic octopuses remains poorly understood, as previous comparative studies have not accounted for the evolutionary history of the lineages involved (see Voight, 2001Voight, , 2002Voight, , 2009).
Here, we study the evolution and diversification of these reproductive traits across benthic octopuses, employing phylogenetic analytical methods that take into consideration patterns of relatedness between different lineages. Phylogenetic methods are currently indispensable to understand patterns of phenotypic diversification and their underlying processes, as well as the direction and magnitude of inferred evolutionary changes (Felsenstein, 1985;Harvey and Pagel, 1991;Rezende and Diniz-Filho, 2012). Accordingly, recent phylogenetic comparative studies have been quite successful in reconstructing the evolution of life-history strategies in cephalopods in response to different environmental pressures (see Lindgren et al., 2012;Ibáñez et al., 2014Ibáñez et al., , 2018Pardo-Gandarillas et al., 2018). In the present work, we used a phylogenetic approach to: (a) reconstruct how hectocotyli and spermatophores have evolved along a molecular phylogeny including 87 species of benthic octopuses, (b) explore the correlated evolution between these traits and body size, and (c) employ variable-rates phylogenetic regression to determine which clades exhibit abnormally high rates of phenotypic evolution in response to selection. Even though results could be interpreted as putative evidence of strong post-copulatory sexual selection, we also discuss alternative adaptive scenarios that might have given rise to the observed differences across lineages.

Dataset and Phylogeny
We performed an extensive literature review to obtain information on the variability of reproductive traits across different lineages of benthic octopuses. For our analyses, we selected descriptive measures that have been extensively studied with relatively standardized protocols and that could, therefore, be readily compared across different studies (see Table 1): mantle length (ML), arm length (AL), ligula length (LL), hectocotylized arm length (HAL), and spermatophores length (SL). Subsequently, we combined this information with a new phylogenetic hypothesis of benthic octopuses encompassing a total of 97 species, including outgroups (Mendeley Datasets: doi: 10.17632/5vkm46hm49.1), that are based on three mitochondrial genes (16S ribosomal RNA, Cytochrome oxidase I, Cytochrome oxidase III) and one nuclear gene (Rhodopsin). A detailed explanation regarding the analyses underlying the phylogenetic reconstruction, estimation of uncertainty and validation of our working phylogeny has been provided elsewhere . Briefly, phylogenetic relationships were inferred from a partitioned matrix (16S, COI + COIII, RHO) with a different substitution model for each gene. This matrix was composed of 97 species, including 88 species from the superfamily Octopodoidea, and two species from the superfamily Argonautoidea, six cirrates and the vampire squid Vampyroteuthis infernalis as outgroups. Bayesian analyses were conducted using MrBayes 3.2 with four chains, each with ten million generations, sampled every 1,000 generations. The first 1,000 trees of each run were discarded as burn-in, and a consensus of the remaining trees was calculated. For simplicity, we have excluded the outgroups, providing the reconstructed relationships for the 87 benthic species in the dataset (excluding Vitreledonella richardi, Figure 2A).

Statistical Analyses
We performed three complementary analyses to determine how the variation in reproductive traits is related to phylogenetic history. First, we performed univariate analyses to estimate the amount of phylogenetic signal (λ) and the mode of evolution (κ) of each trait (ML, HAL, LL, and SL) employing BayesTraits v3 (Meade and Pagel, 2017). Second, we ran three separate variable-rates regression models (Baker et al., 2016) to determine how HAL, LL, and SL vary as a function of body size (ML), and also diagnose in which lineages these structures have diverged more than expected after controlling for size effects. Third, we performed a phylogenetic principal component analysis (PCA) to determine the degree of covariation between reproductive traits and study correlated evolution between them after removing body size effects.
Phylogenetic signal in our univariate analyses was estimated employing Pagel's λ, which quantifies the tendency of closely related lineages to resemble each other in comparison to a Brownian motion model of evolution (Pagel, 1999(Pagel, , 2002: λ = 1 indicates that the distribution of the phenotypic traits along the tips of the phylogeny closely resemble the expectation based on Brownian motion (i.e., high phylogenetic signal), whereas λ = 0 shows that patterns of phenotypic resemblance due to shared phylogenetic history is negligible (i.e., low phylogenetic signal). The mode of phenotypic evolution was estimated using Pagel's κ, which scales the branch lengths between their original values and a single constant, mimicking gradual evolution when κ = 1, and a punctuated model of evolution when κ = 0 (Pagel, 1999(Pagel, , 2002. The posterior distribution of all parameters was visualized in the software Tracer V1.6 (Rambaut et al., 2013). To test if λ and κ estimated values were different from pre-established values, we first estimated the higher posterior distribution of these parameters for each trait in BayesTraits V3, forcing each parameter to have a value of 0 for λ (i.e., no signal) and 1 | Summary of the studied species, with information on their body size (ML, maximum mantle length) and reproductive traits (SL, maximum spermatophore length; LL, maximum ligula length; HAL, maximum hectocotylized arm length; AL, maximum arm length; EL, maximum egg length).

Species
ML (   0 and 1 for κ (i.e., perfectly punctuated or perfectly gradual evolution, respectively) and compared the fit of estimated vs. forced models with log 10 Bayes Factor (BF). The larger the BF value, the better the fit of the estimated model in comparison against the forced one, with BF > 0.5 being generally interpreted as strong support for the estimated model and BF > 1 being considered decisive (Kass and Raftery, 1995). Because these univariate analyses include scaling effects, we expect HAL, LL, and SL to exhibit less signal (i.e., a lower λ) than ML if they evolve faster than this trait due to an evolutionary arms race (i.e., signal is expected to decrease if traits diverge fast in response to selection).
To estimate phenotypic selection on reproductive traits after removing potential scaling effects, we performed three separate variable-rates regressions with log 10 -transformed HAL, LL, and SL as a function of ML. This regression model was recently developed by Baker et al. (2016) and allows the rate of change to vary through the phylogenetic branches and identifies areas of the tree where the rate of evolution departs significantly from background levels (Venditti et al., 2011). With this purpose, this regression method estimates a branch-specific metric V / B that contrasts the expected phenotypic variance V along the branch due to changes in evolutionary rates (i.e., acceleration or deceleration) vs. the expectation attributable to the background evolutionary rate ( B ). Branches in which the amount of estimated phenotypic change doubles the background rate ( V / B > 2) constitute regions of the phylogeny that were likely under positive selection (Baker et al., 2016). We implemented the variable-rates regression model with the phylogenetic independent contrast regression module in BayesTraits, employing the reversiblejump Markov chain Monte Carlo (RJMCMC) to determine whether V / B > 2 results were observed in more than 95% of the posterior distribution. In all Bayesian analyses described above, we ran 20,000,000 iterations via the MCMC method. Parameters were sampled every 1,000 iterations, excluding the first 25% of iterations. The 95% highest posterior density (95% HPD) for each parameter was calculated in Tracer, and all analyses performed in BayesTraits; outgroup and sister groups were excluded.
To determine to what extent the different reproductive traits studied here evolve in tandem, we performed a phylogenetic PCA (Revell, 2009) including log 10 -transformed ML, HAL, LL, and SL. To account for phylogenetic signal, λ was estimated concomitantly with parameters from the PCA. Because we were primarily interested in identifying potentially contrasting evolutionary strategies between lineages, we focused on the second and third principal components (PC2 and PC3) that provide information on differences in morphology or shape after removing size effects embedded in the first principal component.
Finally, to explore the association response to selection on other traits, the correlation between reproductive traits [spermatophore length (SL) and egg length (EL)] and morphological traits [arm length (AL) and HAL] were analyzed using phylogenetic generalized least squares (PGLS) regressions (Pagel, 1999). To account for phylogenetic signal, λ was estimated concomitantly with parameters from the regression model (Pagel, 2002). After reviewing different sources, egg length data was obtained only for 72 species (Table 1).
All traits exhibited a high phylogenetic signal, with λ > 0.75 statistically different from λ = 0, as indicated by BF > 11 in all traits ( Table 2). As hypothesized (see section "Materials and Methods"), λ estimated for reproductive traits were generally lower than values for ML, agreeing with the observation that these traits exhibited more variation across species than ML. Regarding the mode of evolution inferred by κ, calculated in combination with λ in the univariate analyses, estimates for ML, HAL, LL, and SL were intermediate between κ = 0 and 1 and statistically different from those values according to BF estimates (BF > 1.39 for all comparisons) ( Table 2). This implies that all traits evaluated tended to evolve slower than predicted (i.e., evolutionary stasis) in longer branches when compared to shorter branches (Pagel, 1999(Pagel, , 2002. Phylogenetic regressions of reproductive traits as a function of ML indicated that all variables scaled positively with body size, with scaling exponents corresponding to 1.18 for HAL (95% HDP between 0.99 -1.36), 0.73 (0.53 -0.95) for LL, and 0.90 (0.71 -1.10) for SL (Figure 3). Consequently, HAL tends to become disproportionally larger as size increases, whereas SL scales roughly isometrically, and LL is relatively shorter in larger lineages. According to variable-rates phylogenetic regressions controlling for these scaling effects, several regions of the phylogeny exhibited accelerated rates of phenotypic evolution, and met the criterion of V / B > 2 proposed as evidence of positive selection (Figure 2). This was particularly true for HAL and SL, for which we detected selection in a total of 33 and 44 branches, respectively, or roughly 20 to 30% of all branches (Table 3). Interestingly, separate analyses for both traits gave rise to similar qualitative results, suggesting accelerated rates of phenotypic divergence for these traits in Antarctic, and deep-sea octopuses from the family Megaleledonidae (Figures 2B,D). In contrast, evidence of positive selection in LL was limited to only 8 branches, or 4.9% of the total ( Table 3), most of them involving the Cistopus clade ( Figure 2C).
The phylogenetic PCA including log 10 -transformed ML, HAL, LL, and SL strongly supported correlated evolution between these reproductive traits (Figure 4). As expected, PC1 accounted for a substantial fraction of the variance in the original data (70.8%), which could be attributed to a variation associated with body size, whereas the remaining PCs involve phenotypic variation that is independent of size (i.e., "shape" for simplicity; Figure 4). Accordingly, ML loadings in the remaining PCs were very low because most variation in this trait was explained by PC1. After removing the effects of size, PC2, and PC3 combined accounted Numbers in parentheses correspond to the highest posterior density (95% HDP) and represent the 95% credibility confidence estimated with a Bayesian approach, whereas the Bayes factor (BF) provides the weight of the evidence.
for 91.1% of the variance in shape observed across lineages, with loadings indicating that most of the variance in HAL is explained by PC2 and in SL by PC3, with LL falling somewhere in between (Figure 4). Contrasting these results against the outcome of the variable-rates regressions, we can identify two distinct groups (Figure 4), one exhibiting reduced hectocotyli and large spermatophores (low HAL and high SL), and the other with relatively large hectocotyli with small ligulae (high HAL and low LL). As clearly illustrated in Figure 4, results from variablerates regressions performed separately for each reproductive trait provided very consistent results and complementary evidence of positive selection across the same phylogenetic lineages. Spermatophore length did not correlate with egg length (n = 72 species, r = 0.050, 95% HPD between 0.0013 and 0.1082, λ = 0.85, Supplementary Figure S1). Interestingly, the correlation between arm length and hectocotilized arm length was higher than zero (n = 85 species, r = 0.8124, 95% HPD between 0.7946 and 0.8372, λ = 0.75, Supplementary Figure S2).

DISCUSSION
The present results support our original hypothesis on reproductive traits in male benthic octopuses, evidencing that spermatophores, and hectocotyli (arm and ligula) exhibited accelerated rates of evolution, at least in several Antarctic and deep-water lineages (Megaleledonidae and Cistopus), presumably due to sexual selection. All reproductive traits showed a foldrange of morphological variation that is substantially larger than expected based solely on differences in body size (i.e., ML), and variable-rates regression analyses clearly indicated that several lineages tended to deviate substantially from allometric expectations. While our analyses focussed primarily on the variation in reproductive traits after statistically removing the effects of size (effectively using levels of ML divergence between lineages as a standard of comparison), body size may have also been under selection in some lineages within Octopodoidea, as observed in other clades exhibiting sexual dimorphism in size (e.g., Amphitretidae and Tremoctopodidae Jereb et al., 2014). Unfortunately, this possibility cannot be directly tested with our phylogenetic comparative approach in the absence of precise information on most aspects of reproductive behavior of the analyzed species, including competition for mates, mate choice and mating position, as well as their intraspecific body size variation. Additionally, other environmental variables and selective pressures may have contributed to body size evolution of many of these lineages. Nonetheless, it is important to consider that body size may also evolve in response to sexual selection and gametic competition. For instance, we detected a clear positive association between ML and SL, indicating that larger speciesand presumably larger individuals within a species -have bigger spermatophores and consequently more sperm to transfer to females. Accordingly, this same trend was previously described by Voight (2009); therefore, we do not only provide support for such finding within a strict phylogenetic context, but we also detected which groups and lineages deviated from allometric expectations (see below).
Admittedly, while our analyses provided strong evidence of selection in several regions of the phylogeny, some limitations must be highlighted. First, note that our phylogenetic analysis detects regions of the phylogeny with extraordinary rates of evolution in comparison to background rates inferred from the same dataset, which is inherently conservative and can only detect selection in restricted regions of the phylogeny. Therefore, it is possible that we might have missed other evolutionary clades whose phenotypic diversification might be partly explained by selection (while decreasing the V / B > 2 might partly circumvent this problem, it would also increase the type I error). Second, our analyses do not inform specifically on the mechanisms that underlie these results. Consequently, alternative adaptive scenarios must be taken into account to determine the likelihood that observed patterns emerge from sexual selection. In this context, we believe that two possibilities are worth considering: (1) results reflected adaptation to environmental conditions, and/or (2) they partly reflected selection on correlated traits. With regard to the first scenario, the high evolutionary rates (i.e., V / B > 2) in HAL and SL involved primarily lineages of Antarctic and deep-sea octopuses from the family Megaleledonidae. Because ectotherm organisms inhabiting cold waters tend to exhibit lower fecundity, slower growth rates, and larger life-spans (van Voorhies, 1996;Ibáñez et al., 2018), it is plausible to expect that cold-adapted lineages may evolve larger spermatophores compared to warm-water species (Voight, 2009). Moreover, among deep-sea organisms that live in low densities the probability of a mating encounter is reduced (see Hoving et al., 2012), and therefore, selection may favor a high reproductive investment per mating. Consequently, we suggest it is possible that the colonization of cold-waters and deep-sea habitats might partly explain the high evolutionary rates detected in this clade (Megaleledonidae) and their larger spermatophores. Alternatively, it is also possible that some of the patterns detected reflect correlated responses to selection on other traits, which might be particularly true for HAL given its close association with AL (r = 0.81). All lineages detected in the variable-rates regression for HAL exhibited smaller arms than predicted from allometry, though it is not clear exactly which selective pressures might favor smaller arms. Importantly, while these alternative evolutionary scenarios might justify why members of the family Megaleledonidae differed from other groups, they failed to explain the extremely high diversity within this clade and its degree of phenotypic variation (i.e., HAL). Within Antarctic octopods, the family Megaleledonidae is the most diverse, with new species still being discovered (Xavier et al., 2018), and here we show that this highly speciose group also exhibited extremely elevated levels of phenotypic divergence in male reproductive traits (i.e., SL). We contend that the speciation rates observed in this clade in conjunction with the extremely high rates of phenotypic evolution cannot be explained by niche diversification, and likely reflect sexual selection (i.e., "runaway" sexual selection), where the coevolution of female mating preferences and male sexual characters promotes reproductive isolation and foments speciation (see Lande, 1982). The process of sperm competition has been well described in benthic octopuses, highlighting the role of cephalopod behavior in mediating intra-sexual competition (e.g., several males attempting to mate with a female simultaneously; reviewed by Hanlon and Messenger, 1996). Similarly, the occurrence of multipaternity has also been described in some species of octopuses, suggesting that females are able to fertilize eggs with the sperm of multiple males, decreasing the probability of fertilizing high number of eggs with the sperm of single male, as reported for Graneledone boreopacifica (Voight and Feldheim, 2009), Octopus vulgaris (Quinteiro et al., 2011), E. dofleini (Larson et al., 2015), O. minor (Bo et al., 2016), Hapalochlaena maculosa (Morse et al., 2018), and O. oliveri (Ylitalo et al., 2019). Additionally, it has been proposed (but not verified by other authors, nor by our own data) that species with large-sized ligula do not only use this structure for spermatophores transfer, but also to breakdown or modify the position of spermatophores from rival males (Cigliano, 1995;Hanlon and Messenger, 1996;Norman et al., 2004). Nonetheless, this behavior has been described in other taxa (such as insects) that use their copulatory organs to extract the sperm left by other males as a mechanism to counteract sperm competition (Birkhead and Moller, 1999). Finally, cryptic female choice favoring larger spermatophores has been reported in the sepiolid squid Idiosepious paradoxus through postcopulatory behavior (see Sato et al., 2013), and may therefore be taking place in closely related octopods. While these studies leave no doubt that sexual selection is potentially an important factor shaping the evolution of benthic octopuses, the fact that this seems to be particularly the case for members of the family Megaleledonidae remains an open question. It is possible that this group has evolved mating strategies and reproductive habits that exacerbate sperm competition via, for instance, territoriality or female We included log 10 -transformed ML, HAL, LL and SL, removed PC1 that encompassed primarily scaling effects, and worked with the remaining components that indicate differences in "shape" across lineages. (B) Species under positive selection according to variable-rates regression are shown in colors, as in Figures 2, 3. Note that selection was detected in more than a single trait in several lineages. postcopulatory selection. Interestingly, Rocha et al. (2001) speculated that the wide range of sizes of maturing ova described in the megaleledonids Pareledone charcoti and Adelieledone polymorpha could indicate repeated spawning in these species, contrasting with the majority of octopods that are considered terminal spawners. Another potential explanation that is not mutually exclusive corresponds to that, due to the environmental conditions encountered by these Antarctic deep-sea species (i.e., temperature and environmental stability; see Ibáñez et al., 2018), the impact of sexual selection become disproportionally important in this group in comparison to other benthic octopuses. Indeed, the lack of correlation between SL and EL suggest the absence of environmental selection on SL at the poles or at deep water environments.
In other words, we speculate that other factors shaping phenotypic evolution, such as predation, interspecific competition or environmental heterogeneity, may be relatively less important in the Antarctic deep-sea species. Perhaps the combined action of these two phenomena, namely the evolution of exclusive reproductive strategies in this clade in response to specific environmental pressures, may ultimately explain the very strong signal of selection and phenotypic divergence detected across males of this family. Overall, our phylogenetic approach provides some evidence of sexual selection within benthic octopuses, particularly for Megaleledonidae, and a potentially relevant role in their diversification. Detailed studies on different mating behaviors and how they relate with morphological and life-history traits are still necessary to better understand the adaptation of different cephalopod lineages to highly contrasting environments worldwide.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: https://data.mendeley.com/datasets/ 5vkm46hm49/1.

AUTHOR CONTRIBUTIONS
CI, ER, and MP-G conceived the idea, designed the study, analyzed the data, and led the writing of the manuscript. SC, JP-Á, and JC collaborated in literature review, writing, and provided the editorial advice. All authors have read and commented on the manuscript.

FUNDING
This work was partially funded by the FONDECYT research grants 1181153, 11170617, and 1170017 awarded to CI, SC, and ER, respectively. The additional support from the INACH research grant RG 50-18 awarded to MP-G and the CONICYT PIA/BASAL FB0002 to ER are also appreciated.