Thorson's rule, life-history evolution, and diversification of benthic octopuses (Cephalopoda: Octopodoidea)

Abstract

Here, we evaluate the so-called Thorson's rule, which posits that direct-development and larger eggs are favored toward the poles in marine organisms and whose validity been the subject of considerable debate in the literature, combining an expanded phenotypic dataset encompassing 60 species of benthic octopuses with a new molecular phylogeny. Phylogenetic reconstruction shows two clades: clade 1 including species of the families Eledonidae, Megaleledonidae, Bathypolypodidae, and Enteroctopodidae, and clade 2 including species of Octopodidae. Egg size, development mode, and all environmental variables exhibited phylogenetic signal, partly due to differences between the two clades: whereas most species in clade 1 inhabit cold and deep waters, exhibit large eggs and hatchling with holobenthic development, species from clade 2 inhabit tropical-temperate and shallow waters, evolved small eggs, and generally exhibit merobenthic development. Phylogenetic regressions show that egg size exhibits a conspicuous latitudinal cline, and that both egg size and development mode vary with water temperature. Additionally, analyses suggest that egg size is constrained by body size in lineages with holobenthic development. Taken together, results suggest that the variation in egg size and development mode across benthic octopuses is adaptive and associated with water temperature, supporting Thorson's rule in these organisms.

Egg size, holobenthic, life histories, macroevolution, merobenthic, phylogeny

Egg size and development mode constitute crucial aspects of the marine organism's life-history, because they are related to traits such as larval development, dispersal range, hatchling size, juvenile size, and survival that are closely related to fitness according to life-history theory (Laptikhovsky 1998; Voight 1998; Collin 2003; Moran and McAlister 2009; Bownds et al. 2010). Consequently, most present-day geographic patterns of egg size and developmental mode have long been assumed to reflect adaptive responses of organisms to their current environment (Marshall et al. 2012). For example, many polar and deep-sea species exhibit bigger eggs and direct development (without larval stage) in comparison to species from tropical and shallow waters, suggesting that challenging environmental conditions such as low temperatures and primary productivity could select against small eggs and larvae. The logic behind this hypothesis, conceptualized as Thorson's rule (Thorson 1936; 1950; Mileikovsky 1971; Laptikhovsky 2006; Marshall et al. 2012), posits that big eggs can store a greater amount of nutrients required to survive in the deep-sea and polar habitats, in which resource availability is not sufficient for basic biological requirements in the early stages of development. While Thorson's suggestion was highly influential and seemingly supported by earlier studies, its generality has been recently called into question as new evidence across different taxa and exceptions to the rule accumulate in the literature (Young 1990; Clarke 1992; Stanwell-Smith et al. 1999), with some researchers ultimately arguing that the rule “should be laid to rest” (Pearse 1994, p. 26).

Several confounding factors underlie the current debate on the validity of Thorson's rule and, more generally, on the existence of broad biogeographical patterns in life-history strategies across marine invertebrates. Multiple studies are restricted to a limited geographic range, latitudinal trends may vary across taxa (Clarke 1992; Stanwell-Smith et al. 1999) and between geographic regions (Gallardo and Penchaszadeh 2001; Collin 2003) and, as emphasized by Marshall et al. (2012), only a handful of studies have tested observed trends with formal statistical approaches. Furthermore, many of these studies quantify the number or fraction of species that exhibit alternative reproductive strategies at the community level (e.g. Thorson 1950; Jablonski and Lutz 1983; Clarke 1992; Gallardo and Penchaszadeh 2001; Collin 2003; Fernández et al. 2009; Pappalardo and Fernández 2014), ignoring how processes at multiple timescales, such as founder effects, diversification, and recent colonization events, contribute to community composition (Thatje et al. 2005). As a consequence, even if latitudinal trends in life-history strategies happen to be quite general, unraveling the ecological and evolutionary factors that might explain them remains a challenge because species often resemble each other due to shared ancestry, environmental variables change on a correlated fashion and historical patterns of vicariance and colonization are generally unknown. Finally, several studies include reanalysis of older datasets, hence the empirical evidence for or against the rule is actually substantially more limited than a superficial inspection of the literature would suggest (Table (1)).

Table 1.

