Evolutionary drivers of the hump-shaped latitudinal gradient of benthic polychaete species richness along the Southeastern Pacific coast

Introduction

Latitudinal diversity gradients (LDG) have been studied for over two centuries and the search for patterns and their processes remain an active topic in biogeography, macroecology and evolution (Pianka, 1966; Willig, Kaufman & Stevens, 2003; Hillebrand, 2004a; Mittelbach et al., 2007; Tittensor et al., 2010; Kinlock et al., 2018; Pontarp et al., 2019). Several meta-analyses have shown the inter-hemispheric consistency of a monotonic increase in species richness from high latitudes to the tropics, producing a unimodal (canonical) pattern on a global scale (i.e., Willig, Kaufman & Stevens, 2003; Hillebrand, 2004a, b; Costello & Chaudhary, 2017; Worm & Tittensor, 2018; Kinlock et al., 2018; Menegotto & Rangel, 2018, 2019; Rivadeneira & Poore, 2020). However, there is considerable debate regarding the underlying factors that determine the LDG (Rohde, 1992; Willig, Kaufman & Stevens, 2003; Mittelbach et al., 2007; Pontarp et al., 2019).

Synthesis achieved during the last decades recognize that explanations for the canonical LDG (i.e., increasing species richness towards the tropics) can be broadly classified into three major hypotheses categories (Mittelbach et al., 2007; Pontarp et al., 2019): ‘ecological limits’, ‘diversification dynamics’ and ‘time for species accumulation’. The category of ‘ecological limits’ refers to explanations based on present-day environmental conditions via higher productivity, carrying capacity, species coexistence, or niche breadth (Currie, 1991; Rosenzweig, 1995; Evans, Warren & Gaston, 2005). The ‘diversification dynamics’ hypotheses suggest that the LDG are generated by latitudinal variation in diversification rates, with greater speciation and/or lower extinction rates in the tropics (Evans, Warren & Gaston, 2005; Weir & Schluter, 2007). The ‘time for species accumulation’ hypotheses explanations invoke the tropics as sources of evolutionary novelties, where clades originate making the LDG the result of dispersal limitation towards subtropical zones (Wiens & Donoghue, 2004; Wiens & Graham, 2005; Jablonski et al., 2013).

In spite of its apparent generality, the canonical LDG pattern is not matched by many continental and marine taxa (see Platnick, 1991; Kindlmann, Schödelbauerová & Dixon, 2007; Tittensor et al. 2010; Worm & Tittensor, 2018, Kinlock et al., 2018). A growing body of evidence shows that non-canonical LDG (i.e., inverse or hump-shaped trends) are common among marine organisms at a global scale (Chaudhary, Saeedi & Costello, 2016; Woolley et al., 2016; Rivadeneira & Poore, 2020). This same phenomenon has been reported along the Southeastern Pacific coast (0–56° S, hereafter SEP) and also involve widely different taxa such as macroalgae (Santelices, 1982; Santelices & Marquet, 1998; Santelices, Bolton & Meneses, 2009), sponges (Desqueyroux & Moyano, 1987), anthozoans (Häussermann & Försterra, 2005), bryozoans (Moyano, 1991), polyplacophorans (Navarrete et al., 2020), gastropod and bivalves (Herm, 1969; Valdovinos, Navarrete & Marquet, 2003), polychaetes (Lancellotti & Vásquez, 2000; Hernández, Moreno & Rozbaczylo, 2005; Moreno et al., 2006) and different groups of crustaceans (Fernández et al., 2009; Rivadeneira et al., 2011). These inverse trends along the SEP have been shown to not be the result of sampling artifacts (see Valdovinos, Navarrete & Marquet, 2003; Rivadeneira et al., 2011).

