1. Ecology
Download icon

Global gradients in intertidal species richness and functional groups

  1. Jakob Thyrring  Is a corresponding author
  2. Lloyd S Peck
  1. British Antarctic Survey, United Kingdom
  2. Department of Zoology, University of British Columbia, Canada
  3. Arctic Research Centre, Department of Bioscience, Aarhus University, Denmark
  4. Homerton College, University of Cambridge, United Kingdom
  5. Marine Ecology, Department of Bioscience, Aarhus University, Denmark
Research Article
  • Cited 1
  • Views 400
  • Annotations
Cite this article as: eLife 2021;10:e64541 doi: 10.7554/eLife.64541

Abstract

Whether global latitudinal diversity gradients exist in rocky intertidal α-diversity and across functional groups remains unknown. Using literature data from 433 intertidal sites, we investigated α-diversity patterns across 155° of latitude, and whether local-scale or global-scale structuring processes control α-diversity. We, furthermore, investigated how the relative composition of functional groups changes with latitude. α-Diversity differed among hemispheres with a mid-latitudinal peak in the north, and a non-significant unimodal pattern in the south, but there was no support for a tropical-to-polar decrease in α-diversity. Although global-scale drivers had no discernible effect, the local-scale drivers significantly affected α-diversity, and our results reveal that latitudinal diversity gradients are outweighed by local processes. In contrast to α-diversity patterns, species richness of three functional groups (predators, grazers, and suspension feeders) declined with latitude, coinciding with an inverse gradient in algae. Polar and tropical intertidal data were sparse, and more sampling is required to improve knowledge of marine biodiversity.

Introduction

The latitudinal diversity gradient in species richness across ecosystems and various functional groups has been a major research topic that has intrigued scientists since at least Darwin, 1859 and Wallace, 1878. Over time, many hypotheses have been proposed to explain this seemingly general ecological pattern in marine and terrestrial ecosystems. Ecological and physical hypotheses dominate the discussions and include drivers such as habitat area, stability, speciation rates, energy availability, and temperature (Willig et al., 2003; Clarke and Gaston, 2006; Edgar et al., 2017). However, the geographic, functional, and taxonomic generality of latitudinal diversity gradients remain lively debated as unimodal, bimodal, and inverse gradients emerge across clades, habitats, and latitudes (Rivadeneira et al., 2002; Waller, 2008; Chaudhary et al., 2016; Kinlock et al., 2018). In the marine realm, latitudinal gradients in some groups have been shown to be closely related to oceanographic covariates, such as water temperature (Roy et al., 2000), yet mammal richness peaks at high latitudes (Grady et al., 2019), and more species have been reported from polar soft-sediment habitats than at many lower latitudes (Vause et al., 2019).

While gradients in species diversity have received most attention, latitudinal changes across different functional groups (evaluated here by food acquisition, see Materials and methods) remain less studied despite their importance for ecosystem functioning. Macroalgal canopies shelter understory species from environmental stress (Krause-Jensen et al., 2016; Sejr et al., 2021), creating protective microhabitats that increase environmental heterogeneity and biodiversity, thereby maintaining a diversified understory assemblages (Bulleri et al., 2002; Watt and Scrosati, 2013; Piazzi et al., 2018). Suspension feeders are important benthic-pelagic energy couplers (Gili and Coma, 1998), and predation is widely accepted as a central structuring process in the composition and abundance of species (Vermeij, 1987; Stanley, 2008). Latitudinal studies have shown that suspension feeders dominate benthic systems in fully marine environments at high latitudes (Gili and Coma, 1998; De Broyer et al., 2014), while diversity of coastal macroalgae peaks at mid-latitudes (Keith et al., 2014), and predation pressure decreases with latitude and depth from shallow shelves to deep oceans (Taylor and Taylor, 1977; Harper and Peck, 2016). Most functional group studies have focused on a narrow set of taxa, and global-scale assemblage-wide investigations encompassing both hemispheres are rare. Thus, studies demonstrating global patterns in various functional groups are needed to understand spatial patterns, biological interactions, and ecosystem resilience to climate change.

Intertidal shores rank as one of the most studied marine habitats and are often seen as harbingers for the effects of climate change and invasive species. Regional-scale intertidal studies have found richness gradients of gastropods along coastlines in the eastern Pacific Ocean (Rivadeneira et al., 2015; Fenberg and Rivadeneira, 2019). However, no latitudinal diversity gradient of gastropods (Miloslavich et al., 2013) or macroalgae (Konar et al., 2010) was found on a global scale across oceans, and assemblage-wide studies have found missing (Blanchette et al., 2008; Cruz-Motta et al., 2010) or inverse (Griffiths and Waller, 2016) gradients. Conflicting and missing gradients suggest that richness is determined by regional or local (and not global) scale processes through biological interactions and small scale overlapping environmental gradients. For example, high habitat heterogeneity creates environmental stress mosaics that are more important in shaping biological patterns than latitudinal environmental gradients (Helmuth et al., 2006; Jurgens and Gaylord, 2018). However, the large-scale intertidal studies necessary to evaluate global-scale patterns and processes are generally missing (but see Cruz-Motta et al., 2010; Miloslavich et al., 2013).

Global assessments of intertidal biodiversity have been hindered by data scarcity from polar intertidal shores, and a lack of estimates of intertidal geographic areas. In fact, only the extent of intertidal mud flats have been quantified on a large scale (Murray et al., 2019). However, in recent decades intertidal diversity data from both polar regions have advanced (Weslawski et al., 2010; Griffiths and Waller, 2016; Thyrring et al., 2017; Thyrring et al., 2021; Peck, 2018), providing novel opportunities to study biodiversity patterns and variation in specific functional groups on a worldwide scale. In this study, we examined latitudinal diversity gradients in local intertidal species richness (termed α-diversity) and four functional groups from 433 locations throughout the northern and southern hemispheres, from the high Arctic to the Antarctic. Intertidal ecosystems comprises a variety of habitats ranging from tidal flats to mangrove forests, each supporting unique communities and taxa. Therefore, to control for habitat effects, we specifically focus on rocky shores and boulder fields, as they represent the most studied intertidal habitats and are present in all oceans. We aim to examine patterns in intertidal α-diversity and functional groups across 155° of latitude, and test the following hypotheses: (1) α-diversity decreases with increasing latitude, as evident in other ecosystems and (2) biodiversity is predicted by local-scale structuring processes rather than global-scale oceanographic features.

Results

Latitudinal α-diversity gradients

We extracted data from 433 intertidal sites between 74.8°S and 80.5°N (Figure 1). Only 11 sites were located at latitudes above 70°, south and north combined (Supplementary file 1). A generalized additive mixed effect model (GAMM) indicated a non-linear relationship among α-diversity and latitude in the northern hemisphere where α-diversity peaked at mid-latitudes between 30 and 55°N (GAMM, edf = 2.005, p=0.023). While no latitudinal scale environmental driver was significant (Table 1), α-diversity was similar in the tropical and Arctic regions (Figure 2). In the southern hemisphere, a generalized linear mixed effect model (GLMM) indicated that α-diversity was not related to any latitudinal scale environmental drivers (Table 1), and displayed a linear latitudinal pattern with a non-significant slight decline at the highest latitude (Table 1; Figure 2).

Locations of rocky intertidal sampling sites.

Some sites are not visible because of close proximity. Functional diversity data were available in green sites.

Latitudinal patterns in rocky intertidal α-diversity plotted against latitude.

Data are split into southern and northern hemispheres. A linear regression line (southern hemisphere) and a best-fit locally weighted scatterplot smoother (northern hemisphere) was added with 95% confidence intervals to aid visual interpretation.

Table 1
Results of the mixed effect models for rocky intertidal α-diversity.

The estimate, standard error (SE), z-value, and p-value are presented for each variable. The estimated degree of freedom (ed), and p-value are presented for the smoother in the northern hemisphere*.

EstimatesSEz-Valuep-Value
Northern hemisphere
(Intercept)3.0860.28310.92<0.001
Chlorophyll a−0.0040.0880.050.959
Phosphate−0.1630.288−0.560.573
Sea temperature0.0240.0151.640.101
Southern hemisphere
(Intercept)3.1990.4437.225<0.001
Latitude0.0060.0070.8990.369
Chlorophyll a−0.0060.007−0.8290.407
Phosphate−0.2400.234−1.0290.303
Sea temperature −0.0100.016−0.6530.514
  1. *Approximate significance of cubic spline regression smoother: Estimated df = 2.005, p-value=0.023.

Local-scale models revealed significant relationships between α-diversity and the four local-scale variates ice scour, macroalgal canopy cover, salinity, and wave exposure (Table 2). We focused exclusively on canopy-forming algae because canopies provide living space and protection from predation and extreme temperatures, thereby increasing α-diversity and coverage of understory organisms. A positive relationship was identified between α-diversity and salinity and macroalgal cover (Figure 3a,c), with salinity having the strongest effect (R2 = 0.19, 95% CI = 4.82–8.52, Figure 4). Wave exposure and ice scour had negative effects on α-diversity (Figure 3b,d). The non-linear relationship between wave exposure and α-diversity showed that high levels of wave exposure were necessary to affect levels of α-diversity (Figure 3b). Fewer species inhabited ice-scoured shores compared to ice-free shores (95% CI = −0.57 – −0.23, Table 2), although the model only explained around 5% (R2 = 0.05) of the variation due to large differences in the number of species on both ice-scoured and ice-free shores (Figure 3c).

Local-scale relationships between (a) salinity, (b) wave exposure, (c) ice scour, and (d) macroalgal canopy cover on α-diversity.

A best-fit locally weighted scatterplot smoother (LOESS) and 95% confidence intervals (panel b) and boxplots (panels c,d) has been added to aid visual interpretation.

Local-scale mean effect sizes and direction of ice scour, macroalgal cover, and salinity on α-diversity estimated from individual models.

Significance of regression parameters is identified as bootstrapped 95% confidence intervals (error bars) not crossing zero (* indicates significance).

Table 2
Local-scale model summaries with individual models indicated for each model.

Estimated parameters, standard error (SE), bootstrapped 95% confidence intervals (95% CI), z-values, and p-values are reported for the relationship between α-diversity and environmental covariates.

Salinity
R2 = 0.19 (GLM)
EstimatesSE95% CIz-Valuep-Value
(Intercept)−3.2010.891−4.94; −1.47−3.6<0.0001
Salinity6.6760.9464.82; 8.527.06<0.0001
Wave exposure (GAM)
R2 = 0.12
EstimatesSEz-Valuep-Value
(Intercept)*2.4350.02982.4<0.0001
Ice scour (GLM)
R2 = 0.05
EstimatesSE95% CIz-Valuep-Value
(Intercept)2.3530.0552.24; 2.4642.7<0.0001
Ice scour−0.4010.087−0.57; −0.23−4.59<0.0001
Macroalgal cover (GLM)
R2 = 0.14
EstimatesSE95% CIz-Valuep-Value
(Intercept)6.1330.3865.37; 6.8915.9<0.0001
Cover−2.350.547−3.42; −1.28−4.3<0.0001
  1. *Approximate significance of locally weighted scatterplot wave exposure smoother: Estimated df = 8.40, p-value<0.0001.

Latitudinal gradients in functional groups

Four studied functional groups (algae, grazers, predators, and suspension feeders) displayed different latitudinal patterns (Figure 5). Predators, grazers, and suspension feeders declined with latitude in both hemispheres, coinciding with an increase in algal richness (Figure 5). In general, the relative distribution of the functional groups changed faster across latitudes in the southern hemisphere where rapid changes were observed around 20°S, while northern hemisphere communities were more consistent from 25 to 60°N (Figure 5). We were unable to describe functional groups at latitudes 15–25°N as no data were available (indicated as white bars in Figure 5).

Latitudinal variation in the relative number of species from four functional groups in the rocky intertidal between 74.8°S and 80.5°N, where negative values denote southern hemisphere latitudes.

No data were available within the 5° latitudinal bands where bars are missing.