Summary of large-scale studies that discuss the validity of Thorson's rule in the literature

Taxa	Location	Latitude§	Support	Stats	Variable§§	Predictor	Reference	Data§§§
Prosobranch gastropods	North Atlantic	28°N–75°N	Yes	No	Dev (%)		1. Thorson 1950	New
Bivalves	North Atlantic	5°N–30°N	Yes	No	Dev (%)		2. Ockelman 1965	New
Prosobranch gastropods	North Atlantic	4°N–28°N	Yes	No	Dev (%)		3. Mileikovsky 1971	New and Old (Ref. 1)
Prosobranch gastropods	Southern Ocean		Yes	No	Dev (#)		4. Picken 1980	New
Prosobranch gastropods and bivalves	North Atlantic	4°N–28°N	Yes	Yes	Dev (%)	Lat	5. Jablonski and Lutz 1983	New and Old (Ref. 2)
Prosobranch gastropods	North Atlantic	4°N–28°N	Yes	Yes	Dev (%)	Lat	6. Clarke 1992	Old (Refs. 1, 2, 6)
Holothuroids, asteroids and ophiuroids	Southern Ocean		No	No	Dev (%)		7. Pearse 1994	New
Octopuses	Worldwide	50°N–50°S	No	Yes	Egg (sp)		8. Voight 1998	New
Invertebrates	Southern Ocean		No	No	Dev (%)		9. Stanwell-Smith et al. 1999	New
Prosobranch gastropods	SE Pacific vs SW Atlantic	20°S–55°S	Mixed	Yes	Dev (%)	Lat	10. Gallardo and Penchazadhe 2001	New and Old (Ref. 1)
Calyptreid gastropods	Worldwide	50°N–50°S	Yes	Yes	Dev (%)	Lat	11. Collin 2003	New
Bivalves, gastropods and chitons	SE Pacific	18°S–54°S	Yes	Yes	Dev (#)	Lat, Temp	12. Fernandez et al. 2009	New
Anomurans, brachyurans and peracarids	SE Pacific	18°S–54°S	Yes	Yes	Dev (#)	Lat, Temp
Anomurans, brachyurans and peracarids			Yes	Yes		Lat, Temp
Polychaetes	Worldwide	50°N–50°S	Mixed	Yes	Egg (sp)	Temp, Prod	13. Marshall et al. 2012	New
Holothuroids and ophiuroids	Worldwide	50°N–50°S	Yes	Yes	Egg (sp)	Temp, Prod
Calyptreids and saccoglosan gastropods	Worldwide	50°N–50°S	Mixed	Yes	Egg (sp)	Temp, Prod
Bivalves, gastropods and chitons	SE Pacific vs SW Atlantic	18°S–54°S	Yes	Yes	Dev (#, %)	Lat, Temp	14. Pappalardo and Fernandez 2014	New and Old (Refs. 1, 2, 3, 11, 12)
Muricid gastropods	Worldwide	50°N–50°S	Yes	Yes	Dev (sp)	Temp	15. Pappalardo et al. 2014	New

§

Approximate latitudinal range because the exact range is not always reported for each specific taxonomic group.

§§

Response variable, development mode or egg size (% = fraction of the community, # = community-wide species richness, sp = species-specific traits).

§§§

New versus Old refer only to the papers included in the table, to illustrate how results in the literature are partly redundant because data from earlier studies are included in more recent analyses.

Many of these limitations can be readily circumvented with the use of modern phylogenetic statistical tools, and yet there is a paucity of biogeographic studies on marine invertebrate life-history evolution in an explicit phylogenetic context (but see Marshall et al. 2012; Pappalardo et al. 2014). To our knowledge, there are virtually no studies that attempt to disentangle how sets of life-history traits (i.e. reproductive strategies) vary with different environmental variables after controlling for phylogeny, which would shed light not only on general geographic or bathymetric patterns but also on the selective pressures involved. From the original formulation to current interpretations of the rule, for instance, egg size, and development mode are proposed to vary with latitude and water depth due to their impact on water temperature, primary productivity or both (Laptikhovsky 2006: Marshall et al. 2012). To bridge this gap, here we analyze how egg size and development mode vary across benthic octopus lineages in relation to environmental conditions. Benthic octopuses are particularly suitable to address this question because they are widely distributed across the globe, inhabiting deep, and shallow waters with contrasting temperatures and productivity levels. In addition, they exhibit two main life-history strategies: on the one hand, a holobenthic life cycle involving the production of few large eggs that result in well-developed benthic hatchlings (Sweeney et al. 1992; Villanueva and Norman 2008), and on the other hand, a merobenthic strategy (sensu Boletzky 1992) that involves the production of numerous small eggs hatching into free-swimming planktonic paralarvae (Sweeney et al. 1992; Villanueva and Norman 2008). As a result of these contrasting strategies, egg sizes across benthic octopus species range from 1 mm to nearly 40 mm, and fecundity estimates range from 30 to more than 500,000 eggs (Norman 2000; Voight and Drazen 2004).

Our main research questions are: Does egg size and development mode of benthic octopuses follow Thorson's rule? More specifically, how are changes in egg size and development mode related to body size, latitude, and environmental conditions? And finally, how do these results relate to patterns of diversification of benthic octopus lineages inferred from the fossil record? To answer these questions, we first reconstructed the phylogeny of benthic octopuses based on multiple molecular markers and including a large subset of species that remain poorly studied. Subsequently, we combined this information with an updated database of egg size and developmental mode in these organisms, concomitantly with environmental variables such as latitude, bathymetry, temperature, and productivity, to evaluate Thorson's rule within an explicitly evolutionary framework.

Materials and Methods

PHENOTYPIC DATA AND ENVIRONMENTAL VARIABLES

We included a total of 60 octopus species for which gene sequences and information on development mode (merobenthic or holobenthic), maximum egg size (mm), which corresponds to the length of either mature ovarian oocytes or fertilized eggs of reproductive females (these measures are highly correlated for the few species in which both estimates are available, Pearson r = 0.95), and maximum mantle length (mm), was available. For environmental variables, we collected latitude (mean point) and bathymetric distribution (m) from the literature (Dataset S1). For latitude, the mean-point method estimate was used based on Stevens (1989); specifically, the absolute difference between the maximum and minimum value of the latitudinal range was calculated (Fig. S1). Complementarily, habitat water temperature (°C) was obtained from the mean point latitude and maximum depth habitat of each octopus species from the World Ocean Atlas 2013 (Locarnini et al. 2013); we employed the average value of six decades of observations encompassing the following time periods: 1955–1964, 1965–1974, 1975–1984, 1985–1994, 1995–2004, and 2005–2012. Surface data on chlorophyll-a (CHL-A, mg/m³), obtained from the Bio-Oracle database (Tyberghein et al. 2012, www.oracle.ugent.be/) via the DIVA-GIS software (Hijmans et al. 2001), was used as a proxy of productivity, since it has been recognized as an appropriate estimation of maximum photosynthetic rate, and it is shared among all primary producers (Hout et al. 2007).

PHYLOGENETIC RECONSTRUCTION

To reconstruct the phylogeny of the octopus species in our dataset, we first obtained a database of sequences from GenBank for three mitochondrial genes–16S rRNA (16S), Cytochrome Oxidase I (COI), and Cytochrome Oxidase III (COIII)–and one nuclear gene–Rhodopsin (RHO)–for a total of 52 species. In addition, we obtained and added the COI gene sequence for another eight species in the database (Dataset S1), including Octopus vulgaris sensu stricto from France (Amor et al. 2017), with the following protocol. Total DNA was extracted following the saline extraction protocol (Aljanabi and Martinez 1997). The final volume was 25 μl of reaction, containing the following concentrations, 1 × of buffer 10 × (200 mM Tris-HCL (pH 8.4), 500 mM KCL), 2 mM of 50 mM MgCl₂, 0.4 μg/μl of 10 mg/ml BSA, 0.1 μM of each primer (see primers in Allcock et al. 2008), 200 mM of dNTPs, 0.06 U/μl of DNA taq polymerase (Invitrogen) and 2 ng/μl of DNA. After an initial denaturation (1 min at 94°C), the reaction mixtures were subjected to five cycles of 94°C (40 s); 45°C (40 s) and 72°C (60 s). The second phase was composed of 35 cycles of 94°C (40 s); 51°C (40 s), 72°C (60 s) followed by a final extension at 72°C (7 min), using a thermocycler. Double-stranded PCR products were purified and sequenced in both directions using an Automatic Sequencer at Macrogen, Inc. (Seoul, South Korea). Sequences were edited and aligned with Muscle implemented in MEGA 7.0 (Kumar et al. 2016).

Phylogenetic relationships were inferred from a partition matrix (16S, COI+COIII, RHO) (each gene with different substitution model) determined using PartitionFinder v2 (Lanfear et al. 2012) according to the Bayesian Information Criterion (BIC) with the nucleotide alignment divided into three subsets: 16S (GTR+G+I), COI + COIII (GTR+G+I), RHO (HKY85+G+I). This phylogenetic analysis was compared against a matrix including the concatenated dataset (16S+COI+COIII+RHO) using GTR+G+I substitution model by means of Bayes Factor (BF, Kass and Raftery, 1995), implemented in the software Tracer v1.5 (Rambaut and Drummond 2009). Subsequently, we used a phylogenetic hypothesis based on a Bayesian framework using Mr. Bayes v3.2 (Ronquist et al. 2012) fit to the models and partitions identified in the PartitionFinder analysis. This allowed us to include phylogenetic uncertainty in the subsequent comparative analyses. By means of MCMCMC (Metropolis Coupling Monte Carlo Markov Chains), four chains were run using 10,000,000 iterations sampling parameters and trees every 1000 iterations. This analysis was performed at least five times to check the convergence of the chains. The trees were rooted using Vampyroteuthis infernalis (Vampyromorpha) as an outgroup, since this species has been previously described as sister group of the Octopodiformes (Young et al. 1998; Carlini et al. 2001). Finally, after checking the posterior probabilities of trees in BayesTrees v1.3 (Meade 2011), likelihood scores and Effective Sample Size (ESS) were obtained in Tracer v1.5 (Rambaut and Drummond 2009), and the first 1000 trees (10%) were discarded as burn-in.

STATISTICAL ANALYSES

We employed three sets of analysis to estimate how phylogeny and environmental variables relate to variation in egg size and development mode. First, we ran univariate analyses to estimate phylogenetic signal (Pagel's λ) in each of the phenotypic and environmental traits included in our study (Pagel 2002). For quantitative variables, we employed Bayesian inference to estimate parameter λ with the BayesTraits v3 software (Meade and Pagel 2017). This parameter ranges between 0 and 1 and quantifies the amount of phylogenetic signal in the studied trait: when λ = 0 the trait distribution across species is independent of their phylogeny, whereas a λ = 1 indicates that the trait distribution closely matches the expectation according to Brownian motion (Pagel 1997, 1999, 2002). For this analysis, we employed raw latitude, and log₁₀-transformed egg size, mantle length, depth, temperature, and chlorophyll-a, which exhibited a log-normal distribution(Fig. (1)). To determine if the observed signal differs from zero, λ estimates were compared with the marginal likelihoods of λ forced to zero employing the Bayes Factor implemented in the software Tracer. For development mode, which was the only qualitative variable in our dataset (holobenthic vs merobenthic), we employed the phylo.signal.disc algorithm by E. L. Rezende, which compares the number of transitions according to unrestrained parsimony against a null distribution obtained after randomizing the tip data, which should break any phylogenetic structure (Rezende and Diniz-Filho 2012). Our null distribution was based on 1000 replicates, implemented in R v.3.4.4 (R development core 2018).

Figure 1

Phylogeny of benthic octopuses and the distribution of phenotypic traits and environmental data across the studied species (n = 60). All variables are quantitative with the exception of development mode (gray = holobenthic, white = merobenthic), with the size of the circles being proportional to mantle length and egg size.

Second, the association between egg size, development mode, body size (mantle length), and environmental variables was estimated employing Bayesian inference and phylogenetic generalized least squares (PGLS) regressions (Pagel 1997). To account for phylogenetic signal, λ was estimated concomitantly with parameters from the regression model (Pagel 2002). For this analysis, we originally included egg size as the dependent variable and mantle length and development mode as independent variables. Given the limited variance explained by this model, we ran another multiple PGLS regression of egg size versus the different environmental variables (latitude, depth, temperature, and chlorophyll-a), employing absolute instead of raw latitude in this case. Multicollinearity between environmental variables was not an issue according to Pearson pairwise correlations, which ranged between r = –0.366 and 0.571 (n = 60), and variance inflation factor (VIF < 5.0) (Zuur et al. 2010) calculated with the module Regression in BayesTraits. Bayesian analyses were run for 10,000,000 iterations by means of MCMC, sampling parameters every 1000 iterations, eliminating the first 10% of parameters as burn-in, over 500 random trees previously obtained from MrBayes. The 95% highest posterior density interval (HPD) for each parameter was calculated in Tracer software, and we also report the fraction of samples falling above or below 0 (depending on the sign of the mean estimate) as an index of statistical significance analogous to the P-value of frequentist approaches.

Third, we tested how development mode (holobenthic vs merobenthic) is affected by environmental variables employing the likelihood phylogenetic multivariate logistic regression method (PLogReg) proposed by Ives and Garland (2010). This analysis was implemented with the phylolm package (Ho and Ané 2014) executed in R. Confidence intervals (95% CI) were estimated with 1000 iterations of parametric bootstrap. The outgroup V. infernalis was excluded in all analyses described above and, for simplicity, we make no distinction between slope estimates of PGLS and logistic regression and refer to them as β.

Results

The phylogeny resulting from the partition matrix (16S, COI+COIII, RHO) (log₁₀ BF > 5) was better than the another derived from the concatenated matrix based on likelihood scores and ESS >400 (Figs. S2 and S3). The phylogeny from the partitioned dataset is consistent with the current octopus taxonomy (Fig. (1)). The consensus of 9000 phylogenetic trees from MrBayes showed high posterior probability values (>0.95) for most of the nodes. From these trees, two major groups inhabiting very different environments, and exhibiting contrasting life-history strategies, were obtained (Fig. (1)). The first group (hereafter “clade 1″) includes primarily cold water and deep-sea species from the families Bathypolypodidae, Eledonidae, Enteroctopodidae, Megaleledonidae. The second group (“clade 2″) encompasses tropical-temperate and shallow water species from the family Octopodidae (Fig. (1)). Species from clade 1 exhibit substantially large eggs (10.0–17.5 mm) in comparison to the smaller eggs in species from clade 2 (2.9–10.6 mm), and are primarily holobenthic (Fig. (1)).

Partly due to this clear dichotomy between clades, phylogenetic signal was moderate to high in virtually all phenotypic and environmental traits (Table S1). While λ was intermediate for egg size and chlorophyll-a (λ ∼ 0.5), it was considerably high (λ > 0.8) for mantle length and the remaining environmental variables (Table S1). According to 95% HPD and Bayes Factor estimates, phylogenetic signal was higher than 0 for all variables analyzed (BF > 0.5, Table S1). Similarly, development mode also exhibit highly significant signal according to the based on the phylo.signal.disc algorithm by Rezende (P < 0.001; Fig. S4): while 11 transitions are required to obtain the observed distribution of development modes along the octopus phylogeny (Fig. (1)), the median number of transitions across the 1000 replicates employed to build the null distribution was 20 (range between 12–27 transitions). Taken together, these analyses show that there was a high degree of concordance between our phylogenetic reconstruction, life-history strategies, and environmental correlates.

The phylogenetic regression of egg size against mantle length (β = –0.11, 95% HPD between –0.50 and 0.0.29) and development mode (β = –1.04, 95% HPD between –2.26 and 0.16) shows that neither main effects differ statistically from zero, the later possibly reflecting reduced statistical power associated with a highly clustered distribution along the phylogeny (Fig. (1)). Interestingly, however, the interaction between mantle length and development mode was higher than zero (β = 0.68, 95% HPD between 0.12 and 1.25), with the fraction of iterations falling below this value giving rise to an estimated P = 0.011, suggesting that egg size is positively associated with mantle length only in species with holobenthic development (Fig. (2)). Nonetheless, the model fit was moderate (median R² = 0.276) and the resulting λ generally high (median λ = 0.62) and comparable to the univariate estimate for egg size, and therefore a substantial fraction of the residual variation in this trait remains unexplained in this model.

Figure 2

Egg size and development mode. Left. Association between egg size and mantle length as a proxy of body size in benthic octopuses, showing the statistically significant mantle length × development mode interaction (P-value calculated from the parameter distribution across iterations; see Methods). Right. Egg size in relation to development type (gray = holobenthic, white = merobenthic).

Consequently, to facilitate interpretation, we do not control for mantle length and development mode in the PGLS regression between egg size and environmental variables. This analysis shows that, out of the four environmental predictors, only latitude exhibits a positive slope higher than zero (Table (2)), with a median estimate suggesting a 4.5-fold increase in egg size between the Equator (predicted egg size = 3.64 mm) and the highest latitude in the dataset (egg size = 16.4 mm at 63°), everything else being equal (Fig. (3)). Interestingly, when latitude is removed from the model, temperature effects become substantially stronger (β = –0.152, 95% HPD between –0.38 and 0.09) and its median estimate gives rise to a 3.2-fold difference in egg size between species inhabiting the dataset's thermal extremes (predicted egg size = 4.03 mm at 27°C and 12.9 mm at 0.5°C). Given the association between these variables (n = 60, Pearson r = –77, P << 0.0001), this analysis indicates that temperature explains to some degree the latitudinal gradient observed in the full model (Fig. (1)). In these models λ tended to be lower (median λ ∼ 0.32 in both cases), suggesting that a fraction of the residual variation in egg size originally attributed to phylogeny is associated with latitude and temperature.

Table 2.

Phylogenetic generalized least squares (PGLS) analyses to evaluate the multiple regression model between egg size of benthic octopuses and environmental gradients of their habitats

	β	r	P
Intercept	0.384 (–0.221, 1.027)
Latitude	0.0104 (0.004, 0.159)	0.461 (0.424, 0.499)	<0.001
Depth	0.082 (–0.087, 0.256)	0.168 (–0.123, 0.218)	0.177
Temperature	–0.009 (–0.240, 0.194)	–0.013 (–0.188, 0.081)	0.523
Chlorophyll-a	0.033 (–0.090, 0.162)	0.015 (–0.108, 0.194)	0.301

β is the slope from multiple regression and r is the standardized correlation coefficient. The 95% HPD are shown in parentheses, and P-values describe the fraction of iterations with the opposite sign of β, which is analogous to a frequentist test of the null expectation that β = 0.

Figure 3

Relationship between egg size and different environmental variables. Bivariate plots are shown (gray = holobenthic, white = merobenthic), whereas the regression lines represent parameter estimates from PGLS. The continuous line depicts a significant effect whereas the dotted lines are not significant (P-values were calculated from the parameter distribution across iterations; see Methods).

Finally, in the multivariate phylogenetic logistic regression analysis, only temperature was significant related to development mode (Table (3), Fig. (4)). When temperature was removed, the effects of depth became positive and statistically significant (β = 1.420, with 95% CI between 0.481 and 2.211; P = 0.041), suggesting that species from the deep sea are generally holobenthic and this effect may be mediated by differences in temperature in the water column.

Table 3.

Multiple phylogenetic logistic regression analyses to evaluate how development mode is affected by environmental variables

	β	P
Intercept	8.471 (8.116, 8.973)	0.0854
Latitude	0.037 (–0.034, 0.081)	0.3417
Depth	–1.330 (–2.145, –0.397)	0.2555
Temperature	–6.584 (–7.223, –5.610)	0.0103
Chlorophyll-a	0.077 (–0.819, 1.408)	0.9200

Values within parenthesis correspond to 95% confidence intervals, P-values correspond to estimates obtained with parametric bootstrap.

Figure 4

Relationship between development type and environmental predictors. Left. Development mode was significantly affected by temperature, with its frequency distribution shown with the histograms (gray = holobenthic, white = merobenthic). Right. Parameter estimates from the logistic regressions expressed in standardized z-transformed units. The continuous and dotted lines depict significant and nonsignificant effects, respectively (P-values obtained from parametric bootstrap).

Discussion

Our analyses are in agreement with the contention that larger sizes and direct development are favored by selection in colder environments at higher latitudes. In benthic octopuses in particular, our results also indicate a clear dichotomy between two contrasting life-history strategies that are clustered within two well-defined phylogenetic groups (Fig. (1)): whereas clade 1 encompasses species with large egg sizes and hatchlings with holobenthic development, corresponding to a K-strategy in classic life-history theory, species in clade 2 exhibit small eggs sizes and a merobenthic development typically associated with an r-strategy. Taken together, these results provide strong support to the proposition by Thorson (1936) that latitudinal gradients in egg size and development mode reported in many lineages of marine organisms emerge as a response to selective pressures favoring different; life-history strategies.

Our working phylogenetic hypothesis, based on three mitochondrial and one nuclear markers and including 60 species from five families (including specimens from the new family Enteroctopodidae), was highly resolved and is generally consistent with the taxonomy of benthic octopuses and previous phylogenetic reconstructions (Strugnell et al. 2014; Jereb et al. 2014). Recently, Strugnell et al. (2014) investigated phylogenetic relationships of 23 octopus species using four mitochondrial genes and three nuclear genes. The topology of the phylogenetic tree in the present study is similar to that of Strugnell et al. (2014) and others (e.g. Strugnell et al. 2005; Lindgren et al. 2012; Ibáñez et al. 2014). The presence of two well-defined clades with clearly distinct life-history strategies and distribution constitutes a crucial result (Fig. (1)), and yet other phylogenetic patterns deserve attention. Whereas clade 1 depicts phylogenetic relationships between four octopus families and is highly concordant with taxonomic information, clade 2 encompasses species belonging to the family Octopodidae. Within this clade, two monophyletic groups are represented, one composed of Octopus sensu stricto and the other containing species of the genera Amphioctopus and Hapalochlaena, and another two groups are paraphyletic and composed of species of Octopus sensu lato plus Abdopus, Callistoctopus, Cistopus, Macroctopus, Robsonella, and Scaeurgus (Fig. (1)). As in other studies, Octopus is polyphyletic (e.g. Strugnell et al. 2005; Lindgren et al. 2012; Ibáñez et al. 2014), which is probably related to the fact that many Octopus species are poorly described and have not been placed in a genus (sensu Norman and Hochberg 2005). Consequently, we echo the suggestion by previous authors (Gleadall 2004; Kaneko et al. 2011; Ibáñez et al. 2014) that Octopus systematics requires an in-depth deep taxonomic revision, not only for classification purposes but also to understand when, how, and hopefully why these lineages have diversified.

Having shown the high degree of concordance between egg size, developmental mode, and environmental predictors and how these variables are clustered in the phylogeny (Fig. (1)), our analyses raise the question of how the evolution of the two contrasting life-history strategies have contributed to the diversification of these lineages. The diversification of benthic octopuses was estimated to have occurred around Upper Cretaceous (Strugnell et al. 2006; Tanner et al. 2017). Based on the recently discovered an Octopodid fossil (Styletoctopus annae, Lebanon ∼33°N, Fuchs et al. 2009), ancestral octopuses likely lived around 93 Mya during the Cretaceous, in a temperate habitat that might have experienced high temperatures (i.e. sea surface temperatures of 26–35°C, Forster et al. 2007; Litter et al. 2011). The temperature of the deep sea during this period oscillated between 9 and 19°C (Huber et al. 2002), which is substantially warmer than present-day conditions. Based on these observations, our analyses and previous work (Laptikhovsky 1998; Ibáñez et al. 2014), we suggest a common ancestor of benthic octopuses exhibiting small- to medium-sized eggs characteristic of contemporary species with merobenthic development. In the Oligocene, when the Antarctic octopuses originated (Strugnell et al. 2008), the sea floor temperature was 4–5°C and dropped to roughly 2°C in the Miocene (McClain and Hardy 2010), during the radiation of deep-sea octopuses. Given the drastic global oceanographic changes (e.g. temperature, sea level, anoxic events) observed during the Cenozoic (McClain and Hardy 2010), it is plausible that the evolution of large eggs and holobenthic development followed the colonization of the deep sea. Accordingly, phylogenetic evidence suggests that polar and deep-sea octopuses originated from shallow water forms, and rapid diversification in both habitats has been attributed to Southern Ocean cooling during the Miocene (Strugnell et al. 2008, 2011; Ibáñez et al. 2016). The evolution of large eggs in this environment, concomitantly with larger benthic hatchlings with a more restricted distribution (Villanueva et al. 2016), might explain the rapid diversification of Antarctic octopuses and colonization of previously unexplored niches, and partly account for the biogeographic patterns detected in our study.

Succinctly, whereas highly stochastic environments, such as tropical and temperate shelf habitats, favored octopuses who produced numerous small eggs that hatched into free-swimming planktonic paralarvae (i.e. r-selection strategy), colder, and more stable environments favored the evolution of larger eggs and holobenthic development to maximize hatchling survival (K-selection strategy). Importantly, the effect sizes reported here, which support this evolutionary scenario, are far from trivial. While the full PGLS model (Table (2)) suggests a 4.5-fold difference in egg size between species in the equator versus those in the Arctic or the Antarctic, this effect estimated in a linear dimension translates to a 91-fold (4.5 cubed) difference in egg volume, which constitutes a more adequate estimate of maternal investment and nutrient requirements. Perhaps not surprisingly, then, mantle length was associated with egg size in holobenthic species (Fig. (2)), since large-sized species (e.g. Graneledone boreopacifica, Megaleledone setebos, Octopus conispadiceus) have larger eggs. This interaction could represent a macroevolutionary trade-off for holobenthic species since fecundity is negative related to egg size in many cephalopods (Calow 1987), once again supporting the dichotomy between r- versus K-strategists of classic life-history theory across octopus species. Interestingly, experimental evidence suggests that lower temperatures prolong embryonic development periods in octopuses and result in larger eggs (Laptikhovsky 1999; Robinson et al. 2014), hence the emergence of this evolutionary novelty might have been originally mediated by plastic responses to changing environmental conditions. On a similar vein, development rates of planktonic larvae may be highly constrained at cold temperatures, as results for the genus Enteroctopus suggest (Villanueva et al. 2016), which could partly explain the near absence of merobenthic species near the poles.

In summary, our analyses indicate that life-history variation across benthic octopuses is consistent with Thorson's rule. Even though this rule constitutes a recurrent pattern in several marine invertebrate taxa including Annelida, Echinodermata, Crustacea, and Mollusca, and seems to be related to type of development (Thorson 1936, 1950; Mileikovsky 1971; Gallardo and Penchaszadeh 2001; Collin 2003; Fernández et al. 2009; Marshall et al. 2012; Pappalardo et al. 2014), it has been frequently invoked without appropriate phylogenetically informed analyses (e.g. Gallardo and Penchaszadeh 2001; Collin 2003; Fernández et al. 2009). Nonetheless, results of large-scale analyses are often ambiguous (Table (1)), partly because widely used analytical approaches cannot disentangle how evolutionary processes and more recent ecological events contribute to community composition in different localities (Thatje et al. 2005). Here, we show that both the ecogeographic pattern and the explanation proposed by Thorson (1936) hold for benthic octopuses, describe the impact of latitude, temperature, and bathymetric depth on egg size and development mode, and discuss how the emergence of contrasting life-history r- and K-strategies may have evolved in these lineages and contributed to their diversification. Similar approaches in other lineages can shed light, therefore, not only on the generality of different biogeographic patterns, but also how they emerge in evolutionary time.

AUTHOR CONTRIBUTIONS

C.M.I, J.S., and M.C.P-G collected octopuses from different locations and performed laboratory work (e.g. PCRs). C.M.I., R.D.S., E.L.R., C.E.H, E.P., and J.A-L. carried out the phylogenetic reconstruction and the statistical analyses. All authors contributed with manuscript writing and approved the final draft.

ACKNOWLEDGMENTS

We thank Claudio González, Unai Markaida, Cesar Salinas, and Arminda Rebollo for their help with octopus tissue samples. We also thank Samantha Price, Janet R. Voight, and five anonymous reviewers for their comments on an early version of the manuscript. This work was partially funded by FONDECYT #3110152 grant to C.M. Ibáñez, and by CAPES FB 0002–2014 to E.L. Rezende.

DATA ARCHIVING

Sequences generated in this study are under GenBank codes (Table S1). The doi for our data is https://doi.org/10.5061/dryad.07s9f96.

CONFLICT OF INTEREST

The authors declare that they have no conflict of interest.

Associate Editor: J. Huelsenbeck

Handling Editor: Mohamed A. F. Noor