In this context, the biogeography of benthic polychaetes from the SEP has historically received little attention and available studies have been mainly descriptive (see Lancellotti & Vásquez, 1999, 2000). Indeed, a significant gap in knowledge in the whole Eastern Pacific and Chile in particular, is recognized on a recent global study of polychaete biogeography (Pamungkas, Glasby & Costello, 2021). However, Hernández, Moreno & Rozbaczylo (2005) and Moreno et al. (2006) evaluated the latitudinal change in species richness of benthic polychaetes along the coast of Chile (between Arica at 18° S and Cape Horn at 56° S), finding a hump-shaped LDG with a maximum peak around Chiloé Island (42° S) in southern Chile. These authors proposed that this pattern might be determined by geometric constraints and historical events such as eustatic movements, cyclic effects of El Niño-Southern Oscillation, shallowing of the oxygen minimum zone and glacial advances and retreats, all of which have mainly occurred since the Neogene (i.e., during the last 23 Mya). Although these studies represented a significant advance in understanding the mechanisms that may generate the hump-shaped LDG, their analyses used a low spatial resolution (3° latitude bands) to register species occurrences and constrained the study area to the geopolitical limits of Chile (an arbitrary decision shared by the majority of studies in other taxa. i.e., Desqueyroux & Moyano, 1987; Moyano, 1991; Fernández et al., 2009; Rivadeneira et al., 2011; Lee & Riveros, 2012). In these studies (Hernández, Moreno & Rozbaczylo, 2005; Moreno et al., 2006), therefore, the exclusion of a vast portion of the Warm Temperate Southeastern Pacific biogeographic province (sensu Spalding et al., 2007) could generate spurious results derived from truncating the real biogeographic limits of these polychaete species.

Here, we reappraise the patterns and processes related to the LDG reported for benthic polychaetes in Chile, using a greater spatial extent (3–56° S) and resolution (0.5° latitude) than prior analyses. The northern end is defined by the Guayaquil Gulf (limit between the Panamanian and Peruvian–Chilean Provinces, Vegas-Velez, 1980; Boschi, 2000; Robertson & Cramer, 2009; Ibanez-Erquiaga et al., 2018), while the southern limit consist in the end of the South American continent isolated from Antarctica by the Circumpolar Current. Regarding the southern limit, differences in polychaete community composition have been already reported between the Magellan region and Weddell Sea shelves (Montiel et al., 2005).

The existence of non-canonical LDG along the SEP could be the result of a complex interaction of several factors, instead of a single dominant one. To give explanation to multi-variate and non-linear phenomena a statistical framework which incorporates and differentially weights each factor is required. Our approach would be to simultaneously evaluate the relative contribution of several proxies (SST, salinity, shelf area, etc.) which have been previously associated with specific hypothesis to explain LDG, using machine learning methods. The explanatory hypotheses considered can be separated into three categories: ‘ecological limits’, ‘diversification dynamics’, and ‘time for species accumulation’ (Table 1). These categories have been adapted from the conceptual framework used to study the canonical LDG (Mittelbach et al., 2007; Pontarp et al., 2019).

Regarding ‘ecological limits’, we will consider three hypotheses: (1) Seasonal coexistence (Valentine & Jablonski, 2015), which predicts that diversity will decline in areas with abiotic seasonal changes (Diversity 1/∝ SST range and Salinity range); (2) Species-area effect (Chown & Gaston, 2000; Valdovinos, Navarrete & Marquet, 2003), which predicts that diversity will be a direct function of the available area (Diversity ∝ Shelf area); and (3) Energy-dependency (Hawkins, Porter & Diniz-Filho, 2003; Evans, Warren & Gaston, 2005), which predicts that diversity will increase with more energy available in the ecosystem (Diversity ∝ Primary Productivity). In the category of ‘diversification dynamics’, we will consider two hypotheses: (1) Long term climate stability (Dynesius & Jansson, 2000), which predicts that more stable areas will accumulate more species (Diversity 1/∝ Median latitudinal range of the species) and (2) Temperature-dependent speciation (Rohde, 1992; Allen et al., 2006), which predicts a positive association between speciation rates and ambient temperature (Diversity ∝ SST). Finally, in the category of ‘time for species accumulation’, we will assess one hypothesis about Niche conservatism (Wiens & Donoghue, 2004), which predicts that species will tend to retain preferences for the ecological conditions where they originally evolved (Diversity ∝ constancy in SST, salinity, SST range, Median latitudinal range of the species, and Salinity range). As all these explanatory variables may reflect processes operating at both ecological and evolutionary timescales, we have selected them as candidate variables to build a machine learning model to explain the polychaete LDG in the SEP.