In the southern hemisphere, the proportion of predators and grazers declined beyond around 20°S, together forming 6.9% (±6.6 s.d.) of species accounted for between 20 and 55°S, but constituting 18.3% (±15.1 s.d.) of species in communities at latitudes less than 20°S (Figure 5). Suspension feeders were common in low- to mid-latitudes, constituting 34.6% (±9.6 s.d.) of the species present between 20 and 45°S (Figure 5). The relative proportion of algal species increased with latitude, with a small decline at sites between 60 and 65°S (Figure 5). In the northern hemisphere, predators were abundant at low- and mid-latitudes, accounting on average for 18.2% (±5.8 s.d.) of the populations from 0 to 60°N, before declining to 4.7% (±4.1 s.d.) at high latitudes between 61 and 80°N. The number of grazers decreased markedly from around 15°N, but stabilized at 25.4% (±6.4 s.d.) across 25°–80°N (Figure 5). Suspension feeders were the dominant functional group on shores in the northern hemisphere, especially between 25 and 60°N, but decreased to 21.4% (±1.8 s.d.) at high latitudes between 60 and 80°N (Figure 5). Algal species was the most dominant group at high latitudes (Figure 5). In both hemispheres, the assemblage composition at the highest latitudes was substantially different from elsewhere; on these high polar shores, suspension feeders were absent and only a few predators were present, while grazers and algae dominated.

Discussion

Latitudinal α-diversity gradients

Latitudinal diversity gradients have been recognized for centuries, and the generality of this phenomenon has been investigated ever since. Recent state-of-the-art meta-analyses have provided support for latitudinal gradients in various marine and terrestrial ecosystems (Hillebrand, 2004; Kinlock et al., 2018), and latitudinal gradients in assemblage composition and biological interactions have also been identified (Qian and Ricklefs, 2007; Schemske et al., 2009; Harper and Peck, 2016). Although latitudinal diversity gradients are widely quoted as one of the basic laws in ecology (Lomolino, 2004), studies investigating these patterns on a global scale have largely overlooked intertidal ecosystems, despite their global range and importance as unique environments. An intertidal study in the north eastern Pacific Ocean has reported no latitudinal gradients in species richness (Blanchette et al., 2008). It has been suggested that the latitudinal diversity gradient is biogeographically structured (Roy et al., 1994), and since the above investigations were conducted in a temperate biogeographic region, diversity gradients might not occur. For example, while gastropod richness is relatively stable in the temperate region, a steep rise in richness occurs in the tropical region (Fenberg and Rivadeneira, 2019). However, the few published global-scale intertidal studies have revealed contradictory latitudinal patterns with richness of macroalgae (Konar et al., 2010) and small echinoderms peaking at higher latitudes, while large echinoderms peaking at low latitudes (Iken et al., 2010). Consequently, no latitudinal trend in intertidal diversity was demonstrated across 13 marine regions, as overall assemblage richness was similar everywhere despite strong species-specific patterns (Cruz-Motta et al., 2010). While the studies above benefited from analyzing data collected using a standardized protocol across latitudes, they were limited by only including a smaller number of sampling sites (~70 sites), which may have restricted their capacity to demonstrate large-scale trends, or to separate these from locally entrained patterns. Compared to these studies, we expanded the spatial resolution sixfold and present α-diversity data from 433 rocky shore sites across 155° of latitude, yet there was still no support for the latitudinal α-diversity gradient hypothesis. Latitudinal patterns were, however, markedly different between hemispheres, with a mid-latitude peak in the north and a non-significant unimodal trend in the south. We also demonstrated that different functional groups had different patterns (see later discussion).

Northern hemisphere α-diversity displayed a distinct mid-latitude peak, a pattern also reported for open water biodiversity (Chaudhary et al., 2016). The cause behind this remains elusive but hypotheses include a mid-dominance effect where high- and low-latitude species ranges overlap (Powell et al., 2012), and that high temperatures near the equator, and low temperatures and ice cover at the poles, may limit diversity outside of the temperate zone. However, we did not see this pattern in the temperate zone in the southern hemisphere, and we show that water temperatures do not control intertidal α-diversity. Instead, the open water mid-latitudinal richness peak has recently been suggested to be an artifact of limited sampling in less studied tropical regions (Menegotto and Rangel, 2018). This asymmetric sampling holds true for the northern hemisphere rocky intertidal ecosystem as well. While we show data from 11 sites at latitudes above 70° (south and north combined), 44% (n = 193) of all studies were from mid-latitude rocky shores in Europe and North America between 24 and 60°N. Data from low-latitude regions is furthermore limited because rocky shorelines are scarce in the tropics (Fenberg and Rivadeneira, 2019), and the majority of tropical studies have focused on specific taxonomic groups, and not assemblage-wide patterns (e.g., Hartnoll, 1976; Chelazzi and Vannini, 1980; Mille-Pagaza et al., 2002; Konar et al., 2010; Flores-Rodríguez et al., 2014; Rivadeneira et al., 2015; Gbedemah, 2017; Hendrickx et al., 2019, but see Lourenço et al., 2020). Thus, mid-latitudinal rocky shores, especially in the northern hemisphere, are among the most sampled globally, and future tropical collections may change understanding of biodiversity patterns. This disproportionate sampling pattern is also present in other marine habitats, where sampling intensity is highest at mid-latitudes in the northern hemisphere (Menegotto and Rangel, 2018). Polar data were also limited in the southern hemisphere where species diversity was relatively constant across latitudes with a non-significant dip at the highest latitude. A factor that powerfully affects intertidal diversity in the high polar regions is that most intertidal areas are encased in ice for large parts, if not all of the year, which reduces intertidal biodiversity to zero in these areas. Thus, the Arctic permafrost coastline represents around 34% of the world’s coastline (Lantuit et al., 2012), and Antarctica accounts for around 2.7% (45,317 km). However, in the Antarctic only around 12% of that is ice-free in summer (5468 km) and at both poles much less, if any, is ice-free in winter (Peck, 2018). This lack of ice-free intertidal areas in high polar sites restricts the capacity for communities to establish, and limits the development of macroalgal species (Miller and Pearse, 1991; Mystikou et al., 2014). However, high polar shores are seldom studied (Thyrring et al., 2021), and in this study high-latitude Antarctic shorelines were only represented by a single study at a single site, thus this depression in number of species should be interpreted with caution. In ice-free areas, and at lower Antarctic latitudes, previous biodiversity studies have shown intertidal species richness to increase with latitude from the southern Atlantic ocean to the Antarctic Peninsula (Waller, 2008; Griffiths and Waller, 2016). The situation is clearly complex and, in general, more information on the influence of geographic sampling biases across latitudes are, together with more data from more sites, especially the polar and tropical regions, required to understand and model large-scale biodiversity patterns.

Latitudinal diversity gradients have been identified in subtidal habitats (e.g., Roy et al., 2000; Edgar et al., 2017), and the contrasting patterns in biodiversity gradients between intertidal and subtidal ecosystems may originate from differences in the factors controlling α-diversity. Subtidal benthic habitats experience relatively homogenous environmental conditions across small-to-moderate spatial scales, and large-scale latitudinal diversity gradients are therefore less impacted by local factors and reflect larger-scale changes in environmental conditions. For instance, decreasing water temperatures profoundly affect physiological performance of marine ectotherms (e.g., growth, reproduction and activity, muscle performance) and special adaptations are required to live in cold waters (Peck, 2016; Peck, 2018). Latitudinal biodiversity gradients in ectotherms have, therefore, been attributed to changes in water temperatures in subtidal habitats (Roy et al., 2000; Edgar et al., 2017). Our data run contrary to this and indicate that organisms inhabiting the intertidal do not conform to classic LDG patterns, or any of the correlated environmental drivers (e.g., water temperature). Instead, we demonstrate significant effects of local-scale variation in salinity, wave exposure, ice scour, and biogenic habitats (i.e., macroalgal canopy cover), indicating that α-diversity is determined by small-scale processes through biological interactions and overlapping environmental gradients. Indeed, the importance of microhabitats; biological interactions; and local physical parameters on intertidal community assemblage, density, and diversity has been studied extensively on rocky intertidal shores (Paine, 1974; Archambault and Bourget, 1996; Coleman et al., 2006; Scrosati and Heaven, 2007; Piazzi et al., 2018; Sejr et al., 2021), and high spatial variable in intertidal community coverage and structure between closely adjacent areas are well known (Underwood and Chapman, 1996; Benedetti-Cecchi and Cinelli, 1997; Sejr et al., 2021). For example, algal canopies provide food, predation refugium, and a moist understory environment where species sensitive to desiccation stress can survive during low tides (Dayton, 1975), thus the understory community may be rich and distinct from the surrounding open shoreline (Benedetti-Cecchi et al., 2001; Bulleri et al., 2002; Watt and Scrosati, 2013). Canopies and surface topography also shelter organisms from extreme air temperatures, thereby increasing local coverage and diversity, and facilitating recolonization of the substrate by secondary species (Blanchard and Bourget, 1999; Watt and Scrosati, 2013; Ørberg et al., 2018). Moreover, wave exposure decreases diversity, especially at high latitude where floating ice is pushed on the shore by waves, amplifying the impacts of ice scour on exposed coastlines (Scrosati and Heaven, 2007). Even wave splashing and low water timing (i.e., low tides occurring in the tropics during the hottest time of the day are more harmful than low tides in the middle of the night [Helmuth et al., 2002]) have been shown to be important. All these factors interact and change in non-latitudinally related patterns, creating highly localized unique environmental conditions that, in combination, functionally override latitudinal stress or energy gradients. Thus, although the strength and direction of environmental stressors change across latitudes (i.e., from ice scour and extreme sub-zero temperature in polar regions to desiccation and acute heat stress in the tropics), the combined stress experienced by resident populations may roughly balance out across most latitudes and produce no clear gradient. In a recent study in South America, α-diversity generally did not conform to latitudinal diversity gradient predictions for either the Atlantic or Pacific coast (Cruz‐Motta et al., 2020), and we demonstrate similar patterns on a global scale. A conclusion from this research is that instead of focusing on large-scale drivers that implicitly assumed patterns and mechanisms are scale invariant, focus should be on how scales relevant to the organisms affect distribution patterns, as these factors outweigh latitudinal drivers of biodiversity. Notably, intertidal γ-diversity, which may in addition be controlled by regional drivers (e.g., upwelling, ocean currents) may display different latitudinal patterns than described here (Cruz‐Motta et al., 2020), but global γ-diversity patterns remain to be explored in depth. In polar regions where glaciers or a persistent ice foot covers the intertidal, no communities are present (Dayton, 1990; Peck, 2018), and where the intertidal is characterized by a seasonal ice-foot formation, colonization is only possible during spring and summer (Barnes, 1999), although in these areas it is possible that a few species might persist in highly saline brine pools in winter (Clarke and Beaumont, 2020). Any future reduction of glaciers, ice shelves, and ice foots in these regions will therefore permit organisms to colonize the intertidal zone (Węsławski et al., 2011; Kennicutt et al., 2014; Kennicutt et al., 2015; Kennicutt et al., 2019). This will result in increased biomass and α-diversity as is evident from the high Arctic (Weslawski et al., 2010; Kortsch et al., 2012; Thyrring et al., 2021). Given that most studies fail to consider the aspects discussed above, and that data available from tropical and polar regions are sparse or absent, the relative impacts of multiple stressors ranging from local scale to large scale cannot be explored in depth for coastal biodiversity. However, our data show that latitudinal patterns in intertidal biodiversity differ significantly from subtidal studies. They also highlight the need for re-evaluating conservation efforts and climate change predictions for this unique global ecosystem.

Latitudinal gradients in functional groups

In contrast to overall α-diversity, we identified strong latitudinal patterns in assemblage composition. The relative dominance of the four functional groups changed with latitude, with predators displaying the strongest gradient. In the northern hemisphere, the number of predators decline at latitudes above 60°N, supporting local-scale studies showing that although intertidal predators, such as crabs and starfish, live on mid-latitude intertidal shores (Paine, 1974; Jenkins et al., 2008), they are missing or rare on exposed Arctic shorelines (Thorson, 1934; Weslawski et al., 2010; Høgslund et al., 2014). A latitudinal reduction in the number of predators is consistent with the hypothesis of a general reduced level of predation in the polar regions (Aronson et al., 2007; Schemske et al., 2009; Peck, 2018), and similar latitudinal predation clines are found in subtidal benthic habitats (Taylor and Taylor, 1977; Harper and Peck, 2016). Several explanations have been proposed for this latitudinal decline, including the lack of durophagous predators in Antarctica, which has been attributed to the increased solubility of calcium carbonate (CaCO3) with latitude, increasing the cost of shell production (Vermeij, 1987; Watson et al., 2017). The power muscles generation is also strongly affected by temperature and has been proposed as another reason why durophages are absent in high polar regions (Aronson et al., 2007; Peck, 2018), but there is still no full consensus on the processes behind these patterns, and the mechanisms explaining predation gradients remain debated. Predators are important for organizing food webs and community composition. For example, the habitat-forming mussel Mytilus californianus expanded and excluded more than 25 species after removal of its main predator, the starfish Pisaster, in Western North America (Paine, 1974), and in the absence of predators, a tropical seagrass community can support 10 times more species than with predators present (Freestone et al., 2011). Thus, the low proportion of predators at high latitudes indicates that the importance of top-down biological interactions in assemblages decreases with latitude in marine ecosystems, and that species composition in polar regions may primarily be controlled by the physical environment (Barnes, 2002; Schemske et al., 2009; Høgslund et al., 2014; Peck, 2018). However, it should be noted that we only studied benthic intertidal predators collected during low tides. Therefore, we cannot estimate impacts of mobile predators (e.g., fish and birds) or the predation pressure from epibenthic species during submersion of the shoreline, which can be significant (Bertness et al., 1981; Ellis et al., 2007).

Algae formed the most abundant functional group at high latitudes, consisting of both encrusting species surviving in surface depressions, protected from ice scour and extreme temperatures, and larger macrophytes. Some macrophytes found in polar regions are able to vegetatively regenerate tissue after substantial ice damage to fronds and fastholds, and can survive extreme sub-zero air temperatures during emersion (Kiirikki and Ruuskanen, 1996). Indeed air temperature is the prime stressor in polar intertidal systems (Peck et al., 2006; Thyrring et al., 2020), and microhabitats created by macroalgal fronds are important on intertidal shores across latitudes as they provide protection that produces microclimates permitting survival during temperature extremes. The highest latitude shores contain fewer visually obvious species (Waller, 2008), but coralline algae and encrusting species are capable of surviving extreme temperatures in microhabitats, and can survive the winter underneath the intertidal ice foot, likely in protected air pockets or saline brine tidepools reached occasionally and/or tidally by water (Thyrring et al., 2017; Clarke and Beaumont, 2020).

In this study, grazers occurred across all latitudes. Grazing can have powerful effects on algal biomass and distribution at low- to mid-latitudes (Jenkins et al., 2008; Hawkins et al., 2019), and grazing is an important driver of microalgal assemblage structure on King George Island in Antarctica (Zacher et al., 2007). Contrary to this, no patellid limpet has been reported from south Greenland, suggesting grazing is of little importance in this sub-Arctic region (Høgslund et al., 2014). Canopy-forming macroalgae are key structuring species at high latitudes (Ørberg et al., 2018), yet data on grazing impacts in polar intertidal habitats are currently too limited to make conclusions on the ecological implications. Indeed, more quantitative data are needed to fully understand the ecological implications of temporal and spatial changes in functional diversity, but overall our research demonstrates that small-scale local environmental factors dictate the biodiversity patterns observed and overwhelm any large-scale patterns investigated.

Materials and methods

Data collection

Request a detailed protocol

In 2019 we searched the literature using Web of Science, Google Scholar, and Scopus for publications devoted to rocky intertidal local species richness, community composition, assemblage composition or α-diversity across all latitudes from the Antarctic to Arctic. Reports were considered valid when authors presented assemblage-wide data on ‘α-diversity’, ‘richness’, or ‘number of species’ in terms of species lists, tables, or graphs. Specifically, we excluded studies that focused exclusively on algae (e.g. Konar et al., 2010; Gbedemah, 2017), invertebrates (e.g. Hartnoll, 1976; Mille-Pagaza et al., 2002), or individual taxonomic groups. Full species lists were extracted whenever possible, and WebPlotDigitizer was used to extract graphic data (Drevon et al., 2017). Additional data were obtained from the Multi-Agency Rocky Intertidal Network (pacificrockyintertidal.org), the Alaska Ocean Observing System (aoos.org), and the South American Research Group on Coastal Ecosystems (SARCE) via the Ocean Biogeographic Information System (obis.org) website. The supplementary Appendix 1 – data sources provides a complete literature list and is available from the Zenodo repository (Thyrring and Peck, 2021).

In any of the sources used in this study to abstract data, the fraction of species sampled at a given site (inventory completeness) depends on the sampling effort by the collector and local conditions. There are statistical methods to minimize the influence of sampling biases, such as the jackknife estimator of species richness (Smith and Pontius, 2006), which has been developed to estimate regional γ-diversity from subsamples and abundance (Cruz-Motta et al., 2010 Chao and Chiu, 2016). However, we focused on α-diversity obtained from published records, tables, and figures, and these data were often presented as a single value without supplementary information on subsamples or abundances. Thus, we were unable to verify and compare sampling efforts within sites, but instead we collected data from 433 sites to give very large spatial resolution. To further minimize sampling bias among sites, we only considered studies where all species were collected according to described proscribed set of methodologies.