In addition, the study of phylogenetic diversity may shed light into the role of evolutionary processes shaping patterns of species richness (Davies & Buckley, 2011; Fritz & Rahbek, 2012; Eme et al., 2020). For instance, a species-rich area composed by many closely related species (i.e., low phylogenetic diversity) may imply elevated in situ speciation rates of some lineages. On the contrary, a species-rich area composed by many distantly related species (i.e., high phylogenetic diversity) may be associated to immigration influenced by filtering effect acting upon a trait phylogenetically over dispersed. There are many α-phylogenetic diversity indices (Winter, Devictor & Schweiger, 2013), that can be broadly classified into three categories that measures amount of richness, divergence and regularity (Tucker et al., 2017); therefore, the combined used of indices reflecting each one of these three aspects may be more informative of the role of evolutionary processes shaping species richness (Davies & Buckley, 2011; Fritz & Rahbek, 2012; Eme et al., 2020). Despite the fact that these phylogenetic diversity indices are typically based on calibrated molecular phylogenies, the taxonomic structure may be still be used as a coarse proxy of phylogenetic relationships of species in groups lacking robust and complete molecular phylogenies (Soul & Friedman, 2015). In fact, multiple studies in marine taxa have analyzed latitudinal trends in phylogenetic diversity using this taxonomic approach (Woodd-Walker, Ward & Clarke, 2002; Ellingsen et al., 2005; Tolimieri & Anderson, 2010; Rivadeneira et al., 2011; Azovsky, Garlitska & Chertoprud, 2016; Wu, Chen & Zhang, 2016). This opens the possibility to explore latitudinal patterns of phylogenetic diversity in benthic polychaetes along SEP and its connection with the LDG phenomenon.

Our study area encompasses a broad latitudinal gradient of oceanographic conditions spanning more than 50 degrees of latitude, from tropics to sub-polar areas, and hence it is ideal to study the role of ‘ecological limits’ on the LDG (Astorga et al., 2003; Valdovinos, Navarrete & Marquet, 2003; Fernández et al., 2009). At the same time the marine biota at SEP experienced a drastic turnover during the late Cenozoic, driven by tectonic, physiographic, climatic, and oceanographic shifts (Herm, 1969; Rivadeneira & Marquet, 2007; Kiel & Nielsen, 2010; Villafaña & Rivadeneira, 2014; Rivadeneira & Nielsen, 2017). These alterations induced latitudinal differences in the diversification trends that created the inverse LDG characterizing present-day bivalves and gastropods along the SEP (Kiel & Nielsen, 2010) and may have also driven the LDG of benthic polychaetes (Moreno et al., 2006). Thus, SEP is also an ideal region to study the role of evolutionary and historical processes (i.e., ‘diversification dynamics’, and ‘time for species accumulation’) shaping the LDG. Therefore, the specific goals of our study are: (1) to re-evaluate the existence of a hump-shaped LDG, including the entire Peruvian and Magellanic biogeographic provinces along the SEP, and (2) to jointly evaluate different hypotheses (by significance assessment of proxies), which could explain the formation of a hump-shaped LDG. Our study confirms the existence of a non-canonical LDG, robust to sampling bias, and suggests that the underlying explanations may be linked to processes occurring at evolutionary timescales.

Materials and Methods

Database

This study relies on a database which: (a) includes 643 polychaete species found between 3° S (southern Ecuador) and 56° S (Cape Horn, Chile) inhabiting the continental shelf (i.e., ≤200 m depth) (Fig. 1A), (b) summarizes information compiled from exhaustive reviews of literature (ca. 1,000 bibliographic references), museum collections and field expeditions, integrating more than 160 years of studies over the SEP, (c) was constructed using georeferenced data from geographic information systems (Datum WGS84), (d) was enriched by intensive field samplings conducted over the last decade for updating records and correcting biogeographical biases, particularly in southern Chile where diversity surveys historically had been scarce, specially through the research cruises CIMAR Fiordos 13, 14, 15, 17 and 18 (Comité Oceanográfico Nacional, Chile), and (e) was complemented with existing records from the OBIS database (OBIS, 2021). Taxonomy was cross-validated with the World Register of Marine Species (WoRMS Editorial Board, 2021). Altogether, such information currently constitutes the most complete database of benthic polychaetes in this region of the world. The reported latitudinal ranges of all species used in this study are provided in the Data S1.

We first determined the maximum and minimum latitude of geographic distribution of each species and assigned species presence into 0.5° latitudinal bands using a range-through approach (i.e., assuming a continuous geographic distribution between their range limits), as commonly done in marine macroecological studies (Roy et al., 1998; Tomašových et al., 2016). OBIS records were used to estimate the latitudinal limits of species beyond our study area along the eastern Pacific, from Antarctica to Alaska. We also examined whether differences in the sampling intensity may potentially mask the overall pattern, that is to say, whether areas of higher species richness may be the result of higher sampling effort. In order to test this potential confounding factor, we used 57 soft-bottom benthic assemblages sampled in the coastal shelf and collated from a literature survey. Since the sampling effort was not the same (i.e., different sampling devices and number of samples) we used the number of individuals reported for each species and each site to carry out a rarefaction analysis. This method estimates the expected species richness for a similar number of individuals sampled (Sanders, 1968; Gotelli & Colwell, 2001). In our case, we set E_s = 20 (the site with a smaller number of individuals). The information containing the local occurrences and abundances of 152 species in the 57 local assemblages is presented in the Data S2. We estimated the mean rarefied species richness per 0.5° latitudinal bin and correlated it with the species richness estimated with the range-though approach. In the absence of severe spatial bias in sampling intensity, the correlation between rarefied and regional-scale species richness should be positive and significant. However, patterns of species richness may not be necessarily similar across spatial scales (Gray, 2002; Rivadeneira, Fernández & Navarrete, 2002; Hillebrand, 2004a, 2004b; Edgar et al., 2017), and in fact the relationship between local and regional species richness may be non-linear, i.e., at higher regional species richness local assemblages may be saturated (Ricklefs, 1987; Srivastava, 1999; Rivadeneira, Fernández & Navarrete, 2002). The contrary trend, i.e., a linear relationship between local and regional species richness, suggests that local assemblages may be ‘open’ to the dispersal of species from the regional pool (Ricklefs, 1987). We explored this idea by carrying out a second-order polynomial OLS regression, where the significance of the terms can be used to support the saturation (i.e., significant quadratic component) or unsaturation (significant linear component, but non-significant quadratic term) of local assemblages.

Results

Our compiled data set shows that the species diversity of benthic polychaetes along the SEP presents a hump-shaped LDG (Fig. 1B). The maximum richness (299 species) occurs at 42° S, which is two-fold higher than in southern Ecuador (3° S) where richness was 152 species. A secondary peak is located ca. 53° S, and then species richness drops abruptly towards 55° S, with species richness levels similar to those observed in Perú and southern Ecuador. The evaluation of potential sampling artifacts showed that rarefied species richness in local assemblages has the same trend observed when using the range-through approach (r = 0.49, p = 0.009, n = 29 bins, Fig. 1B), and no evidence of saturation was detected, as the linear but not the quadratic term of the local vs. regional OLS regression was significant (p = 0.008 and 0.179, respectively).

The random forest model showed a high accuracy predicting the LDG (Figs. 2A and 3A), with a pseudo-r² = 0.95. The model underestimated the species richness in the most extreme latitudinal bins (Fig. 2B), but overall spatial autocorrelation was not detected at any spatial distance (Fig. 2C). Out of the seven predictors, four were significant (median latitudinal range, SST, salinity and SST range; ranked according to its importance, see Table 2). For a reference of the latitudinal variation of the selected variables see Figs. 3B, 3C. Partial dependence shows a hump-shaped response of species richness to SST (Fig. 4A), with maximum values at ca. 12 °C.

A more pronounced hump-shaped response of species richness was observed for SST range (Fig. 4B). Salinity and median latitudinal range showed a monotonic negative relationship with species richness, which reaches maximum values at the lowest salinity levels (Fig. 4C) and bins with lower median latitudinal ranges (Fig. 4D).