We validated taxonomic status by checking all the names in the World Register of Marine Species (WoRMS; http://www.marinespecies.org). Whenever species lists were available, we allocated species to functional groups (Figure 1). Species were allocated to one of four functional groups based on food acquisition; algae (all primary producers), grazers (including scrapers), predators (including scavengers), and suspension feeders. Finally, distribution of algae, grazers, predators, and suspension feeders were analyzed and presented as the relative percentage .

Statistical analysis

Request a detailed protocol

The full dataset included latitude, longitude, country, and five latitudinal-scale oceanographic drivers (chlorophyll a, nitrate, phosphate, sea surface temperature, and silicate) available from NASA Earth Observations (NEO; https://neo.sci.gsfc.nasa.gov) and US National Ocean and Atmospheric Administration (NOAA) World Ocean Atlas (https://www.nodc.noaa.gov) as 1° latitudinal × 1° longitudinal cells. For each of the 433 α-diversity sites, we extracted values of chlorophyll a, nitrate, phosphate, sea surface temperature, and silicate from the 1° cell (latitudinal × longitude) in which the site was found. If data was unavailable from the specific cell, we used the value from the nearest available cell within the same latitudinal band (data was unavailable from the specific 1° cell for 5% of the sites).

Relationships between co-variates were assessed using Pearson’s correlation coefficients (Zuur et al., 2010). This revealed collinearity between nitrate and phosphate (Pearson r = 0.94), silicate and temperature (Pearson r = 0.73), and oxygen and temperature (Pearson r = 0.99). Based on the correlation coefficients and the variance inflation factor (VIF; threshold <10) we excluded nitrate, oxygen, and silicate from models to avoid inference from correlation between covariates (Montgomery and Peck, 1992; Zuur et al., 2010) The final models therefore included latitude, chlorophyll a, sea surface temperature, and phosphate. Latitudinal α-diversity patterns in the northern and southern hemispheres were analyzed using R (R Development Core Team, 2019). Initial data exploration revealed a non-linear residual pattern in the northern hemisphere, and we therefore analyzed this data using a GAMM with a negative binomial distribution (θ = 9; Link function = log) from the gamm4 R package (Wood and Scheipl, 2020). The smoother for latitude was fitted with a cubic spline regression (cs), and the optimal amount of smoothing was selected using cross-validation (Zuur et al., 2014). The southern hemisphere α-diversity pattern was analyzed using a GLMM with a negative binomial distribution (θ = 15.58; Link function = log). In both models, negative binomial distributions were fitted to ensure acceptable residual patterns and avoid over-dispersion. We included the factor ‘country’ as a random intercept. The final models were validated by plotting residuals versus fitted values, versus each covariate in the model, and versus each covariate not in the model (Zuur et al., 2014; Zuur and Ieno, 2016). Bubble plots, and variograms from the R package gstat (Gräler et al., 2016), were used to inspect model residuals for spatial correlations, and showed no indication of spatial autocorrelation in either model.

To investigate the general influence of local environmental stressors (stressors relevant to species at site-level), data on four local-scale environmental drivers (ice scour, macroalgal cover, salinity, and wave exposure) were extracted from the published literature (see Materials and method, described above) when available through text, tables, or direct author correspondence. We focused on ‘macroalgal cover’ of canopy-forming algae (i.e., non-canopy algae were excluded from this dataset) because canopies create understory living space and a predation refugia that increase α-diversity (Benedetti-Cecchi et al., 2001; Bulleri et al., 2002; Watt and Scrosati, 2013). This data extraction resulted in four individual datasets, which were analyzed separately as most papers only investigated one of the oceanographic drivers; effects of salinity and wave exposure were separately analyzed using fitted standardized input variables (scale from 0 to 1) to estimate effect sizes of data originating from different scales, depending on the literature source (Schielzeth, 2010). We analyzed salinity effects using a generalized linear model (GLM) with a Poisson distribution, while ice scour and macroalgal cover effects were analyzed using a GLM with a negative binomial distribution because α-diversity was characterized by over-dispersion (determined after visual inspection of residual patterns), but showed no zero-inflation (Hilbe, 2011). Wave exposure effects were initially analyzed using a GLM, but model residuals showed a non-linear pattern, and the data were re-analyzed using a generalized additive model (GAM) with a negative binomial distribution to ensure acceptable residual patterns (Zuur, 2012). All models were finally validated by inspection of standardized residual patterns plotted against fitted values. Effect sizes for GLM models were estimated using bias-corrected parametric bootstrap methods (10,000 iterations) (Zuur and Ieno, 2016), and the effect of the variable was considered significant when 95% CI did not overlap zero.

Data availability

Data used in this study are freely available from the Multi-Agency Rocky Intertidal Network (MARINe; pacificrockyintertidal.org), the Ocean Biogeographic Information System (OBIS; obis.org) and the Alaska Ocean Observing System (AOOS; https://gulf-of-alaska.portal.aoos.org). Publications obtained from the literature search can be found on the Zenodo data repository (Thyrring & Peck, 2021).

The following data sets were generated
    1. Thyrring J
    2. Peck LS
    (2021) Zenodo
    Literature and data for "Global gradients in intertidal species richness and functional groups".
    https://doi.org/10.5281/zenodo.3942061
The following previously published data sets were used

References

    1. Barnes DKA
    (1999) The influence of ice on polar nearshore benthos
    Journal of the Marine Biological Association of the United Kingdom 79:401–407.
    https://doi.org/10.1017/S0025315498000514
    1. Barnes DKA
    (2002) Polarization of competition increases with latitude
    Proceedings of the Royal Society of London. Series B: Biological Sciences 269:2061–2069.
    https://doi.org/10.1098/rspb.2002.2105
    1. Clarke A
    2. Gaston KJ
    (2006) Climate, energy and diversity
    Proceedings of the Royal Society B: Biological Sciences 273:2257–2266.
    https://doi.org/10.1098/rspb.2006.3545
  1. Book
    1. Darwin CR
    (1859)
    The Origin of Species by Means of Natural Selection, or the Preservation of Favoured Races in the Struggle for Life
    John Murray.
  2. Book
    1. Dayton PK
    (1990)
    Polar benthos
    In: Smith W, editors. Polar Oceanography. Academic Press. pp. 631–685.
  3. Book
    1. De Broyer C
    2. Koubbi P
    3. Griffiths HJ
    4. Raymond B
    5. Udekem d’Acoz C
    6. Van de Putte AP
    7. Danis B
    (2014)
    Biogeographic Atlas of the Southern Ocean
    In: De Broyer C, Koubbi P, Griffiths H. J, Raymond B, Udekem d’Acoz C, Van de Putte AP, Danis B, David B, Grant S, Gutt J, Held C, Hosie G, Huettmann F, Post A, Ropert-Coudert Y, editors. Scientific Committee on Antarctic Research. International Science Council. pp. 3–7.
    1. Hendrickx ME
    2. Salgado-Barragán J
    3. Cordero-Ruiz M
    (2019)
    Moluscos litorales (Bivalvia, Gastropoda, Polyplacophora, cephalopoda) de playas rocosas de la región de Guaymas, golfo de california, méxico = littoral mollusks (Bivalvia, Gastropoda, Polyplacophora, cephalopoda) from rocky beaches in the area of Guaymas, gulf of California, Mexico
    Geomare Zoológica 1:51–88.
  4. Book
    1. Lomolino MV
    (2004)
    Conservation biogeography
    In: Lomolino M. V, Heaney L. R, editors. Frontiers of Biogeography: New Directions of the Geography of Nature. Sinauer Associates, Inc. pp. 293–296.
    1. Mille-Pagaza S
    2. Carrillo-Laguna J
    3. Pérez-Chi A
    4. Sánchez-Salazar ME
    (2002)
    Abundancia y diversidad de los invertebrados litorales de isla socorro, archipiélago revillagigedo, méxico
    Revista De Biologia Tropical 50:97–105.
    1. Miller KA
    2. Pearse JS
    (1991)
    Ecological studies of seaweeds in McMurdo Sound, Antarctica
    American Zoologist 31:35–48.
  5. Book
    1. Montgomery DC
    2. Peck EA
    (1992)
    Introduction to Linear Regression Analysis
    Wiley-Interscience.
  6. Software
    1. R Development Core Team
    (2019) R: A language and environment for statistical computing
    R Foundation for Statistical Computing, Vienna, Austria.
  7. Book
    1. Wallace AR
    (1878)
    Tropical Nature and Other Essays
    Macmillan Publishers Limited.
  8. Book
    1. Zuur AF
    (2012)
    A Beginner’s Guide to Generalized Additive Models with R
    Highland Statistics Ltd.
  9. Book
    1. Zuur AF
    2. Saveliev AA
    3. Ieno EN
    (2014)
    A Beginner’s Guide to Generalised Additive Mixed Models with R
    Highland Statistics Ltd.

Decision letter

  1. Detlef Weigel
    Senior and Reviewing Editor; Max Planck Institute for Developmental Biology, Germany
  2. Phillip Fenberg
    Reviewer
  3. Lisandro Benedetti-Cecchi
    Reviewer; Universita di Pisa, Italy

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

This study addresses an important topic of broad ecological interest and provides important insights into the role of local-scale processes in shaping patterns of species diversity, in order to (i) assess if there is a global latitudinal diversity gradient (using α diversity) of rocky shore organisms and its functional groups and, (ii) whether there are any large scale or local environmental predictors of richness patterns. The strength of this paper is the global coverage of studies analyzed, showing for the first time that rocky shore richness does not appear to peak in the tropics – in contrast to many other studies of marine and terrestrial systems.

Decision letter after peer review:

[Editors’ note: the authors submitted for reconsideration following the decision after peer review. What follows is the decision letter after the first round of review.]

Thank you for submitting your work entitled "Global gradients in intertidal species richness and functional diversity" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and a Senior Editor. The following individual involved in the review of your submission has agreed to reveal their identity: Phillip Fenberg (Reviewer #2).

Our decision has been reached after consultation between the reviewers. Based on these discussions and the individual reviews below, we regret to inform you that your work will not be considered further for publication in eLife.

Your manuscript asks important questions related to the global distribution of intertidal marine biodiversity. Both reviewers and editors agreed that the global scope of the manuscript was interesting and that it could be a significant addition to the literature. Nonetheless, one of the external reviewers found important gaps in the statistical treatment of the data, and both external reviewers identified gaps in the literature review that can potentially affect your dataset and the outcomes arising from its analysis. We think however that a new improved version (addressing comments below) could be resubmitted to eLife. Please note that it would be treated as a new submission, but that we would try to recruit the same outside experts for evaluation.

Reviewer #1:

This manuscript tests for consistent latitudinal gradients in species richness, including in the species richness of some specified functional groups. The question is interesting and important, and has not been investigated extensively in rocky shore systems as far as I can tell (but see below). The authors' approach is to extract raw richness values (counts of species) from the literature, and then test for relationships between richness and a range of explanatory variables including latitude, using GLMs, GLMMs, and GAMs. I am not certain, but I think only univariate relationships are fitted to the data. The northern and southern hemispheres are tested separately as well as together, but different coastlines (e.g., E Pacific vs W Pacific vs W Atlantic) appear to be lumped together for analysis. Spatial structure (I think) is accounted for by establishing 5-degree "bins" and then putting random effects on those bins. I don't know if sites in the same latitudinal bin but different parts of the world (e.g., Chile vs Australia) are put in the same bin or different bins, but I think they are the same. Some recommendations I have for improvement of the study are related to the following:

1. One of the paper's stated aims is to study functional diversity. However, functional diversity (i.e., the breadth of functional space occupied by a biota, or the number of distinct functions represented by a biota) is not actually investigated in the paper. Rather, the authors analyze overall species richness, and then analyze species richness separately for a few different functional groups.

2. One previous study, by Cruz-Motta et al. 2010, is cited by the manuscript but its approach and findings are not discussed at all, which is surprising since the underlying questions are very similar. Another of which I'm aware of is Rivadeneira et al. (2015, DOI 10.1111/geb.12328) focus on gastropods only, but is also a very extensive study of the eastern Pacific latitudinal gradient. That paper is not cited in the current manuscript. The paper does not clearly explain how the current study improves upon those earlier studies or provides complementary insights, nor is there any explicit comparative discussion of the findings.

3. Presumably the species lists considered represent a very heterogeneous body of studies that differ substantially in their sampling effort and potentially their focal taxa (i.e., some studies may count all species in some unit of area, and others only species from certain taxa). Was there any consistency applied in study selection, such as: a set of focal taxa was identified by the authors, and a study was only included if it identified species comprehensively from each of those focal taxa? Sampling effort would no doubt have varied substantially as well, but there does not appear to be any standardization for sampling effort, either through the application of a richness estimator such as Chao1 or the jackknife, or rarefaction, or including some measure of sampling effort as a covariate in the analysis. Lastly, this is technically a meta-analysis, but meta-analytical methods (which would account for, say, differences in sampling variance among replicates in estimates of richness) do not appear to be employed. Given a large number of study sites, I doubt this would make a huge difference, relative to the other issues, so it is less of a concern to me than the other points.

4. There is a lot of important information missing from the description of the methods.

a. For example, how the random effects structure was imposed isn't entirely clear, although my best guess is that a random effect was placed on every 5 degrees latitudinal band. This isn't really an ideal way of coping with spatial autocorrelation, given that adjacent latitudinal bands are also likely to have some autocorrelation – and dealing with it is complicated by the fact that one of the explanatory variables is latitude (one way of modeling spatial autocorrelation is to use spatial coordinates as a trend surface, so in a sense, at least a specific form of the latitudinal component of spatial autocorrelation is in fact what is being tested for).

b. A "passion distribution" (I presume the authors mean Poisson?) was used to test the salinity model, but a negative binomial error distribution was used for the other model fits. This is confusing because the error distribution should depend on the behavior of the response variable (richness), not the explanatory variables. Moreover, the authors then talk about checking for normality and homoscedasticity which is confusing since a negative binomial distribution of residuals (presumably what was modeled) is neither normal nor homoscedastic.

c. It appears but was not clearly stated that only univariate relationships were tested. Why was there no model selection conducted investigating the potential for multiple variables to independently or interactively influence richness? Also, what are the statistical relationships among the potential explanatory variables? To what extent could variables that have statistically meaningful relationships with richness be coincidental, and due to the fact that they are spatially correlated with other variables that play a more causal role?

d. There are references in the text to global versus "regional" analyses, but what exactly these regions are isn't clear. Is it just northern and southern hemisphere data analyzed separately?

My suggestion for addressing the methods issues would be a more comprehensive presentation (or at least investigation and detailed description) of the model fit diagnostics, such as spatial correlograms of both residuals and random effects estimates for the different model fits. Structure in either would indicate something important missing. I would encourage the authors to revisit model diagnostics for models with non-normal error structures. Especially for things like the negative binomial, interpreting standard model diagnostic plots is extremely difficult and non-intuitive, and I'm concerned from the text about how much the authors dug into this. If the authors aren't sure where to start here, one option is to simulate data from a fitted model, and compare the distribution of the residuals from the model fitted to the real data with what the residuals look like when data are simulated to conform exactly to model assumptions. R package DHARMa will do that. Obviously, though that's not the only way forward. If there's evidence of lack of fit, the authors might consider some models with multiple explanatory variables. I am also surprised the authors didn't consider a model accounting for the potential effects of different coastlines. Along those lines, one way to account for "region" is to include it as a categorical factor in the original analysis, rather than do a completely different analysis. For example, E Pacific, W Pacific, E Atlantic, W Atlantic, etc. Hemisphere effects could be considered in this way also. If not, I would want to see that there is no such structure in the residuals (e.g., if one color-coded residual or random-effects estimates according to the coastline, for example, they would not be grouped together).

5. My overall take on the Discussion was that the study did not permit very strong conclusions to be drawn – for instance, the authors conclude that there is a mid-latitude peak in richness in mid-latitudes in the northern hemisphere but not the southern, and they attribute this to ice scour in the far north and desiccation and temperature at the equator. However, the authors don't clearly link this causal attribution to their analyses – the temperature effect looks non-significant, and the ice scour effect small in magnitude, which does not seem to support the authors' conclusion. My broad recommendation for the Discussion is for the causal inferences to be linked more concretely to the results of the analyses. A more thorough set of analyses that addresses some of the concerns that I raised in the previous points might make this somewhat easier.

Reviewer #2:

I really enjoyed reading this paper as it is very much within my area of expertise. The largely held belief that the latitudinal diversity gradient is a common pattern held across all major ecosystems is mostly substantiated by the literature. However, as you correctly state in your paper, it is rarely studied at the global scale in rocky shore systems. If it has been, then it is usually done using a coarse scale (e.g. 5 latitudinal bins) that also includes shallow sub-tidal environments, and/or using species richness based on overlapping geographic ranges (but see Rivadeneira et al. 2015 along the entirety of the eastern Pacific…below). Your findings suggest that there is no clear latitudinal gradient using α diversity on a global scale. While this is true based on your existing dataset, I do have one major comment that should help in your revision of the article:

There is a clear sampling gap in tropical latitudes (based on your figure 1). While this is partly because the rocky shore is sparse within tropical latitudes (e.g. Fenberg and Rivadeneira 2019; Ecology Letters), you are missing some regions in your dataset that clearly have papers on the α diversity of tropical rocky shore ecosystems. For example, along the eastern Pacific coast, you have no data points along the mainland Pacific coast of Mexico. Please have a look at Figure 3b in Rivadeneira et al. 2015 Global Ecology and Biogeography. This paper (which I am a co-author on) is based on rocky shore gastropods using literature sources across the whole of the eastern Pacific coast. While there are fewer examples of sites within the tropics compared to the temperate regions, there is a clear increase in species richness within the tropical latitudes (and this is just for gastropods). This runs counter to your argument that there is no gradient in α diversity, so you could include it in your paper as a potential outlier. Related to the above: see the following papers: the Revillagigedo Islands in Mexico, in Mille-Pagaza et al. 2002 ("Abundancia y diversidad de los invertebrados litorales de isla Socorro,Archipiélago Revillagigedo, México") they find 161 species of marine invertebrates (across a few different sites). The paper is written in Spanish, so it may have escaped your search terms. Here are a couple of other papers showing high diversity along the pacific rocky coast of Mexico:

Hendrickx et al. 2019: "Moluscos litorales (Bivalvia, Gastropoda, Polyplacophora, Cephalopoda) de playas

rocosas de la región de Guaymas, golfo de California, México" find 113 species.

Flores-Rodríguez, P., et al. (2014) Mollusks of the Rocky Intertidal Zone at Three Sites in Oaxaca,

Mexico. Open Journal of Marine Science, 4, 326-337.

These are just a few examples of Mexican papers that you have missed (that seem to fit your search criteria).

You also have missed tropical papers along the coast of eastern Africa: Tanzania: Hartnoll Estuarine and Coastal Marine Science (1976)

Somalia: Chelazzi and Vannini (1979) "Zonation of Intertidal Molluscs on Rocky

Shores of Southern Somalia".

And it appears that you have missed papers from New South Wales (your example from Figure 1 appears to be in Queensland): Underwood 1981 "Structure of a rocky intertidal community in New South Wales: patters of vertical distribution and seasonal changes".

Benkendorff and A.R. Davis (2002): "Identifying hotspots of molluscan species richness on rocky intertidal reefs"

These are just a few examples of tropical papers that seem to fit your search criteria but are missing in your analysis. The point I am trying to make here is that while there may indeed be fewer papers in the tropical regions, there are definitely more to be found. While this may not change the overall results of your paper, I think you should adjust your search criteria as you are likely missing quite a few from tropical rocky shores. Some of the papers will be in different languages (especially the Spanish papers from Mexico and Latin America), but they should not be discounted in my opinion because some of them clearly show high levels of α diversity within the tropical latitudes. In conclusion, I really do like this paper and I feel like it will make a valuable contribution – but I feel like the dataset is incomplete and would benefit from a more thorough search for papers within tropical rocky shores.

[Editors’ note: further revisions were suggested prior to acceptance, as described below.]

Thank you for submitting your article "Global gradients in intertidal species richness and functional groups" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by Detlef Weigel as the Senior and Reviewing Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Phillip Fenberg (Reviewer #1); Lisandro Benedetti-Cecchi (Reviewer #2).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

As the editors have judged that your manuscript is of interest, but as described below that additional experiments are required before it is published, we would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). First, because many researchers have temporarily lost access to the labs, we will give authors as much time as they need to submit revised manuscripts. We are also offering, if you choose, to post the manuscript to bioRxiv (if it is not already there) along with this decision letter and a formal designation that the manuscript is "in revision at eLife". Please let us know if you would like to pursue this option. (If your work is more suitable for medRxiv, you will need to post the preprint yourself, as the mechanisms for us to do so are still in development.)

Summary:

This study addresses an important topic of broad ecological interest and provides important insights into the role of local-scale processes in shaping patterns of species diversity, aiming to (i) assess if there is a global latitudinal diversity gradient (using α diversity) of rocky shore organisms and its functional groups and, (ii) whether there are any large scale or local environmental predictors of richness patterns. The strength of this paper is the global coverage of studies analyzed, showing for the first time that rocky shore richness does not appear to peak in the tropics – in contrast to many other studies of marine and terrestrial systems. These outcomes are not specific for rocky intertidal systems, with an increasing number of studies showing that the search for global ecological patterns may be elusive. While sampling in the tropics and the polar regions is poor (acknowledged by the authors), this should be viewed as a call for further research in these regions – not as a weakness of the paper per se. There are also some reservations on how the analysis has been conducted, including the lack of standardization of sampling effort and other details (e.g., size of sampling units) to derive a comparable measure of diversity across sites.

Essential revisions:

The latitudinal gradient of diversity has been studied and confirmed in many aquatic and terrestrial habitats and species across the globe. In the vast majority of cases, richness increases towards the tropics. Using an impressive global dataset of latitudinal diversity gradients in 433 rocky intertidal assemblages of algae and invertebrates from the Arctic to the Antarctic, Thyrring and Peck show that rocky shore ecosystems may not follow this general pattern. The authors show that there is no clear latitudinal gradient for rocky shore organisms using alpha diversity - as posited by prevailing theories - although some functional groups exhibit contrasting patterns. Diversity within functional groups of predators, grazers and filter-feeders decreased towards the poles, whereas the opposite was observed for macroalgae. Correlation with environmental drivers highlighted the importance of local-scale processes in driving spatial patterns of diversity in rocky intertidal assemblages. The paper is well written and the many of the analyses are well done, but there is the concern, which the authors acknowledge, that sampling within tropical latitudes is sparse and needs to be carefully considered when interpreting the results of this paper.

1. The relevant data to standardize species richness may not be available from the primary literature. However, it should be possible to employ relevant standardization methods within the 5{degree sign} latitudinal bands in which the data have been aggregated. An analysis based on standardized data, at least for the more data-rich latitudinal bands, must be added.

2. Employ models that allow assessing unimodality, which is stated but untested. At the bare minimum, a quadratic relationship with latitude should be included in the GLMM. As implemented here, the GLMM employed to relate diversity to latitude can only detect linear trends, but not unimodal patterns and the mid-latitude peak suggested by LOESS for the northern hemisphere. To provide a formal test for unimodality, models with or without a quadratic term could be contrasted using standard model comparison procedures. Alternatively, GAM could be used to evaluate nonlinear effects.

3. Clarify whether p-values are relevant or not. As is, it is confusing. For example, the legend of Table 1 mentions p-values, but these are not reported. Materials and methods indicate that 95% confidence intervals are used to take decisions on null hypotheses, suggesting that p-values are not used in the analysis (lines 436-439). Nevertheless, p-values are reported in Table 2.

4. Provide a rationale for distinguishing between canopy and other algal forms (the distinction is compelling, but it is not explained).

5. We like the conclusion on the importance of local-scale processes. This should be placed in the context of previous studies that have quantified patterns and processes at multiple scales reaching the same conclusion.

6. We could not access the data repository indicated in ref. 91, so we could not assess whether the analysis may have missed potentially relevant papers.

7. Provide the number of studies available for each band in an Appendix.

8. The analysis on macroalgae (e.g., Figure 5) distinguished between canopy and no-canopy algae. This is probably correct, but the rationale for this distinction has not been provided. I think some context is needed, especially to clarify the role of algal canopies in maintaining diversified understory assemblages.

9. The overall conclusion that more studies are needed to assess the magnitude and influence of physical and biotic drivers across multiple scales is important and appropriate. However, many studies have examined processes across scales in rocky intertidal systems (including canopy-removal and limpet-exclusion experiments) and many descriptive studies have quantified variation across multiple spatial scales, emphasizing the importance of small-scale variability in pattern of distribution, abundance and diversity of species on rocky shores. A more thorough discussion of this literature (e.g., Underwood & Chapman, Benedetti-Cecchi, Denny, Coleman, Martins, Fraschetti, etc) would be welcome.

10. In the abstract please say something about more sampling being needed in the tropics. Perhaps line 34.

11. Line 85: What do you mean by "inadequate estimates of intertidal areas"? Inadequate in what way? And in what sense are you talking about area?

12. Line 96: Replace "controlled" by "predicted".

13. Figure 2: where are the R2 values on this plot? How do you get an R2 from a LOESS fit…?

14. Salinity is a local and a regional variable in your analyses, please briefly explain why (Figure 3 and 4).

15. Figure 5: hard to see the box plots, consider making them a different colour or shade.

https://doi.org/10.7554/eLife.64541.sa1

Author response

[Editors’ note: the authors resubmitted a revised version of the paper for consideration. What follows is the authors’ response to the first round of review.]

Reviewer #1:

This manuscript tests for consistent latitudinal gradients in species richness, including in the species richness of some specified functional groups. The question is interesting and important, and has not been investigated extensively in rocky shore systems as far as I can tell (but see below). The authors' approach is to extract raw richness values (counts of species) from the literature, and then test for relationships between richness and a range of explanatory variables including latitude, using GLMs, GLMMs, and GAMs. I am not certain, but I think only univariate relationships are fitted to the data. The northern and southern hemispheres are tested separately as well as together, but different coastlines (e.g., E Pacific vs W Pacific vs W Atlantic) appear to be lumped together for analysis. Spatial structure (I think) is accounted for by establishing 5-degree "bins" and then putting random effects on those bins.

I don't know if sites in the same latitudinal bin but different parts of the world (e.g., Chile vs Australia) are put in the same bin or different bins, but I think they are the same.

Yes, sites across all longitudes within in 5° latitudinal bin was put together to analysis the generality of rocky shores LDGs across all oceans. The method of including sites from all longitudes in this way follows previous investigations on LDG in both terrestrial and marine systems [1–6]. We have highlighted this method in the Materials and methods section:

“Sites across all longitudes was binned within corresponding 5° latitudinal bands” (Manuscript line 403-404).

Some recommendations I have for improvement of the study are related to the following:

1. One of the paper's stated aims is to study functional diversity. However, functional diversity (i.e., the breadth of functional space occupied by a biota, or the number of distinct functions represented by a biota) is not actually investigated in the paper. Rather, the authors analyze overall species richness, and then analyze species richness separately for a few different functional groups.

We agree the term ‘functional diversity’ was inappropriate when we actually analyse patterns in various functional groups. We now refer to “functional groups” throughout the text, and have change the title accordingly:

“Global gradients in intertidal species richness and functional groups”

2. One previous study, by Cruz-Motta et al. 2010, is cited by the manuscript but its approach and findings are not discussed at all, which is surprising since the underlying questions are very similar. Another of which I'm aware of is Rivadeneira et al. (2015, DOI 10.1111/geb.12328) focus on gastropods only, but is also a very extensive study of the eastern Pacific latitudinal gradient. That paper is not cited in the current manuscript. The paper does not clearly explain how the current study improves upon those earlier studies or provides complementary insights, nor is there any explicit comparative discussion of the findings.

We agree, the original discussion did not reflect the relevant literature (such as, Cruz-Motta et al. (2010), Konar et al. (2010), Rivadeneira et al. (2015) [7–9]), or explained how we improved upon these. In the revised version we provide a comprehensive discussion about the current knowledge on global intertidal diversity gradients, and highlight how we expand on this and provide new information:

“Although latitudinal diversity gradients are widely quoted as one of the basic laws in ecology [43], studies investigating these patterns on a global scale have largely overlooked intertidal ecosystems, despite their global range and importance as unique environments. […] We also demonstrated different functional groups had different gradients (see later discussion).” (Manuscript lines 204-223).

We have also rephrased parts of the introduction, so it now more correctly describes the current knowledge on the contradicting intertidal gradients:

“Regional-scale intertidal studies have found richness gradients of gastropods along coastlines in the eastern Pacific Ocean [7,15]. However, no latitudinal diversity gradient of gastropods [16] or macroalgal [8] was found on a global-scale across oceans, and assembly-wide studies have found missing [9,11,12] or inverse [17] gradients.” (Manuscript lines 73-76).

3. Presumably the species lists considered represent a very heterogeneous body of studies that differ substantially in their sampling effort and potentially their focal taxa (i.e., some studies may count all species in some unit of area, and others only species from certain taxa). Was there any consistency applied in study selection, such as: a set of focal taxa was identified by the authors, and a study was only included if it identified species comprehensively from each of those focal taxa? Sampling effort would no doubt have varied substantially as well, but there does not appear to be any standardization for sampling effort, either through the application of a richness estimator such as Chao1 or the jackknife, or rarefaction, or including some measure of sampling effort as a covariate in the analysis.

Sampling efforts and inventory completeness are important drivers for the number of species found in any survey. Studies investigating regional richness, or γ-diversity, should, as the reviewer suggest, apply richness estimators, such as Chao1 [18] or jackknife [19]. These methods rely on subsamples and abundances for estimating regional richness. However, we analyse a-diversity from different sites using data obtained from text, figures or tables, which was often only available in the form of a single value. Thus, any supplement information on abundance, sampling protocol or sampling effort was often missing. We acknowledge this caveat, and in order to minimize data biases we collected data from 433 sites to maximize spatial resolution, and minimise the influence from any odd studies. To further minimise sampling bias among sites, we only considered studies where all species was collected, according to the described methodologies. Thus, we excluded studies that focused exclusively on algae, invertebrates or individual taxonomic groups (see also our response to Reviewer 2 below). We have added this information in the Material and Method section:

“In any of the sources used in this study to abstract data, the fraction of species sampled at a given site (inventory completeness) depends on the sampling effort by the collector and local conditions. […] To further minimize sampling bias among sites, we only considered studies where all species were collected according to described proscribed set of methodologies. (Manuscript lines 380-389).

Lastly, this is technically a meta-analysis, but meta-analytical methods (which would account for, say, differences in sampling variance among replicates in estimates of richness) do not appear to be employed. Given a large number of study sites, I doubt this would make a huge difference, relative to the other issues, so it is less of a concern to me than the other points.

4. There is a lot of important information missing from the description of the methods.

a. For example, how the random effects structure was imposed isn't entirely clear, although my best guess is that a random effect was placed on every 5 degrees latitudinal band. This isn't really an ideal way of coping with spatial autocorrelation, given that adjacent latitudinal bands are also likely to have some autocorrelation – and dealing with it is complicated by the fact that one of the explanatory variables is latitude (one way of modeling spatial autocorrelation is to use spatial coordinates as a trend surface, so in a sense, at least a specific form of the latitudinal component of spatial autocorrelation is in fact what is being tested for).

We agree the statistical analysis needed much revision, and we thank the reviewer for the constructive suggestions. Accordingly, we re-analysed the data and looked for spatial autocorrelation using Variograms and bubble plots. The data actually showed clear autocorrelation, and accordingly all GLMMs were fitted with a Matérn spatial covariance structure with longitudinal and latitudinal values as spatial variables. We still nested sites in 5° latitudinal bins and used these as a random factor to account for dependency among sites nested within bands. The final GLMMs were validated by plotting standardized residual patterns plotted against fitted values, and using residual diagnostic tools form the DHARMa R package [28]. Please see the fully revised Materials and methods section for more details.

b. A "passion distribution" (I presume the authors mean Poisson?) was used to test the salinity model, but a negative binomial error distribution was used for the other model fits. This is confusing because the error distribution should depend on the behavior of the response variable (richness), not the explanatory variables. Moreover, the authors then talk about checking for normality and homoscedasticity which is confusing since a negative binomial distribution of residuals (presumably what was modeled) is neither normal nor homoscedastic.

We apologies for the typo. Indeed, we meant a Poisson distribution, and this has been corrected in the text.

The local models (ice-scour, macroalgal cover, salinity and wave exposure) were analysed as 4 separate models because the available data varied among the four drivers. Consequently, we had four datasets with individual response variables (richness), and the error distribution was therefore different in each analysis to ensure valid models were used. We have clarified this in the Material and Method section to avoid any confusion:

“Because the available data varied among the four parameters, this data extraction resulted in four individual datasets, which were analysed separately as most papers only investigated one of oceanographic drivers” (Manuscript lines 422-424).

c. It appears but was not clearly stated that only univariate relationships were tested. Why was there no model selection conducted investigating the potential for multiple variables to independently or interactively influence richness? Also, what are the statistical relationships among the potential explanatory variables? To what extent could variables that have statistically meaningful relationships with richness be coincidental, and due to the fact that they are spatially correlated with other variables that play a more causal role?

We have clarified how the data exploitation process was conducted before the analysis. We tested for collinearity among environmental covaries, using a variance inflation factor (VIF) analysis, which showed no collinearity. Since environmental drivers may interact (e.g. temperature and salinity), the original full model included interactions among all covariates. However, there were no significant interactions, and we therefore excluded interactions in the final model [29]. Please see the substantially revised Materials and methods section.

d. There are references in the text to global versus "regional" analyses, but what exactly these regions are isn't clear. Is it just northern and southern hemisphere data analyzed separately?

The regional analyses were meant as an analysis focusing on small-scale drivers (e.g. salinity, wave exposure) that changes on a small spatial scale and among sites. We have rephrased the text, and now the “regional” analyses are called “local-scale” to symbolise its drivers relevant to specific sites. This has been changed throughout the text and we have expanded the describing in the Materials and Method section:

To investigate the general influence of local environmental stressors (stressors relevant to species at site-level), data on four local-scale environmental drivers (ice-scour, macroalgal cover, salinity and wave exposure) were extracted from the published literature (see method described above) when available through text, tables or direct author correspondence” (Manuscript lines 419-422).

My suggestion for addressing the methods issues would be a more comprehensive presentation (or at least investigation and detailed description) of the model fit diagnostics, such as spatial correlograms of both residuals and random effects estimates for the different model fits. Structure in either would indicate something important missing. I would encourage the authors to revisit model diagnostics for models with non-normal error structures. Especially for things like the negative binomial, interpreting standard model diagnostic plots is extremely difficult and non-intuitive, and I'm concerned from the text about how much the authors dug into this. If the authors aren't sure where to start here, one option is to simulate data from a fitted model, and compare the distribution of the residuals from the model fitted to the real data with what the residuals look like when data are simulated to conform exactly to model assumptions. R package DHARMa will do that. Obviously, though that's not the only way forward. If there's evidence of lack of fit, the authors might consider some models with multiple explanatory variables. I am also surprised the authors didn't consider a model accounting for the potential effects of different coastlines. Along those lines, one way to account for "region" is to include it as a categorical factor in the original analysis, rather than do a completely different analysis. For example, E Pacific, W Pacific, E Atlantic, W Atlantic, etc. Hemisphere effects could be considered in this way also. If not, I would want to see that there is no such structure in the residuals (e.g., if one color-coded residual or random-effects estimates according to the coastline, for example, they would not be grouped together).

5. My overall take on the Discussion was that the study did not permit very strong conclusions to be drawn – for instance, the authors conclude that there is a mid-latitude peak in richness in mid-latitudes in the northern hemisphere but not the southern, and they attribute this to ice scour in the far north and desiccation and temperature at the equator. However, the authors don't clearly link this causal attribution to their analyses – the temperature effect looks non-significant, and the ice scour effect small in magnitude, which does not seem to support the authors' conclusion. My broad recommendation for the Discussion is for the causal inferences to be linked more concretely to the results of the analyses. A more thorough set of analyses that addresses some of the concerns that I raised in the previous points might make this somewhat easier.

We have, as suggested, rewritten large parts of the discussion to better link our results and the discussion. We hope out revised discussion more clearly link our work and previous studies. Please see the revised Discussion section.

Reviewer #2:

I really enjoyed reading this paper as it is very much within my area of expertise. The largely held belief that the latitudinal diversity gradient is a common pattern held across all major ecosystems is mostly substantiated by the literature. However, as you correctly state in your paper, it is rarely studied at the global scale in rocky shore systems. If it has been, then it is usually done using a coarse scale (e.g. 5 latitudinal bins) that also includes shallow sub-tidal environments, and/or using species richness based on overlapping geographic ranges (but see Rivadeneira et al. 2015 along the entirety of the eastern Pacific…below). Your findings suggest that there is no clear latitudinal gradient using α diversity on a global scale. While this is true based on your existing dataset, I do have one major comment that should help in your revision of the article:

There is a clear sampling gap in tropical latitudes (based on your figure 1). While this is partly because the rocky shore is sparse within tropical latitudes (e.g. Fenberg and Rivadeneira 2019; Ecology Letters), you are missing some regions in your dataset that clearly have papers on the α diversity of tropical rocky shore ecosystems. For example, along the eastern Pacific coast, you have no data points along the mainland Pacific coast of Mexico. Please have a look at Figure 3b in Rivadeneira et al. 2015 Global Ecology and Biogeography. This paper (which I am a co-author on) is based on rocky shore gastropods using literature sources across the whole of the eastern Pacific coast. While there are fewer examples of sites within the tropics compared to the temperate regions, there is a clear increase in species richness within the tropical latitudes (and this is just for gastropods). This runs counter to your argument that there is no gradient in α diversity, so you could include it in your paper as a potential outlier. Related to the above: see the following papers: the Revillagigedo Islands in Mexico, in Mille-Pagaza et al. 2002 ("Abundancia y diversidad de los invertebrados litorales de isla Socorro,Archipiélago Revillagigedo, México") they find 161 species of marine invertebrates (across a few different sites). The paper is written in Spanish, so it may have escaped your search terms. Here are a couple of other papers showing high diversity along the pacific rocky coast of Mexico:

Hendrickx et al. 2019: "Moluscos litorales (Bivalvia, Gastropoda, Polyplacophora, Cephalopoda) de playas

rocosas de la región de Guaymas, golfo de California, México" find 113 species.

Flores-Rodríguez, P., et al. (2014) Mollusks of the Rocky Intertidal Zone at Three Sites in Oaxaca,

Mexico. Open Journal of Marine Science, 4, 326-337.

These are just a few examples of Mexican papers that you have missed (that seem to fit your search criteria).

You also have missed tropical papers along the coast of eastern Africa: Tanzania: Hartnoll Estuarine and Coastal Marine Science (1976)

Somalia: Chelazzi and Vannini (1979) "Zonation of Intertidal Molluscs on Rocky

Shores of Southern Somalia".

And it appears that you have missed papers from New South Wales (your example from Figure 1 appears to be in Queensland): Underwood 1981 "Structure of a rocky intertidal community in New South Wales: patters of vertical distribution and seasonal changes".

Benkendorff and A.R. Davis (2002): "Identifying hotspots of molluscan species richness on rocky intertidal reefs"

These are just a few examples of tropical papers that seem to fit your search criteria but are missing in your analysis. The point I am trying to make here is that while there may indeed be fewer papers in the tropical regions, there are definitely more to be found. While this may not change the overall results of your paper, I think you should adjust your search criteria as you are likely missing quite a few from tropical rocky shores. Some of the papers will be in different languages (especially the Spanish papers from Mexico and Latin America), but they should not be discounted in my opinion because some of them clearly show high levels of α diversity within the tropical latitudes. In conclusion, I really do like this paper and I feel like it will make a valuable contribution – but I feel like the dataset is incomplete and would benefit from a more thorough search for papers within tropical rocky shores.

We agree that there is a rich literature on diversity from intertidal shores in various tropical regions, including the very nice examples listed above. Unfortunately, the majority of tropical studies focus on specific taxonomic groups, or they only investigate either the flora or fauna. This is also the case for the studies suggested above, specifically:

Hendrickx et al. 2019 – invertebrates

Flores-Rodríguez, P., et al. (2014) – molluscs

Rivadeneira et al. (2015) – gastropods

Fenberg and Rivadeneira (2019) – gastropods

Mille-Pagaza et al. (2002) –invertebrates

Benkendorff & Davis (2002) – molluscs

Chelazzi and Vannini (1979) – invertebrates

Hartnoll (1976) –invertebrates

Benkendorff & Davis (2002) – molluscs

We are grateful for the introduction to the paper by Underwood (1981), which data we now include (e.g. see figure 1). Unfortunate, the other papers are unsuitable for our study as we focus on assembly-wide a-diversity, and therefore cannot include papers on specific taxa or groups (i.e. only algae or animals). We clearly had not highlighted the selection criteria in the original method sections, and we have made substantial changes to the Material and Method section to avoid confusion: For example:

“Reports were considered valid when authors presented assembly-wide data on “a-diversity”, “richness” or “number of species” in terms of species lists, tables or graphs. Specifically, we excluded studies that focused exclusively on algae [26,49], invertebrates [50,51] or individual taxonomic groups [e.g. 24,25,51–53,87]. Full species lists were extracted whenever possible, and WebPlotDigitizer was used to extract graphic data [90].” (Manuscript lines 370-375).

References:

1. Malizia A, Blundo C, Carilla J, Osinaga Acosta O, Cuesta F, Duque A, et al. Elevation and latitude drives structure and tree species composition in Andean forests: Results from a largescale plot network. PLoS One. 2020;15: e0231553. doi:10.1371/journal.pone.0231553

2. Menegotto A, Rangel TF. Mapping knowledge gaps in marine diversity reveals a latitudinal gradient of missing species richness. Nat Commun. 2018;9: 4713. doi:10.1038/s41467-01807217-7

3. Chaudhary C, Saeedi H, Costello MJ. Bimodality of latitudinal gradients in marine species richness. Trends Ecol Evol. 2016;31: 670–676. doi:10.1016/j.tree.2016.06.001

4. Mateo RG, Broennimann O, Normand S, Petitpierre B, Araújo MB, Svenning J-C, et al. The mossy north: an inverse latitudinal diversity gradient in European bryophytes. Sci Rep. 2016;6: 25546. doi:10.1038/srep25546

5. Valdovinos C, Navarrete SA, Marquet PA. Mollusk species diversity in the Southeastern Pacific: why are there more species towards the pole? Ecography 2003;26: 139–144.

doi:10.1034/j.1600-0587.2003.03349.x

6. Roy K, Jablonski D, Valentine JW, Rosenberg G. Marine latitudinal diversity gradients: tests of causal hypotheses. Proc Natl Acad Sci U S A. 1998;95: 3699–702.

doi:10.1073/PNAS.95.7.3699

7. Rivadeneira MM, Alballay AH, Villafaña JA, Raimondi PT, Blanchette CA, Fenberg PB. Geographic patterns of diversification and the latitudinal gradient of richness of rocky intertidal gastropods: the ‘into the tropical museum’ hypothesis. Glob Ecol Biogeogr.

2015;24: 1149–1158. doi:10.1111/geb.12328

8. Konar B, Iken K, Cruz-Motta JJ, Benedetti-Cecchi L, Knowlton A, Pohle G, et al. Current patterns of macroalgal diversity and biomass in northern hemisphere rocky shores. Thrush S, editor. PLoS One. 2010;5: e13195. doi:10.1371/journal.pone.0013195

9. Cruz-Motta JJ, Miloslavich P, Palomo G, Iken K, Konar B, Pohle G, et al. Patterns of spatial variation of assemblages associated with intertidal rocky shores: a global perspective. PLoS One. 2010;5: e14354. doi:10.1371/journal.pone.0014354

10. Lomolino M V. Conservation biogeography. In: Lomolino M V., Heaney LR, editors. Frontiers of Biogeography: New Directions of the Geography of Nature. Sinauer Associates, Inc; 2004. pp. 293–296.

11. Fenberg PB, Menge BA, Raimondi PT, Rivadeneira MM. Biogeographic structure of the northeastern Pacific rocky intertidal: the role of upwelling and dispersal to drive patterns. Ecography 2015;38: 83–95. doi:10.1111/ecog.00880

12. Blanchette CA, Miner CM, Raimondi PT, Lohse D, Heady KEK, Broitman BR.

Biogeographical patterns of rocky intertidal communities along the Pacific coast of North

America. J Biogeogr. 2008;35: 1593–1607. doi: 10.1111/J.1365-2699.2008.01913.X

13. Schoch GC, Menge BA, Allison G, Kavanaugh M, Thompson SA, A. Wood S. Fifteen degrees of separation: Latitudinal gradients of rocky intertidal biota along the California Current. Limnol Oceanogr. 2006;51: 2564–2585. doi:10.4319/lo.2006.51.6.2564

14. Iken K, Konar B, Benedetti-Cecchi L, Cruz-Motta JJ, Knowlton A, Pohle G, et al. Largescale spatial distribution patterns of echinoderms in nearshore rocky habitats. PLoS One. 2010;5: e13845. doi:10.1371/journal.pone.0013845

15. Fenberg PB, Rivadeneira MM. On the importance of habitat continuity for delimiting biogeographic regions and shaping richness gradients. Ecol Lett. 2019;22: 664–673.

doi:10.1111/ele.13228

16. Miloslavich P, Cruz-Motta JJ, Klein E, Iken K, Weinberger V, Konar B, et al. Large-scale spatial distribution patterns of gastropod assemblages in rocky shores. PLoS One. 2013;8:

e71396. doi:10.1371/journal.pone.0071396

17. Griffiths HJ, Waller CL. The first comprehensive description of the biodiversity and biogeography of Antarctic and Sub-Antarctic intertidal communities. J Biogeogr. 2016;43: 1143–1155. doi:10.1111/jbi.12708

18. Yang W, Ma K, Kreft H. Geographical sampling bias in a large distributional database and its effects on species richness-environment models. J Biogeogr. 2013;40: 1415–1426.

doi:10.1111/jbi.12108

19. Smith CD, Pontius JS. Jackknife estimator of species richness with S-PLUS. J Stat Softw. 2006;15: 1–12. doi:10.18637/jss.v015.i03

20. Chao A, Chiu C-H. Nonparametric Estimation and Comparison of Species Richness. eLS.

Chichester, UK: John Wiley & Sons, Ltd; 2016. pp. 1–11.

doi:10.1002/9780470015902.a0026329

21. Gbedemah ST. Current patterns in intertidal macro-algal diversity and zonation of two sites on Ghana’s coast. J Oceanogr Mar Res. 2017;5: 159. doi:10.4172/2572-3103.1000159

22. Hartnoll RG. The ecology of some rocky shores in tropical East Africa. Estuar Coast Mar 1976;4: 1–21. doi:10.1016/0302-3524(76)90002-5

23. Mille-Pagaza S, Carrillo-Laguna J, Pérez-Chi A, Sánchez-Salazar ME. Abundancia y diversidad de los invertebrados litorales de isla Socorro, Archipiélago Revillagigedo, México. Rev Biol Trop. 2002;50: 97–105.

24. Flores-Rodríguez P, Flores-Garza R, García-Ibáñez S, Torreblanca-Ramírez C, GaleanaRebolledo L, Santiago-Cortes E. Mollusks of the rocky intertidal zone at three sites in Oaxaca, Mexico. Open J Mar Sci. 2014;04: 326–337. doi:10.4236/ojms.2014.44029

25. Chelazzi G, Vannini M. Zonation of intertidal molluscs on rocky shores of Southern Somalia. Estuar Coast Mar Sci. 1980;10: 569–583. doi:10.1016/S0302-3524(80)80076-4

26. Hendrickx ME, Salgado-Barragán J, Cordero-Ruiz M. Moluscos litorales (Bivalvia, Gastropoda, Polyplacophora, Cephalopoda) de playas rocosas de la región de Guaymas, golfo de California, México = Littoral mollusks (Bivalvia, Gastropoda, Polyplacophora, Cephalopoda) from rocky beaches in the area of Guaymas, Gulf of California, Mexico.

GeomareZoológica. 2019;

27. Benkendorff K, Davis AR. Identifying hotspots of molluscan species richness on rocky intertidal reefs. Biodivers Conserv. 2002;11: 1959–1973. doi:10.1023/A:1020886526259

28. Hartig F. DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models. R package version 0.3.2.0. https://CRAN.R-project.org/package=DHARMa.

29. Ieno EN, Zuur AF. A Beginners’s Guide to Data Exploration and Visualisation with R. Newburgh, United Kingdom: Highland Statistics Ltd.; 2015.

[Editors’ note: what follows is the authors’ response to the second round of review.]

Essential revisions:

1. The relevant data to standardize species richness may not be available from the primary literature. However, it should be possible to employ relevant standardization methods within the 5{degree sign} latitudinal bands in which the data have been aggregated. An analysis based on standardized data, at least for the more data-rich latitudinal bands, must be added.

Following a complete re-analysis of the original data, we no longer aggregate a-diversity within 5° latitudinal bands. The use of 5° latitudinal bands provided a too coarse environmental resolution; local conditions vary substantially between sites in different oceans/regions despite being located within the same 5° of latitude. We therefore extracted environmental data from 1° x 1° latitudinal x longitudinal cells, and allocated the corresponding environmental value to sites located within each cell (please see the much revised method section, and our reply to point 2 below). This approach provides a much more detailed dataset.

Consequently, we do not standardize data across 5° latitudinal bands, and in the discussion, we highlight that patterns of local a-diversity cannot be compared with patterns in regional diversity where other processes (e.g. upwelling, ocean currents) may be important driving factors. We do agree that more studies on global patterns in regional intertidal diversity is needed, but this is beyond the scope of our contribution.

2. Employ models that allow assessing unimodality, which is stated but untested. At the bare minimum, a quadratic relationship with latitude should be included in the GLMM. As implemented here, the GLMM employed to relate diversity to latitude can only detect linear trends, but not unimodal patterns and the mid-latitude peak suggested by LOESS for the northern hemisphere. To provide a formal test for unimodality, models with or without a quadratic term could be contrasted using standard model comparison procedures. Alternatively, GAM could be used to evaluate nonlinear effects.

We thank the reviewers for the good and constructive advice on the statistical protocol.

We acknowledge the issue of a non-linear pattern, and we have re-analyzed the dataset. The residuals distribution in the northern hemisphere dataset did indeed indicate a non-linear pattern, and we re-analysed these data using GAMM. The southern hemisphere dataset is suitable for linear models (e.g. no residuals patterns or other indications), and these data have been analyzed using a GLMM. The whole data exploration process was redone.

Specifically, we started with a full dataset that included latitude, and five latitudinal-scale oceanographic drivers (chlorophyll a, nitrate, phosphate, sea surface temperature and silicate). However, when looking at relationships between co-variates using Pearson’s correlation coefficients, high collinearity was observed between nitrate and phosphate (Pearson r = 0.94), silicate and temperature (Pearson r = 0.73), and oxygen and temperature (Pearson r = 0.99). Based the variance inflation factor (VIF; threshold <10 (Montgomery & Peck, 1992)) we excluded nitrate, silicate and oxygen from the models. In the reduced models was there no collinearity. Thus the final models included chlorophyll a, phosphate and sea surface temperature.

The GAMM and GLMM was fitted with a negative binomial distribution to ensure acceptable residual patterns and to avoid over-dispersion. Both models included the factor ‘country’ as a random effect. For the GAMM, a cubic spline regression smoother was fitted for latitude, and the optimal amount of smoothing was selected using cross-validation (Zuur et al., 2014). The final model was validated by plotting residuals versus fitted values, versus each covariate in the model, and versus each covariate not in the model, following the protocal decribed by (Zuur et al., 2014; Zuur & Ieno, 2016). Bubble plots and variogram were used to inspect model residuals for spatial correlations. They showed no indication of spatial autocorrelation among sites. Final model validation showed no residuals patterns or violations, and we used both models to analyse latitudinal patterns.

Because of these extensive changes, the entire ‘statistical analysis’ section has been rewritten. Please read the revised method section.

3. Clarify whether p-values are relevant or not. As is, it is confusing. For example, the legend of Table 1 mentions p-values, but these are not reported. Materials and methods indicate that 95% confidence intervals are used to take decisions on null hypotheses, suggesting that p-values are not used in the analysis (lines 436-439). Nevertheless, p-values are reported in Table 2.

In the revised manuscript we present p-values for all analyses

4. Provide a rationale for distinguishing between canopy and other algal forms (the distinction is compelling, but it is not explained).

We have explained why we only focus on canopy forming. Please see a detailed response to this issue under point 8 below.

5. We like the conclusion on the importance of local-scale processes. This should be placed in the context of previous studies that have quantified patterns and processes at multiple scales reaching the same conclusion.

We have revised large sections of the discussion to better place our work in context of previous studies on small-scale processes. Please see the revised discussion.

6. We could not access the data repository indicated in ref. 91, so we could not assess whether the analysis may have missed potentially relevant papers.

We will make the data repository public upon acceptance of the manuscript. We have submitted the list of data manuscripts/databases for the reviewers to see if we missed any papers they find relevant. All data will be freely available upon acceptance.

7. Provide the number of studies available for each band in an Appendix.

This has been done, and there is a reference in the Results section:

“Only 11 sites were located at latitudes above 70°, south and north combined (Appendix 1)” (Lines 102-103)

8. The analysis on macroalgae (e.g., Figure 5) distinguished between canopy and no-canopy algae. This is probably correct, but the rationale for this distinction has not been provided. I think some context is needed, especially to clarify the role of algal canopies in maintaining diversified understory assemblages.

We have clarified the important role of intertidal algal canopies in the introduction:

“Macroalgal canopies shelter understory species from environmental stress (Krause-Jensen et al., 2016; Sejr et al., 2021), creating protective microhabitats that increase environmental heterogeneity and biodiversity, hereby maintaining a diversified understory assemblages (Bulleri et al., 2002; Watt & Scrosati, 2013; Piazzi et al., 2018).” (lines 54-58)

In the result section:

“We focused exclusively on canopy-forming algae because canopies provide living space and protection from predation and extreme temperatures, thereby increasing a-diversity and coverage of understory organisms” (lines 118-120)

In the method section:

“We focused on ‘macroalgal cover’ of canopy-forming algae (i.e. non-canopy algae were excluded from this dataset) because canopies creates understory living space and a predation refugium that increase a-diversity (Benedetti-Cecchi et al., 2001; Bulleri et al., 2002; Watt & Scrosati, 2013).” (lines 468-470)

9. The overall conclusion that more studies are needed to assess the magnitude and influence of physical and biotic drivers across multiple scales is important and appropriate. However, many studies have examined processes across scales in rocky intertidal systems (including canopy-removal and limpet-exclusion experiments) and many descriptive studies have quantified variation across multiple spatial scales, emphasizing the importance of small-scale variability in pattern of distribution, abundance and diversity of species on rocky shores. A more thorough discussion of this literature (e.g., Underwood & Chapman, Benedetti-Cecchi, Denny, Coleman, Martins, Fraschetti, etc) would be welcome.

We have added a more in-depth discussion of the importance of small-scale variability on rocky shores. Please see the revised discussion.

10. In the abstract please say something about more sampling being needed in the tropics. Perhaps line 34.

We have added the following sentence addressing sampling in the tropics:

“Polar and tropical intertidal data were sparse, and more sampling in tropical regions is required to improve knowledge of marine biodiversity. (lines 36-37).

11. Line 85: What do you mean by "inadequate estimates of intertidal areas"? Inadequate in what way? And in what sense are you talking about area?

What we meant to say was, that the geographic area (in km2) of intertidal zones hasn’t been estimated on a global scale. In 2019, Murray et al. mapped for the first time the extent of tidal flats, but the extent of other intertidal habitats remain unknow.

We have clarified this in the text:

“Global assessments of intertidal biodiversity have been hindered by data scarcity from polar intertidal shores, and a lack of estimates of intertidal geographic areas. In fact only the extent of intertidal mud flats have been quantified on a large scale (Murray et al., 2019).” (lines 84-86)

12. Line 96: Replace "controlled" by "predicted"

Done.

13. Figure 2: where are the R2 values on this plot? How do you get an R2 from a LOESS fit…?

This was an unfortunate typo. There are no R2 values presented on Figure 2. We have updated the figure legend:

“Figure 2: Latitudinal patterns in rocky intertidal a-diversity plotted against latitude. Data are split into southern and northern hemisphere. A linear regression line (southern hemisphere) and a best-fit locally weighted scatterplot smoother (northern hemisphere) was added with 95% confidence intervals to aid visual interpretation” (lines 958-961)

14. Salinity is a local and a regional variable in your analyses, please briefly explain why (Figure 3 and 4).

We have given this much thought, and we have decided to exclude salinity as a latitudinal factor. We originally included this parameter as salinity usually reaches maximum values around latitudes 20°N and 20°S, and decreases toward high latitudes. However, the latitudinal data is obtained from satellite data, and coastal salinity is too strongly affected by coastal processes (such as precipitation and melt water run off; Valiela et al., 2012; Meire et al., 2017; Sejr et al., 2017; Mouginot et al., 2019; Covernton & Harley, 2020; Duarte et al., 2020; Navarro et al., 2020). Thus, satellite derived salinity data may provide limited ecological relevant data, and we wanted to avoid making non-logical conclusions on the impacts of salinity, which is a key abiotic driver for ecological patterns (as we show on the local-scale analysis)

15. Figure 5: hard to see the box plots, consider making them a different colour or shade.

We agree. The box plots have been highlighted and superimposed on the data points (see revised Figure 3).

References:

Benedetti-Cecchi, L., Pannacciulli, F.G., Bulleri, F., Moschella, P.S., Airoldi, L., Relini, G. & Cinelli, F. (2001) Predicting the consequences of anthropogenic disturbance: large-scale effects of loss of canopy algae on rocky shores. Marine Ecology Progress Series, 214, 137– 150.

Bulleri, F., Benedetti-Cecchi, L., Acunto, S., Cinelli, F. & Hawkins, S.J. (2002) The influence of canopy algae on vertical patterns of distribution of low-shore assemblages on rocky coasts in the northwest Mediterranean. Journal of Experimental Marine Biology and Ecology, 267, 89– 106.

Covernton, G. & Harley, C. (2020) Multi-scale variation in salinity: a driver of population size and structure in the muricid gastropod Nucella lamellosa. Marine Ecology Progress Series, 643, 1– 19.

Duarte, C.M., Rodriguez-Navarro, A.B., Delgado-Huertas, A. & Krause-Jensen, D. (2020) Dense Mytilus Beds Along Freshwater-Influenced Greenland Shores: Resistance to Corrosive Waters Under High Food Supply. Estuaries and Coasts, 43, 387–395.

Krause-Jensen, D., Marbà, N., Sanz-Martin, M., Hendriks, I., Thyrring, J., Carstensen, J., Sejr, M.

& Duarte, C. (2016) Long photoperiods sustain high pH in arctic kelp forests. Science Advances, 2, e1501938.

Meire, L., Mortensen, J., Meire, P., Juul-Pedersen, T., Sejr, M.K., Rysgaard, S., Nygaard, R., Huybrechts, P. & Meysman, F.J.R. (2017) Marine-terminating glaciers sustain high productivity in Greenland fjords. Global Change Biology, 23, 5344–5357.

Montgomery, D.C. & Peck, E.A. (1992) Introduction to Linear Regression Analysis, WileyInterscience.

Mouginot, J., Rignot, E., Bjørk, A.A., van den Broeke, M., Millan, R., Morlighem, M., Noël, B., Scheuchl, B. & Wood, M. (2019) Forty-six years of Greenland Ice Sheet mass balance from 1972 to 2018. Proceedings of the National Academy of Sciences of the United States of America, 116, 9239–9244.

Murray, N.J., Phinn, S.R., DeWitt, M., Ferrari, R., Johnston, R., Lyons, M.B., Clinton, N., Thau, D. & Fuller, R.A. (2019) The global distribution and trajectory of tidal flats. Nature, 565, 222– 225.

Navarro, J.M., Détrée, C., Morley, S.A., Cárdenas, L., Ortiz, A., Vargas-Chacoff, L., Paschke, K., Gallardo, P., Guillemin, M.L. & Gonzalez-Wevar, C. (2020) Evaluating the effects of ocean warming and freshening on the physiological energetics and transcriptomic response of the Antarctic limpet Nacella concinna. Science of the Total Environment, 748, 142448.

Piazzi, L., Bonaviri, C., Castelli, A., Ceccherelli, G., Costa, G., Curini-Galletti, M., Langeneck, J., Manconi, R., Montefalcone, M., Pipitone, C., Rosso, A. & Pinna, S. (2018) Biodiversity in canopy-forming algae: Structure and spatial variability of the Mediterranean Cystoseira assemblages. Estuarine, Coastal and Shelf Science, 207, 132–141.

Sejr, M.K., Mouritsen, K.N., Krause-Jensen, D., Olesen, B., Blicher, M.E. & Thyrring, J. (2021) Small scale factors modify impacts of temperature, ice scour and waves and drive rocky intertidal community structure in a Greenland fjord. Frontiers in Marine Science, 7, 607135. Sejr, M.K., Stedmon, C.A., Bendtsen, J., Abermann, J., Juul-Pedersen, T., Mortensen, J. & Rysgaard, S. (2017) Evidence of local and regional freshening of Northeast Greenland coastal waters. Scientific Reports, 7, 13183.

Valiela, I., Camilli, L., Stone, T., Giblin, A., Crusius, J., Fox, S., Barth-Jensen, C., Monteiro, R.O., Tucker, J., Martinetto, P. & Harris, C. (2012) Increased rainfall remarkably freshens estuarine and coastal waters on the Pacific coast of Panama: Magnitude and likely effects on upwelling and nutrient supply. Global and Planetary Change, 92–93, 130–137.

Watt, C.A. & Scrosati, R.A. (2013) Regional consistency of intertidal elevation as a mediator of seaweed canopy effects on benthic species richness, diversity, and composition. Marine Ecology Progress Series, 491, 91–99.

Zuur, A.F. & Ieno, E.N. (2016) A protocol for conducting and presenting results of regression-type analyses. Methods in Ecology and Evolution, 7, 636–645.

Zuur, A.F., Saveliev, A.A. & Ieno, E.N. (2014) A beginner’s Guide to Generalised Additive Mixed models with R, Highland Statistics Ltd., Newburgh, United Kingdom.

https://doi.org/10.7554/eLife.64541.sa2

Article and author information

Author details

  1. Jakob Thyrring

    1. British Antarctic Survey, Cambridge, United Kingdom
    2. Department of Zoology, University of British Columbia, Vancouver, Canada
    3. Arctic Research Centre, Department of Bioscience, Aarhus University, Silkeborg, Denmark
    4. Homerton College, University of Cambridge, Cambridge, United Kingdom
    5. Marine Ecology, Department of Bioscience, Aarhus University, Silkeborg, Denmark
    Contribution
    Conceptualization, Data curation, Formal analysis, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    jakyrr57@bas.ac.uk
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1029-3105
  2. Lloyd S Peck

    British Antarctic Survey, Cambridge, United Kingdom
    Contribution
    Conceptualization, Supervision, Investigation, Visualization, Writing - review and editing
    Competing interests
    No competing interests declared

Funding

Danmarks Frie Forskningsfond (7027-00060B)

  • Jakob Thyrring

Natural Environment Research Council

  • Lloyd S Peck

H2020 Marie Skłodowska-Curie Actions (797387)

  • Jakob Thyrring

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

The authors thank Brenda Konar, Takashi Noda, Tan Koh Siang, Robin Elahi, Valdivia Nelson, and Megan Dethier for providing intertidal data. We gratefully acknowledge the Valdivia Nelson’s projects FONDECYT (grant #1141037) and CONICYT-PIA (ANILLO; grant #ART1101) for providing data from Chile and Antarctic, respectively. This study further utilized data collected by the South American Research Group on Coastal Ecosystems (SARCE) sponsored by TOTAL foundation and data collected by the Multi-Agency Rocky Intertidal Network (MARINe): a long-term ecological consortium funded and supported by many groups. Please visit pacificrockyintertidal.org for a complete list of the MARINe partners responsible for monitoring and funding these data. Data management has been primarily supported by BOEM (Bureau of Ocean Energy Management), NPS (National Parks Service), The David and Lucile Packard Foundation, and United States Navy. This work is a contribution to The future of Arctic biodiversity in a climate change era’ project.

Senior and Reviewing Editor

  1. Detlef Weigel, Max Planck Institute for Developmental Biology, Germany

Reviewers

  1. Phillip Fenberg
  2. Lisandro Benedetti-Cecchi, Universita di Pisa, Italy

Publication history

  1. Received: November 2, 2020
  2. Accepted: March 18, 2021
  3. Accepted Manuscript published: March 19, 2021 (version 1)
  4. Version of Record published: April 8, 2021 (version 2)

Copyright

© 2021, Thyrring and Peck

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 400
    Page views
  • 67
    Downloads
  • 1
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Ecology
    Lei Guo et al.
    Research Article Updated

    The Varroa destructor mite is a devastating parasite of Apis mellifera honeybees. They can cause colonies to collapse by spreading viruses and feeding on the fat reserves of adults and larvae. Amitraz is used to control mites due to its low toxicity to bees; however, the mechanism of bee resistance to amitraz remains unknown. In this study, we found that amitraz and its major metabolite potently activated all four mite octopamine receptors. Behavioral assays using Drosophila null mutants of octopamine receptors identified one receptor subtype Octβ2R as the sole target of amitraz in vivo. We found that thermogenetic activation of octβ2R-expressing neurons mimics amitraz poisoning symptoms in target pests. We next confirmed that the mite Octβ2R was more sensitive to amitraz and its metabolite than the bee Octβ2R in pharmacological assays and transgenic flies. Furthermore, replacement of three bee-specific residues with the counterparts in the mite receptor increased amitraz sensitivity of the bee Octβ2R, indicating that the relative insensitivity of their receptor is the major mechanism for honeybees to resist amitraz. The present findings have important implications for resistance management and the design of safer insecticides that selectively target pests while maintaining low toxicity to non-target pollinators.

    1. Ecology
    2. Evolutionary Biology
    Hanna M Bensch et al.
    Research Article

    Living with relatives can be highly beneficial, enhancing reproduction and survival. High relatedness can, however, increase susceptibility to pathogens. Here, we examine whether the benefits of living with relatives offset the harm caused by pathogens, and if this depends on whether species typically live with kin. Using comparative meta-analysis of plants, animals, and a bacterium (nspecies = 56), we show that high within-group relatedness increases mortality when pathogens are present. In contrast, mortality decreased with relatedness when pathogens were rare, particularly in species that live with kin. Furthermore, across groups variation in mortality was lower when relatedness was high, but abundances of pathogens were more variable. The effects of within-group relatedness were only evident when pathogens were experimentally manipulated, suggesting that the harm caused by pathogens is masked by the benefits of living with relatives in nature. These results highlight the importance of kin selection for understanding disease spread in natural populations.