The three phylogenetic diversity indices showed different latitudinal trends and their correlation to species richness varied in sign and magnitude (Fig. 5). While PD_SES was positively correlated to species richness (r_Pearson = 0.51, p < 0.0001, Fig. 5A), MPD_SES did not show a significant relationship (r_Pearson = 0.04, p = 0.69, Fig. 5B), and VPD was negatively correlated with species richness (r_Pearson = −0.49, p < 0.0001, Fig. 5C). PD_SES values were not different than expected by the null models for most of the latitudinal gradient (Figs. 5A, 5B), except by areas around northern Chile and Perú (north of 20° S). Observed MDP_SES values were significantly lower than expected in two large areas centered around Chiloé (42° S), and northern Chile and Perú (ca. 15° S). For both PD_SES and MDP_SES these negative values suggested a phylogenetic clustering pattern. Blomberg’s K was only significant for the median latitudinal range, SST and primary productivity (Table 2).

Discussion

Our study supports the existence of a hump-shaped LDG for benthic polychaetes along the SEP, increasing the number of exceptions to the canonical pattern of higher species richness towards the tropics, both in the SEP and other regions (i.e., Moyano, 1991; Valdovinos, Navarrete & Marquet, 2003; Willig, Kaufman & Stevens, 2003; Hillebrand, 2004a, 2004b; Kindlmann, Schödelbauerová & Dixon, 2007; Santelices, Bolton & Meneses, 2009; Rivadeneira et al., 2011). Particularly, it is congruent with a recently published global study on latitudinal gradients of polychaetes (Pamungkas, Glasby & Costello, 2021). Our research also confirms previous results for benthic polychaetes along SEP (Hernández, Moreno & Rozbaczylo, 2005; Moreno et al., 2006) based on partial datasets of a more restricted scale and lower spatial resolution (ca., 3° of latitude), and it seems robust to possible sampling artifacts, as the overall shape of the LDG still holds after accounting for differences in sampling in local assemblages. Therefore, the recorded pattern in this study reflects a robust biogeographic pattern for benthic polychaetes in the SEP.

Our analyses reveal that the LDG of polychaete species is well explained by a reduced subset of predictors that could be attributed to processes operating mainly, but not exclusively, at evolutionary timescales (Table 1). Overall, these significant variables (median latitudinal range, SST, salinity, and SST range) are proxies for the predictions of four (Long term climate stability, Niche conservatism, Temperature-dependent speciation, and Seasonal coexistence) of the six ‘specific hypothesis’, corresponding to all three hypotheses categories proposed in previous literature (Mittelbach et al., 2007; Pontarp et al., 2019). It is important to note that all the selected predictors are also linked to ‘time for species accumulation’ category, and particularly to the hypothesis of ‘niche conservatism’ (Wiens & Donoghue, 2004). This hypothesis is independently supported by the detection of significant phylogenetic signal (Blomberg’s K) in the two most important predictors (median latitudinal range and SST). In addition, the LDG is also correlated to spatial patterns of phylogenetic diversity, measured according to two different indices (PD_SES and VPD).

In the context of a non-canonical LDG, as our case, the ‘niche conservatism’ hypothesis assumes that clades are originated outside of the tropics. Although we lack paleontological or phylogenetic information to fully test this idea, the prevalence of the median latitudinal range as the most important predictor of the species richness, being also the predictor with a high degree of phylogenetic conservatism, strongly supports the idea of an extra-tropical origin for many polychaetes species along the SEP. Areas of high species richness are associated to latitudinal bins with species with narrower latitudinal ranges, as predicted by the ‘Rapoport’s rule’ (Stevens, 1989; Hernández, Moreno & Rozbaczylo, 2005). Our results reinforce the idea that the Southern Chile fjord area, particularly the Chiloé region, may act as a source of new species of marine polychaetes along SEP creating a hotspot of diversity (Moreno et al., 2006).

A possible historical explanation for the origin of this hotspot is related with the fact that during the Last Glacial Maximum (LGM) the Patagonian Ice Sheet covered a large portion of the fjords probably pushing many species from the canals into the glacial coastal line (Davies et al., 2020). Indeed, it has been shown that the populations of the southern bull-kelp (Durvillaea antarctica) on the fjord area correspond to recent, probably post glacial, re-colonization (Fraser et al., 2010), implying the unavailability of this region during the LGM. Chiloé Island corresponds to the northern most coastal area where the ice sheet reached sea level, indeed the north-west part of the island has remained ice-free since the last 35,000 years (Davies et al., 2020). Therefore, the waters around Chiloé Island, especially in the north-west corner, might have worked as marine coastal refugia preserving species with previous more extended distributions towards the south.

In addition, the relationship between the LDG and phylogenetic diversity indices also suggest the importance of local diversification processes (i.e., origination and extinction, Davies & Buckley, 2011; Fritz & Rahbek, 2012; Eme et al., 2020). Thus, the significant phylogenetic clustering found in the fjords area (i.e., many species are phylogenetically closely related) suggests that this area may have not only acted as a potential glacial refugia, but also a hotspot of in situ diversification, as seen in the fossil record of marine mollusks of the region (Kiel & Nielsen, 2010).

Oceanographically, bidirectional dispersion from this area could be facilitated because precisely at this latitude the West Wind Drift diverges into the Humboldt and Cape Horn currents; running north and south, respectively (Strub et al., 2019). Biologically, the extent of the dispersion might be limited due to the strong niche conservatism detected on latitudinal range and median SST experienced for the species. Additionally, the West Wind Drift might work as a colonization pathway of the SEP from the western Pacific. Altogether, in the same way as the canonical LDG shows a decline in diversity while moving away from the tropics; on this case the decline occurs while moving away from Chiloé Island.

On another hand, the existence of phylogenetic clustering in areas of low species richness such as northern Chile and Perú is more consistent with a negative net diversification trend, produced by high extinction rates coupled to species sorting (via selective survival). Paleontological studies have shown the existence of high extinction rates of marine forms during the late Neogene–Pleistocene along the Peruvian and northern Chile coast (Herm, 1969; Kiel & Nielsen, 2010), which are linked to the Neogene onset of the modern Humboldt Current (Dekens, Ravelo & McCarthy, 2007). In particular, extinction rates of bivalves and marine vertebrates at that time were phylogenetically clustered (Rivadeneira & Marquet, 2007; Villafaña & Rivadeneira, 2014), likely associated to a strong environmental filtering (e.g., thermal tolerance). If similar processes affected polychaetes at SEP, a high extinction rate with higher survivorship of particular clades via environmental filtering process, may lead to an impoverished species richness highly phylogenetically clustered, as it characterizes the northern Chile and Perú regions.

Despite the prevalence of evolutionary-based explanations for the LDG of polychaetes, the lack of phylogenetic signal of the SST range and salinity suggest that processes operating at ecological timescales (i.e., ‘ecological limitations’) may be also important. In fact SST range is much higher in vast areas of Perú and northern Chile, likely related to the effect of ‘El Niño’ events, which governs the inter-annual variability of the sea water temperature (Shaffer et al., 1999). Associated with this and congruent with the selection of SST as a predictor of our model, a recent study also refers to SST as an important factor to globally shape the polychaetes LDG (Pamungkas, Glasby & Costello, 2021).

Conclusions

The SEP provides a natural laboratory for marine biogeography studies given the high diversity of shapes of the LDG in marine taxa. In addition to the existence of canonical LDG (Ojeda, Labra & Muñoz, 2000; Astorga et al., 2003) and inverse LDG (Santelices & Marquet, 1998; Fernández et al., 2009; Rivadeneira et al., 2011), our results validate the existence of a hump-shaped LDG in marine benthic polychaetes centered ca. Chiloé Island (42° S). More importantly, the same conceptual framework used to study the canonical LDG phenomenon could also be used to investigate non-canonical LDGs. The latitudinal patterns of distribution of benthic polychaetes needs to be evaluated in other geographical areas, e.g., central and northeastern Pacific coast, in order to fully understand the generality of the non-canonical LDG, and whether the underlying processes are the same in other regions. In the absence of suitable fossil record for polychaetes, molecular phylogenies of selected clades may shed more direct evidence of the importance of evolutionary processes shaping the LDG in this group. Future studies may be also evaluating the role of dispersal of taxa from other areas of the southern ocean (e.g., Glasby & Alvarez, 1999; Glasby, 2005) shaping the ‘Chiloé hotspot’. An integration of functional traits of species, and phylogenetic and phylogeographic analysis may help to test and validate the importance of ecological and evolutionary determinants of the diversity of benthic polychaetes.

Supplemental Information