Integrating laboratory experiments and biogeographic modelling approaches to understand sensitivity to ocean warming in rare and common marine annelids

Among ectotherms, rare species are expected to have a narrower thermal niche breadth and reduced acclimation capacity and thus be more vulnerable to global warming than their common relatives. To assess these hypotheses, we experimentally quantified the thermal sensitivity of seven common, uncommon, and rare species of temperate marine annelids of the genus Ophryotrocha to assess their vulnerability to ocean warming. We measured the upper and lower limits of physiological thermal tolerance, survival, and reproductive performance of each species along a temperature gradient (18, 24, and 30 °C). We then combined this information to produce curves of each species’ fundamental thermal niche by including trait plasticity. Each thermal curve was then expressed as a habitat suitability index (HSI) and projected for the Mediterranean Sea and temperate Atlantic Ocean under a present day (1970–2000), mid- (2050–2059) and late- (2090–2099) 21st Century scenario for two climate change scenarios (RCP2.6 and RCP8.5). Rare and uncommon species showed a reduced upper thermal tolerance compared to common species, and the niche breadth and acclimation capacity were comparable among groups. The simulations predicted an overall increase in the HSI for all species and identified potential hotspots of HSI decline for uncommon and rare species along the warm boundaries of their potential distribution, though they failed to project the higher sensitivity of these species into a greater vulnerability to ocean warming. In the discussion, we provide some caveats on the implications of our results for conservation efforts.


Introduction
The recognition that organisms are unevenly distributed on Earth and that rarity is more frequent than commonness has laid the foundation for the empirical investigation of the causes of rarity and commonness and their consequences on large-scale ecological processes, species biogeography, and biodiversity conservation (Kunin and Gaston 1997). Rarity has been an important driver of species disappearance in past extinctions (Harnik et al. 2012;Calosi et al. 2019), and it is a widely used criterion in conservation biology and ecosystem management to assess the current levels of vulnerability and extinction risk of species (IUCN 2001). Narrow geographic ranges, restricted habitat specificity, small population size, and a combination of these factors make rare species more at risk of decline or extinction compared with their common relatives (Rabinowitz 1981;Calosi et al. 2008a). The common-rare species paradigm not only defines the structure, dynamics, and function of ecosystems but is Communicated by Donovan P German. also useful for projecting the current and future sensitivity of taxa and the impending reorganisation of biodiversity patterns in response to ongoing climate changes (Violle et al. 2017).
Among the attributes of rarity, geographic range size is considered a good predictor of species susceptibility to environmental perturbations (Harnik et al. 2012;Staude et al. 2020). The strict correspondence between the extent of the spatial occurrence of a species and the ecological requirements under which it can thrive gave origin to the niche breadth-range size hypothesis (Brown 1984), which postulates that species with broader ecological niches can persist across a wider range of environmental conditions compared to species with narrower ecological niches. This hypothesis bears important implications for our understanding of the mechanisms underpinning commonness and rarity (Calosi et al. 2010;Slatyer et al. 2013) and provides a framework for the prediction of the redistribution of biodiversity under climate change scenarios (Sunday et al. 2011(Sunday et al. , 2012. Such implications are particularly relevant for ectotherms, which constitute the vast majority of species on Earth (Longhurst 2010), and for which temperature is a major abiotic driver of their past, present, and future distribution via its pervasive effect on basic physiological mechanisms, such as thermal tolerance and energy metabolism (Gagné et al. 2020;Bennet et al. 2021).
In recent decades, considerable effort has been devoted to finding a mechanistic link between the thermal niche breadth and range size of marine species to predict their distribution and abundance changes under ocean warming. The inclusion of species-specific physiological constraints in nichebased models has helped improve our ability to generate more accurate predictions of species range shifts and species losses in the face of ocean warming (e.g., Kearney and Porter 2009;Cheung et al. 2015;Reygondeau 2019). However, quantitative information on the geographic occurrence and distribution of species is often insufficient for most taxa and is disproportionally deficient for marine species. For many marine invertebrates, for example, assessments of rarity based solely on the geographic range size is challenged by patchy distributions, overdispersion, and the colonisation of microrefugia that may be difficult to detect (Chapman 1999;Benkendorff and Przeslawski 2008). Accordingly, whether the applicability and predictive power of the concepts of commonness and rarity also hold true in these cases has yet to be tested on a wider array of taxa whose geographic distribution is not comprehensively characterised. Furthermore, most niche-model studies estimate the fundamental thermal niche of a species based on its mean basal thermotolerance alone, and thus these models do not account for acclimation through phenotypic plasticity that can broaden the thermal requirements necessary for a species to thrive (Sunday et al. 2011). Further improvements to the predictive power of niche-based models can thus be obtained by integrating, within the common-rare paradigm, the assumptions of those climate-based hypotheses that link a species' thermal physiology and geographic range size to its acclimation ability to predict its sensitivity to ocean warming. Among these hypotheses, the climatic variability hypothesis (CVH) posits that the breadth of the thermal tolerance window of a species correlates positively with its latitudinal extent (Stevens 1989). An extension of the CVH asserts that widespread species or species experiencing higher levels of climatic variability (e.g., common species) should have a broader physiological thermal tolerance and greater acclimation capacities (i.e., Brattstrom rule) when compared to less extensively distributed species or species colonising more stable environments (e.g., rare species) (Calosi et al. 2008a,b;Bozinovic et al. 2011;Magozzi et al. 2015).
In this study, we tested the possibility of using the concept of commonness and rarity to assess the relative levels of thermal sensitivity of marine ectotherms with differing biogeographies. To do so, we used seven species of temperate, interstitial annelids of the genus Ophryotrocha Claparède and Mecznikow, 1869 (Dorvilleidae) characterised by different geographic range sizes and abundances as case studies (Prevedelli et al. 2005;Simonini et al. 2009): O. labronica La Greca and Bacci, 1962, O. japonica Paxton and Åkesson, 2010, O. adherens Paavo et al., 2000, O. puerilis Claparède and Mecznikov, 1869, O. diadema Åkesson, 1976, O. hartmanni Huth, 1933, and O. robusta Paxton and Åkesson, 2010 We specifically assessed whether the common or rare status of these species was a good predictor of their vulnerability to ocean warming. To achieve our goal, we first experimentally measured the physiological thermal tolerance limits, survival, and reproductive performance of the studied species along a temperature gradient, i.e., 18, 24, and 30 °C. We then used the generated experimental data to build the fundamental thermal niche of each species following the thermal habitat suitability index (HSI) concept developed by Helaouët and Beaugrand (2009), taking into account the thermal plasticity of the measured traits. Finally, we spatially projected the HSI values onto a present-day  scenario for ocean temperatures and predicted their changes under mid-(2050-2059) and late-(2090-2099) 21st Century scenarios of ocean warming. Two representative concentration pathways were chosen, RCP2.6 and RCP8.5, representing a 'strong mitigation' and 'business-as-usual' climate change scenario, respectively (IPCC 2013). Furthermore, to identify the putative ecological mechanisms underpinning the differences between common and rare species, we explored the validity of the niche breadth-range size hypothesis in the studied species and assessed whether their physiological acclimation ability was related to their level of commonness or rarity as measured by the extent of their geographic range size (Brattstrom rule). Tests were performed on laboratory strains of each species originating from subtidal populations collected in the Mediterranean region, a climate change hotspot, where the average sea surface temperature is predicted to increase by 0.93-1.24 °C over the 2010-2039 period and by 3.45-4.37 °C over the 2010-2099 period under the RCP8.5 scenario (IPCC 2013). Ocean warming is causing the 'tropicalization' of the southeast Mediterranean basin due to the migration of alien tropical species from the Red Sea and the 'meridionalization' of the northern parts of the basin due to the increased occurrence of indigenous thermophilic species, jeopardizing the persistence of rare, endemic, and colder-tolerant species (Coll et al. 2010).

Biology and thermal habitats of the studied species
The seven studied species have a similar ecology, morphology, reproductive biology, and dispersal ability. They are all interstitial generalist grazers found in organic-rich benthic habitats and have an average adult body length that ranges between 3 and 12 mm (O. adherens and O. puerilis, respectively). They reproduce semicontinuously, laying clutches of eggs inside a protective matrix or tube that are cared for by the parents until the eggs hatch, and the embryos undergo direct development (Simonini et al. 2009(Simonini et al. , 2010. Occurrence data for these species have been obtained mainly from the Northern Hemisphere. The coasts of the Shetland Islands (Scotland, UK, 60° 20′ 49.05′′ N, 1° 14′ 8.38′′ W; O. hartmanni) and Oahu (Hawaii,USA,21° 19′ 0.23′′ N,157° 57′ 13.04′′ W;O. labronica) are the most northerly and southerly latitude locations, respectively, where at least one of these species was found (Simonini et al. 2009). All species have been recorded in the subtidal zone (first 60-90 cm below the lower tide level) of shallow marine and brackish coastal waters (Simonini et al. 2010), where they showed different patterns of distribution, prevalence, and abundance across thermal regions and habitats (Prevedelli et al. 2005;Simonini et al. 2009). For example, the distribution range of O. labronica, the most ubiquitous and abundant of the seven species, extends from the temperate coasts of the North Atlantic Sea to the subtropical coasts of the Red Sea (Simonini et al. 2009) and comprises locations, where minimum and maximum annual water temperatures can range between 5 and 30 °C, respectively (Massamba-N'Siala et al. 2011). In contrast, the known occurrence of the rare species O. robusta and O. diadema is limited to Mediterranean-like climates, specifically to locations, where water temperatures do not exceed 27.5 °C.

Assessment of species commonness-rarity status
We identified the level of commonness or rarity of the studied species using two of the three criteria suggested by Rabinowitz (1981): geographic range size and local population abundance. We measured geographic range size as the extent of occurrence (EOO) (Gaston and Fuller 2009), which was calculated with the tool GeoCAT (Bachman et al. 2011) using information on the known global distribution of the studied species in the Northern Hemisphere (see Online Resource 1, Fig. S1; Online Resource 2, Table S1). We calculated local population abundance as the mean maximum number of individuals recorded in each sampled locality (Online Resource 2, Tables S1-S2). We then coded species as common (c), uncommon (u), or rare (r) following a modified version of the 'proportion of species' method described by Gaston (1994). Species with an EOO and local abundance higher than the third quartile of the EOO and abundance distribution were considered common; O. labronica and O japonica. Species with an EOO or local abundance between the first and third quartiles of the EOO and abundance distribution were categorised as uncommon: O. adherens and O puerilis. Species with an EOO or local abundance below the first quartile of the EOO and abundance distribution were considered rare: O. hartmanni, O diadema, and O. robusta (Online Resource 1, Fig. S2).

Specimen maintenance and experimental design
The specimens used in our study came from laboratory cultures maintained for approximately ten generations (approx. For each species, we produced sixty reproductive pairs (parental generation) by pairing sexually mature individuals that were haphazardly collected from the laboratory cultures and by isolating them in a 20 mL glass bowl under the same conditions as the cultures of origin. When the first egg mass was laid and the eggs hatched (F1 generation), the parents were removed and the development of the F1 individuals was followed for approximately 20 d in O. labronica, O. japonica, O. diadema, and O. robusta, 10 d in O. adherens, 35 d in O. puerilis, and 25 d in O. hartmanni. These times corresponded to the late juvenile stage that precedes the full development of oocytes in the coelomatic cavity of the individuals and were used as a temporal reference to randomly assign broods to one of the three acclimation temperatures: 18, 24, and 30 °C, for a total of twenty broods per acclimation temperature, thus ensuring that exposure to the acclimation temperatures started at the same developmental stage under all treatments in all species. Two of the acclimation temperatures, 18 and 24 °C, were chosen within the thermal range normally experienced by all the tested species in their natural environment, while 30 °C was the highest temperature recorded in the site colonised by the common species O. labronica (Massamba-N'Siala et al. 2011). Exposure to 18 and 30 °C was achieved through the increase or decrease in temperature, respectively, from the culture temperature of 24 °C at a rate of 1 °C h −1 (Massamba-N'Siala et al. 2012) using a programmable incubator (Everlasting, AM SLIM 701 TNV, Towcester, UK). At the end of this gradual pre-exposure, for each acclimation temperature, a number of F1 individuals were randomly taken from the broods and assigned to one of two experiments aimed at determining the physiological thermal tolerance limits (Experiment 1) and measuring the survival and reproductive performance (Experiment 2) of the studied species. Each of these experiments is described in detail in the following sections (see also Online Resource 1, Fig. S3 for a schematic representation of the experimental design). During exposure to the acclimation temperatures, the experimental replicates were isolated in 20 mL glass bowls placed inside plastic containers and kept in programmable incubators (Everlasting, AM SLIM 701 TNV, Towcester, UK) with a light-dark regime set at 12:12 h. Salinity (35 ‰) was kept constant throughout the experiments. Artificial sea water was made by dissolving artificial sea salt (Reef Crystals, Instant Ocean, St. Blacksburg, VA, USA) in distilled water. Sea water in the bowls was changed every 3 days, and individuals were fed on the same day ad libitum with a solution of spinach minced in sea water (Massamba-N'Siala et al. 2012).

Experiment 1: Determination of physiological thermal tolerance limits
Upper and lower physiological thermal tolerance limits were measured using as metric maximum and minimum critical temperature limits (CT max and CT min ), respectively (Lutterschmidt and Hutchison 1997). CT max and CT min were defined as the temperature at which we detected specific, sequential behavioural endpoints used as proxies to define the gradual physiological impairment caused by exposure to suboptimal and lethal temperatures (Massamba-N'Siala et al. 2012. Endpoints were measured using a dynamic method (Lutterschmidt and Hutchison 1997), specifically by continuously observing the behaviour of individuals exposed to a constantly increasing or decreasing temperature for CT max and CT min, respectively, per unit of time (1 °C min −1 ), starting at the acclimation temperature to which individuals had been exposed for 7 days (Massamba-N' Siala et al. 2012). Tests were performed on a subset of twenty individuals for CT max and twenty individuals for CT min randomly taken from the broods after an exposure of 7 days to the assigned acclimation temperature. We tested a maximum of five individuals per round inside a borosilicate glass bowl (diam. = 3 cm, depth = 2 cm) immersed in a computer-controlled recirculating ethylene glycol bath (VWR-I462-7017, Radnor, PA, USA). During each round, we continuously measured the temperature of the seawater inside the glass bowl with a digital thermometer (Testo 922, 2-channel Thermocouple Thermometer, Milano, Italy) equipped with a precision fine wire thermocouple (T/C type K, Testo). The size of each individual was measured before starting each round by counting the number of chaetigers, i.e., the segments bearing bristles (Massamba-N'Siala et al. 2012). The temperatures at which the loss of locomotor control, onset of spasms, and death (lethal temperature) occurred were measured on the same individual as (sequential) endpoints to define CT max . In contrast, the temperatures at which the onset of spasms and chill coma occurred were the (sequential) endpoints measured to define CT min (Massamba-N'Siala et al. 2012. More specifically, loss of locomotor control was identified as the occurrence of reversible arrhythmicity in locomotor activities. Lethal temperature was defined as the temperature at which individuals no longer responded to poking and showed no recovery after being cooled back to the corresponding acclimation temperature. The onset of spasms was defined as when spasmodic, uncontrolled contraction of the entire body started. Finally, a chill coma was defined as a reversible status, where motionless individuals did not respond to prodding. We quantified the physiological acclimation ability of the studied species following the approach used by Calosi et al. (2008a). More specifically, we calculated the absolute difference between the mean temperature at which an endpoint (a proxy for CT max or CT min ) was measured after acclimation at 24 and 18 °C (mean CT max/min [24 °C] -mean CT max/min [18 °C]).

Experiment 2: Determination of survival and reproductive performance
Individual survival and reproductive performance were measured for each species in twenty pairs per acclimation temperature during 90 days of exposure to the assigned acclimation temperature. Pairs were formed by haphazardly selecting and matching sexually matured individuals taken from different broods within the same acclimation temperature to avoid inbreeding (Massamba-N'Siala et al. 2011). We checked each pair daily to record the presence or death of each individual and the number of eggs laid, which were 1 3 counted and discarded. We calculated survival as the percentage of days that an individual was alive over the total duration of the experiment (90 days) and reproductive performance as the total fecundity, which was defined as the total number of eggs per chaetiger laid by a pair (Massamba-N'Siala et al. 2011).

Determination of a species' fundamental thermal niche
Following the theoretical approach suggested by Helaouët and Beaugrand (2009; see also Sunday et al. 2012;Pörtner et al. 2007), we used the experimental data on CT max , CT min , and total fecundity to determine the shape of the fundamental thermal niche of the studied species. Curves were built on four thermal ranges-for critical tolerance, tolerance, growth, and reproduction-defined by the mean upper and lower limits for physiological thermal tolerance and reproduction measured at each acclimation temperature and quantified in the temperature dimension as an HSI. This index ranges from 0 (when the environment is unsuitable for the species) to 1 (when the environment is optimal for the species). First, we bounded the thermal niche of each species using the maximum mean value of the upper lethal temperature and the minimum mean temperature at which chill comas occurred from each acclimation temperature. The HSI at these two designated points was set at 0.01, while an HSI below or beyond this range was set as 0. Following the same approach, we set the critical tolerance range (HSI = 0.1) at the mean value of upper lethal temperature and chill coma averaged across acclimation temperatures. We then retrieved the lower and upper limits of the tolerance range (HSI = 0.3) using the minimum mean temperature at the onset of spasms (CT min ) and the maximum mean temperature at the loss of locomotor control, respectively, measured across acclimation temperatures. We set the range for growth (HSI = 0.5) between the maximum mean temperature at the onset of spasms (CT min ) and the minimum mean temperature at the loss of locomotor control. We obtained the limits for the reproduction range (HSI = 0.8) from the curve of reproductive performance derived from the data on total fecundity. We assigned an HSI value of 1 to the temperature at which the maximum number of eggs was estimated to be produced. Finally, we applied a shape-preserving, piecewise cubic interpolation using all defined HSI values to quantify the full extent of the thermal niche from − 4 to 42.7 °C (mean minimum and maximum global annual sea surface temperatures) every 0.1 °C. The interpolated value at a query point was based on a shape-preserving, piecewise cubic interpolation of the values at 35 neighbouring grid points using the function interp1 from MATLAB's statistical toolbox. If a threshold point had not been quantified experimentally (e.g., total fecundity for O. puerilis (u) and O. hartmanni (r)), we extrapolated the HSIs by interpolation using the other existing points.
We calculated the breadth of the fundamental thermal niche of the species as the difference between the maximum and minimum temperatures that defined each of the four ranges used to build the curves of the fundamental thermal niche.

Thermal habitat suitability maps
We extracted the yearly average sea surface and bottom temperature from the IPSL-CM5-MR (Institut Pierre Simon de Laplace), GFDL-ESM2 M (Geophysical Fluid Dynamics Laboratory), and MPI-ESM (Max Plank Institute) Earth system models (ESMs) from the Coupled Model Intercomparison Project Phase 5 (CMIP5) website for two emission scenarios, RCP2.6 and RCP8.5, and for the reference , mid-(2050-2059) and late-(2090-2099) 21st Century periods. We computed average values for each geographic cell on regular 0.5° × 0.5° grids between the ESMs for each year and each decade. Then, using the fundamental niche previously computed, we retrieved the corresponding HSI values averaged for each grid cell using sea surface and bottom temperatures. In this way, we calculated HSI values for each species in the North Atlantic and Mediterranean regions at the different periods for the two RCPs. Since the studied species are associated with coastal environments, HSI values located in the open ocean domain were set to 0 following the biogeographic raster of Reygondeau et al. (2013). In addition, to better visualise the changes in HSI under the two RCPs, we plotted maps to show the mathematical difference between the HSI values projected in the future scenarios and those calculated for the reference period.
Following the approach described above, we also calculated HSI values for each species based on real temperature data for the North Atlantic and Mediterranean regions using the World Ocean Atlas 2018 database (Locarnini et al. 2019) and for the Mediterranean Sea using the MEDAR/ MEDATLAS database (Fichaut et al. 2003). These maps showed a good match with those obtained for the reference period  used in our projections (Online Resource 3, Appendix 1; Online Resource 1, Figs. S4-S5; Online Resource 2, Table S3).

Statistical analyses
We tested the effects of species status (Status; three levels: c, u, r), species identity (Species(Status); seven levels, nested in Status), acclimation temperature (T; two levels: 18, 24 °C), and their interactions on CT max and CT min with generalised linear models (GLMs; family Gaussian, identity link), using body size as a covariate. This analysis was performed for only two of the three acclimation temperatures, because none of the rare species and only one uncommon species, O. adherens, survived a 7 days exposure to 30 °C. We then performed a second analysis with GLMs (family Gaussian, identity link), using body size as a covariate, to test for the effect of species identity (Species; three levels), acclimation temperature (T; three levels: 18, 24, 30 °C), and their interaction on the CT max and CT min of the three species that survived the exposure to 30 °C: O. labronica (c), O. japonica (c), and O. adherens (u). For this latter analysis, there was no replication for the species Status.
Due to technical problems, all the specimens of O. puerilis (u) and O. hartmanni (r) died before being paired for the measurement of survival and reproductive performance. Consequently, we used GLMs to assess the effect of species identity (Species; five levels), acclimation temperature (T; three levels: 18, 24, 30 °C), and their interactions on survival (family quasi-Poisson, log link) and total fecundity (family Gaussian, identity link) for the remaining species, as there was no replication for species Status.
Post hoc tests were performed with Tukey's test using the R package multcomp. Since multiple endpoints, used to define physiological traits (CT max or CT min ), or multiple life-history traits (survival and total fecundity) were measured on the same individual, we applied the Benjamini-Hochberg correction for false positives to all P values.
Preliminary analyses showed that none of the measured traits exhibited a phylogenetic signal (λ ≠ 0, P < 0.05), which was estimated using Pagel's λ using the 'phylosig' function in the R package phytool (for details, see Online Resource 3, Appendix 2; Online Resource 1, Fig.  S6; Online Resource 2, Table S4). Accordingly, we performed final analyses without applying any phylogenetic correction.
We investigated the relationship between the physiological acclimation ability and geographic range size for each proxy of CT max and CT min using a phylogenetic generalised least squares multiple regression with the R package ape, using body size as a covariate. Using the same technique, we assessed the relationship between niche breadth and geographic range size by fitting linear models for the breadth of the fundamental thermal niche as a function of the EOO for each of the four performance ranges previously identified: critical tolerance, tolerance, growth, and reproduction. Assumptions of normality and homogeneity of variance were checked on the residuals with the Shapiro-Wilk test and Levene test, respectively, and were met for all the regression models.
We performed all statistical analyses using R software, version 4.0.0 (RStudio Team 2020; see Online Resource 4 for R scripts).

Physiological upper thermal tolerance limits
Mean CT max (± SE) ranged between 33.19 ± 0.03 and 41.82 ± 0.02 °C in the group of common species (O. labronica and O. japonica), between 27.28 ± 0.79 and 35.70 ± 0.05 °C in the group of uncommon species (O. adherens and O. puerilis), and between 29.33 ± 0.04 and 36.58 ± 0.04 °C in the group of rare species (O. hartmanni, O diadema, and O. robusta) (Fig. 1a-c). Common species had a significantly higher mean CT max than their uncommon and rare congeners both at 18 and 24 °C and for all endpoints (Status by T (18,24 °C) : maximum deviance = 182.48, P < 0.001; Table 1a-c, Fig. 1a-c. See Online Resource 2, Table S5 for pairwise comparisons). Uncommon species showed a slightly but significantly lower mean CT max compared to the rare species when measured as the loss of locomotor control (Fig. 1a) and lethal temperature (Fig. 1c) following exposure to 18 °C. In all other cases, uncommon and rare species showed comparable mean CT max (Fig. 1a-c).
We found significant interspecific differences for all the endpoints used to define CT max (Species(Status) by T (18,24 °C) : maximum deviance = 144.67, P < 0.001; Table 1a-c; Fig. 2a-c). Mean CT max generally increased with increasing acclimation temperatures for all species when measured as the loss of locomotor control (Fig. 2a) and for all species except O. japonica (c) when measured as the onset of spasms (Fig. 2b). This trend was less consistent in the measurement of lethal temperature, for which three out of the five rare-uncommon species showed comparable levels of mean CT max at 18 and 24 °C (Fig. 2c). Mean CT max was the highest in O. labronica (c), followed by O. japonica (c) (Fig. 2a-c. See Online Resource 2, Table S6 for pairwise comparisons and Table S7 for mean raw values). In general, O. diadema (r) and O. robusta (r) showed lower mean CT max than O. japonica (c). Ophryotrocha hartmanni (r) and O. puerilis (u) stood out for having the lowest mean CT max of all species (Fig. 2a-c). Finally, O. adherens (u) had among the lowest mean CT max when measured as the loss of locomotor control (Fig. 2a) and had mean CT max between those measured for O. japonica (c) and O. robusta (r) when measured as the onset of spasms and the lethal temperature (Fig. 2a-c).
The lowest and highest mean (± SE) CT max ranged between 34.34 ± 0.14 °C (O. adherens (u)) and 42.66 ± 0.13 °C (O. labronica (c)), respectively, following exposure to 30 °C (Online Resource 1, Fig. S7; Online Resource 2,  Table S8 for pairwise comparisons) and lethal temperature (Species: deviance = 68.63, P < 0.001; Table 2c; Online Resource 1, Fig.  S8). Finally, mean CT max also increased with increasing acclimation temperature when 30 °C was included in the analysis. The only exception to this pattern was found in the loss of locomotor control in O. japonica (c), for which no significant differences were detected between CT max measured at 24 and 30 °C (Online Resource 1, Fig. S7).

Physiological lower thermal tolerance limits
Mean (± SE) CT min ranged between − 1.69 ± 0.02 and 8.08 ± 0.02 °C in the group of common species, − 2.07 ± 0.05 and 8.53 ± 0.07 °C in the group of uncommon species, and − 2.40 ± 0.01 and 8.20 ± 0.02 °C in the group of rare species (Fig. 1d, e). The common species had slightly lower mean CT min compared to uncommon species when measured as the onset of spasms and had slightly lower mean CT min compared to both uncommon and rare species when measured as chill coma following exposure to 24 °C. No significant difference in mean CT min was found between these three groups in any other cases (Status by T (18,24 Fig. 2d,e). Mean CT min generally increased with increasing acclimation temperatures in all species. An exception to this trend was observed in O. labronica (c), which showed comparable mean CT min at 18 and 24 °C when measured as the chill coma (Fig. 2e). In general, O. adherens (u) and O. hartmanni (r) showed the highest mean CT min , while O. puerilis (u) showed the lowest mean CT min (Fig. 2d,e; Online Resource 2, Table S6 for  pairwise comparisons and Table S7 for mean raw values). The remaining species showed intermediate mean CT min that ranked differently across acclimation temperatures depending on the endpoint considered (Fig. 2d,e).
The lowest and highest mean (± SE) CT min ranged between − 0.63 ± 0.23 °C (O. japonica (c)) and 12.45 ± 0.23 °C (O. adherens (u)), respectively, following exposure to 30 °C (Online Resource 1, Fig. S7; Online Resource 2, Table S7 for mean raw values). Ophryotrocha adherens (u) also had the highest mean CT min following exposure to 30 °C (Species by T (18,24,30 °C) : maximum deviance = 72.94, P < 0.001; Table 2d,e. See Online Resource 2, Table S8 for pairwise comparisons). Finally, mean CT min generally decreased with increasing acclimation temperatures in all three species and for all endpoints when exposure Fig. 1 Difference in a-c upper and d-e lower physiological thermal tolerance limits (CT max and CT min , respectively) among common (c, green border), uncommon (u, orange border) and rare (r, dark-red border) species of the genus Ophryotrocha following 7 d of exposure to 18 and 24 °C (white and grey fill, respectively). The boundary of the box closest to zero indicates the 25th percentile, a horizontal line and a black diamond symbol inside the box mark the median and mean, respectively, and the boundary of the box farthest from zero indicates the 75th percentile. Whiskers above and below the box indicate the 10th and 90th percentiles. Points above and below the whiskers indicate outliers. Capital letters indicate differences between species according to their common, uncommon, and rare status within acclimation temperature. An asterisk indicates differences between acclimation temperatures within status. Statistical significance: P < 0.05. Sample size: 14 ≤ N ≤ 20 1 3 to 30 °C was included in the analyses (Online Resource 1, Fig. S7).

Survival and reproductive performance
At 18 °C, uncommon and rare species showed a comparable mean survival to O. japonica (c) (100% at all acclimation temperatures) and a higher mean (± SE) survival compared to O. labronica (c) (76.78 ± 6.41%) (Species by T (18,24,30 °C) : deviance = 509.30, P < 0.001; Fig. 3a, Table 3a. See Online Resource 2, Table S9 for mean raw values and Table S10 for pairwise comparisons). At 24 °C, the mean survival of O. diadema (r) and O. robusta (r) was comparable to that measured in O. japonica (c), as well as comparable to that measured at 18 °C, while the mean survival of O. adherens (u) was lower than that measured at 18 °C (49.79 ± 2.28%). Survival of O. labronica (c) at 24 °C showed mean values comparable to those measured at 18 °C and were intermediate among those measured for O. adherens (u) and all the other species. At 30 °C, O. japonica (c) and O. labronica (c) (97.89 ± 1.77%) had comparable mean survival, and both common species had a significantly higher mean survival than uncommon and rare species; mean values for survival in the latter two groups ranged between 14.67 ± 2.57 and 11.32 ± 1.53% in O. robusta (r) and O. adherens (u), respectively (Fig. 3a).
Ophryotrocha robusta (r) was the most fecund species at 18 and 24 °C (60.39 ± 3.93 and 67.57 ± 10.16 eggs chaetigers −1 , respectively), followed by O. adherens (u) and O. labronica (c) at 18 °C and by O. adherens (u), O. Table 1 Results of the analysis of deviance for the effect of species Status (common, uncommon, rare), species identity (Species(Status)) and acclimation temperature (T: 18 and 24 °C), and their interactions on the upper (a-c) and lower (d-e) physiological thermal tolerance limits (CT max and CT min , respectively) of the species of the genus Ophryotrocha Degrees of freedom (Df), deviance residuals (Dev. Res.), degrees of freedom residuals (Df Res.), deviance (Dev.), and P-value with Benjamini-Hochberg correction for false positives (P(BH)) are provided , and both common species were on average more fecund than the uncommon and rare species. More specifically, O. diadema (r) and O. robusta (r) did not reproduce, while O. adherens (u) had an extremely low mean total fecundity (0.56 ± 0.20 eggs chaetigers −1 ) (Fig. 3b).

Acclimation ability, thermal niche breadth, and range size
Physiological acclimation ability was not related to the geographic range size of the studied species for all the analysed endpoints (Online Resource2, Table S11). Similarly, we found no significant relationship between the breadth of the fundamental thermal niche and range size in all the performance ranges considered (Online Resource 2, Table S12).

Fundamental thermal niches and thermal habitat suitability
The fundamental thermal niche of the studied species was represented by a unimodal, left-skewed curve, where the thermal habitat suitability (HSI) increases with increasing temperatures until it reaches a peak at the optimal temperature for reproduction and steeply descends thereafter (Fig. 4). Common species extended their HSI towards higher temperatures for all functions and maximised their reproductive performance under warmer conditions (25-27 °C) compared with uncommon and rare species (21-23 °C), thereby intensifying the asymmetry of their curves (Fig. 4). In contrast, the left part of the curves unfolded within a narrower range of temperatures and had no recognisable patterns in variation between species. Projections based on the fundamental thermal niche of the studied species showed that HSI values higher than 0.8 ranged between the 30th and 45th parallel latitudes when mapped in the reference period (Online Resource 1, Figs. S9-S10). An overall increase in the HSI was observed for all species under both RCP scenarios and was much more marked poleward, especially for O. labronica (c) and O. robusta (r) (Fig. 5 and Online Resource 1, Fig. S11). In both the reference and future scenarios, uncommon and rare species tended to possess higher HSI values (0.9-1) in many more cells than their common relatives, particularly at the lowest latitude of their potential distribution (

Discussion
Our study sheds light on a number of fundamental physiological and ecological mechanisms explaining the commonness and rarity of temperate marine ectotherms, showing that uncommon and rare species are more physiologically sensitive to elevated temperatures than common species due to their lower physiological heat tolerance. Our simulations also suggest that uncommon and rare species may be more prone to declines along the warm boundaries of their potential distribution. However, the simulations fail to project these species' higher thermal sensitivity into their overall greater vulnerability to ocean warming. Our results are consistent with the known biogeography of the studied species and support the use of the dichotomous commonness and rarity concept to assess macrophysiological patterns of thermal sensitivity in ectotherm species. In addition, our results indicate that the predictive power of this conceptual framework is less accurate when dealing with intermediate levels of rarity, for example, the uncommon versus rare species in this case.

Upper thermal tolerance: a driving factor in defining biogeographic differences
We find strong support for the idea that common species possess a greater physiological heat tolerance than their uncommon and rare congeners at all acclimation temperatures and for all lethal and sublethal endpoints measured. In addition, common species outperformed uncommon and rare species at the highest acclimation temperature, showing higher survival and reproductive capacities. Similar results were obtained by Thibault et al. (2020), who found that O. japonica (c) was able to persist for two generations at an elevated temperature (28 °C) thanks to the beneficial adjustment of its energy metabolism, as characterised using a targeted metabolomics approach, while O. robusta (r) faced a lethal increase in energy requirements before the first generation of viable offspring could be produced. In contrast, no clear physiological or ecological patterns linked to the biogeographical status of the studied species were found when cold tolerance was considered or when survival and reproductive performance were measured and compared within the 18 to 24 °C range. Our observations provide a mechanistic explanation of the known global, regional, and local geographic distribution of the studied species. For example, the greatest Fig. 3 Difference in a survival and b total fecundity between species of the genus Ophryotrocha measured during 90 d of exposure to 18, 24 and 30 °C (white, grey, and black fill, respectively). Common (c), uncommon (u), and rare (r) species are indicated with green, orange, and dark-red borders, respectively. See Fig. 1 for the description of the box plots. Capital letters indicate differences among species within acclimation temperature. Lowercase letters indicate differences among acclimation temperatures within the same species.   Fig. 4 Relationship between temperature and the habitat suitability index (HSI) defining the fundamental thermal niche of the common (green lines), uncommon (orange lines), and rare (dark red lines) species of the genus Ophryotrocha. Curves were interpolated across different levels of physiological impairment (CT min and CT max ) and reproductive optima. CC = chill coma, LT = lethal temperature, OS = onset of spasms (CT min ), LLC = loss of locomotor control, T optimal = optimal temperature for reproduction. Subscripts min, mean and max refer to the minimum, mean or maximum values among the means calculated at each acclimation temperature for the traits considered heat tolerance of O. labronica (c) accounts for its broadest geographic range size and its ability to persist in tidal and subtidal waters, where temperatures as high as 30 °C have been recorded (Massamba-N'Siala et al. 2011). Not surprisingly, O. labronica (c) has been the only species to date found in subtropical waters, specifically along the coasts of the Northern Red Sea (Simonini et al. 2009). Ophryotrocha japonica (c), a species originally from the Fig. 5 Maps of the difference in the thermal habitat suitability (HSI) between the mid-(2050-2059) and late-(2090-2099) 21st Century scenario, and the present day  scenario under the RCP8.5 scenario in the North Atlantic and Mediterranean regions, obtained for the common (c), uncommon (u), and rare (r) species of the genus Ophryotrocha. Black arrows indicate areas characterized by HSI decreases North Pacific coasts, is the other most heat tolerant species. Since it was first recorded in the Mediterranean Sea, this species was able to colonise the same thermal regions as O. labronica (c) and now has the second largest range size extent among the studied species (Simonini et al. 2009). In contrast, O. hartmanni (r), the least heat tolerant species, is characterised by very low densities in subtidal waters (1-7 indiv.; Online Resources 2, Table S1) but is among the most abundant macrobenthic species in the sediments of fish farms on the North Atlantic coasts (Pereira et al. 2004), suggesting its preference for colder waters. Similarly, O. puerilis (u) is among the most physiologically sensitive species to elevated temperatures and it capitalizes on thermal habitats that better match its preference for colder temperatures. For example, by colonising organic-rich sediments at higher depths (Taboada et al. 2017) or by reaching higher local densities during the coldest months of the year (Prevedelli et al. 2005).
Of course, several factors other than temperature can drive the spatial and temporal distribution of marine species (Gagné et al. 2020). Different sensitivities to other abiotic factors (such as salinity, oxygen, pH), biotic interactions, capacity and limitation to dispersal, habitat preference, and colonisation history can prevent a species from fully realising its fundamental thermal niche, contributing to a mismatch between the physiological thermal tolerance and distribution range of a species (i.e., realised niche, Pulliam 2000;Arribas et al. 2012;Sánchez-Fernández et al. 2012). For example, O. diadema (r), a species initially reported in the Pacific coasts, and O. robusta (r), a species endemic to the Mediterranean Sea (Simonini et al. 2009), display a greater physiological heat tolerance among rare species, but their known distribution is limited to only two and eight known localities, respectively (Online Resource 1, Fig.  S1). For these species, upper thermal tolerance limits are expected to be major determinants of their colonisation success in a new place. However, their recent colonisation history (for O. diadema) or the limited chance of dispersion of the species of the genus Ophryotrocha, which usually passively occurs through the ballast water of ships (Simonini et al. 2009), coupled with the absence of thermal microrefugia along the dispersal trajectories due to the patchy occurrence of suitable habitats, may play an important role in determining the rarity level of these annelid species.

Macrophysiological and macroecological patterns for common and rare species
Over the past decades, an increasing number of studies have attempted to use thermal tolerance limits to search for macrophysiological and macroecological patterns that could inform species vulnerability to climate change (e.g., Calosi et al. 2008aCalosi et al. , 2010Bozinovic et al. 2011;Sunday et al. 2011Sunday et al. , 2012Bennett et al. 2021). Most documented patterns have produced hypotheses linking the organismal thermal physiology and geographic range size. Among these hypotheses, the climatic variability hypothesis (CVH) and its corollaries (e.g., Brattstrom rule) predict a positive association between a species' geographic range size, physiological thermal tolerance, and acclimation capacities (Bozinovic et al. 2011;Magozzi and Calosi 2015). In diving beetles of the genus Deronectes, for example, heat and cold thermal tolerance were positively correlated with range size, and restricted species were foreseen to be most at risk from global warming due to their reduced acclimation capacities (Calosi et al. 2008a(Calosi et al. , 2010. Our results provide partial support for the CVH: a broader range size was associated with a greater physiological heat tolerance in the two common species, but we report no quantitative link between the species range extent and heat tolerance within the group of uncommon and rare species. In addition, no relationship appears to exist between a species' range size and its physiological acclimation ability, and upper and lower thermal tolerance limits increase and decrease with increasing acclimation temperatures, respectively, in all species. The fact that we used strains originating from subtidal temperate populations that are adapted to seasonal thermal changes may explain this result. The existence of a relationship between a species' thermal tolerance breadth and its range size is another mechanism suggested to explain differences in the distribution of common and rare species (Brown 1984). The niche breadth-range size hypothesis has recently received empirical validation (Slayter et al. 2013;Stuart-Smith et al. 2017), and numerous large-scale studies have demonstrated the tight relationship existing between the width of the physiological thermal tolerance window and the extent of geographic range boundaries across latitudinal gradients (e.g., Calosi et al. 2010;Compton et al. 2007;Sunday et al. 2011). In our study, the breadth of the fundamental thermal niche failed to predict the extent of the range size. Our result contributes to the ongoing debate around whether the niche breadth-range size hypothesis is universally valid (e.g., Hirst et al. 2017). A mismatch between the scale at which niche breadth and range size are measured could explain why our data provided no support for this hypothesis (Kambach et al. 2019).
The comparable or higher thermal sensitivity that uncommon species show compared to rare species may (at least in part) help explain why our results provide weak or no support for the tested hypotheses (CVH, niche breadth-range size hypothesis). On one hand, this highlights the limitations of the approach used to categorise species into discrete biogeographical groups within a continuum of forms, commonness and rarity, which are relative concepts. Rare-common species cut-off points can, in fact, change depending on the level of rarity of other congeneric species, the spatial and temporal scale of investigation, or following the inclusion of other criteria (Gaston 1994;Flather and Sieg 2007), potentially leading to different conclusions on the relative differences in thermal sensitivity between uncommon and rare species. Currently, there are only a few continuous indices of rarity, and the most used indices (in the terrestrial realm) are based on exhaustive distributions (mostly from expert range maps) that are not available for marine invertebrates. In fact, for this group of metazoans, gathering information on the criteria used for categorising common and rare species (e.g., range size) is particularly challenging. The use of more flexible categorisation approaches could help increase the accuracy with which closely related species are arranged into intermediate levels of species rarity. This could be achieved by including metrics that are generally considered less rigorous in predicting species distribution, such species prevalence, i.e., the presence/absence ratio to the total number of sites collected (Jiménez-Valverde et al. 2009). For example, the inclusion of information on prevalence among the criteria used to define the common-rare status of our studied species would have made the distinction less marked between O. adherens (u), O. diadema (r), and O. robusta (r), all of which had a prevalence lower than or close to the first quartile (0.08) that defined the uncommon-rare species cut-off point for this trait (Online Resource 2, Table S2). These observations were confirmed by the occurrence data obtained for the studied species in the Mediterranean Sea, a region subjected to a more extensive and systematic decadal collection effort and thus less biased by unbalanced sampling compared to global distribution records (Simonini et al. 2009). On the other hand, we cannot discard the possibility that species thermal physiology alone cannot provide a comprehensive explanation of intermediate levels of rarity. As discussed above, temperature is not the only factor influencing the geographic distribution of uncommon and rare species. The role of other abiotic factors (e.g., Sunday et al. 2014), interspecific interactions (e.g., Godsoe and Harmon 2012), and the number of dispersal opportunities (e.g., Grantham et al. 2003;Arribas et al. 2012) is worth considering in future investigations of the mechanisms explaining rarity.

Habitat suitability shift in common and rare species under ocean warming: conclusions and caveats
The unimodal, left-skewed curves of the fundamental thermal niche of the studied species resembled the thermal performance curves used to describe the thermal sensitivity of biological performances in ectotherms (Schulte et al. 2011). The level of asymmetry of these curves, as well as the position of the thermal optimum relative to the mean of the environmental temperature, have important implications for predicting species vulnerability to ocean warming. In fact, species exposed to thermal conditions that are close to or above their thermal optimum for performance are expected to have narrower or no safety margins to cope with increasing temperatures and, consequently, be more at risk of population declines under ocean warming (Deutsch et al. 2008). In our study, this prediction should have been relevant for O. puerilis (u), O. adherens (u), and O. hartmanni (r) due to the closer proximity of their thermal optimum to their upper thermal limits for reproduction, as confirmed (see below) by our simulations.
The spatial projections of the HSI captured the overall latitudinal range, where this group of annelids is found. This correspondence is better represented when the fundamental niche is defined by the thermal range for reproduction (0.8 < HSI < 1) rather than the range for physiological thermal tolerance (0.1 < HSI < 0.5). The former is in fact expected to provide a more accurate approximation of the conditions that shape the ecological success of a species (Helaouët and Beaugrand 2009). In both ways, O. labronica (c) is the only species shown to be a thermal-range conformer: it exhibits a close match between its latitudinal range and its thermal tolerance (Sunday et al. 2012). The idea of the presence of a 'hyper' performing species within the annelid genus Ophryotrocha echoed what was observed in the widespread diving beetle Deronectes latus Stephens, 1829 (Calosi et al. 2008a(Calosi et al. , 2010 and the gammarid Gammarus duebeni Lilljeborg, 1852 (Gaston and Spicer 2001). Both D. latus and G. duebeni are the most widespread and physiologically tolerant species of their genus, supporting Gaston and Spicer's (2001) idea of the existence of species that are jacks-of-all-trades and masters-of-all, further confirmed by observations on the thermal metabolic performance and tolerance plasticity of the invasive decapod Asian shrimp Palaemon macrodactylus Rath-bun, 1902 (Magozzi andCalosi 2015). Conversely, all of the other species investigated in our study did not meet their potential distribution range. The mismatch between the geographic range of a species and its potentially occupied range based on its thermal physiology can reduce the accuracy with which projections that are based only on temperature can generate predictions on range shifts (Aubry et al. 2017).
Our simulations foresee two trends of HSI shifts extensively documented in temperate marine ectotherms in response to ocean warming: (i) range expansion towards higher latitudes and (ii) range contraction at lower latitudes (Cheung et al. 2009;Morley et al. 2018). These patterns varied depending on the common-rare status of the studied species. In the reference decades used in this study, uncommon and rare species commonly lived closer to their thermal optimum for reproduction and at many more sites than common species. This pattern is projected to intensify at the end of the twenty-first Century, particularly at the lowest latitudes of the Atlantic and Mediterranean regions, where uncommon and rare species are projected to live at their optimal temperatures for performance in most localities, or above it in some regions at the warm edges of their potential distribution under RCP8.5. Here, we identified hotspots of potential population decline for O. adherens (u), O. puerilis (u), and O. hartmanni (r), the three species that showed narrower upper safety margins for niche adjustment. Some of these areas, such as the Levantine Sea, are already known for their high occurrence of extinction events of endemic species (Rilov 2016). It must be noted, however, that such a scenario seems to be potentially relevant only for O. adherens, the only species recorded in this area. Indeed, our projections suggest that ocean warming may have an overall beneficial effect on the persistence of these annelid species, as indicated by the increased HSI predicted for all the species across their potential distribution. Such predictions may be more optimistic, as they were built on the thermal responses of strains subjected to laboratory conditioning over multiple generations. Exposure to acute, severe conditions in subtidal habitats can in fact exert a positive influence on an organism's thermal performance via the activation of mechanisms regulating heat stress responses (Giomi et al. 2016). Consequently, we cannot rule out the possibility that our studied species may be able to withstand warming waters more efficiently in their natural habitat compared to an experimental setting. Nevertheless, the magnitude of these considerations must be carefully examined in light of the limitations of the methodological approach used and the ecology of the studied species. First, our simulations are based on scenarios of ocean warming that rely on yearly averaged temperatures that do not accurately predict the thermal conditions on which coastal species set their thermal optima (Boersma et al. 2016;Stuart-Smith et al. 2017). This is particularly true in subtidal habitats, where temperature can vary greatly at a short temporal scale and where the expected increase in the frequency, intensity, and duration of extreme heat events may pose an additional challenge to the persistence of less tolerant species. This highlights the importance of capturing the heterogeneity of thermal conditions at a finer scale to identify the presence of microrefugia for thermally sensitive species (Potter et al. 2013). However, temperature records that are more relevant to the scale or habitat under investigation are often unavailable or not easily accessible to the scientific community. Second, poorer dispersers are expected to have a lower probability of reaching new suitable habitats and expanding their geographic range towards higher latitudes (Thomas et al. 2004). Species of the genus Ophryotrocha can actively disperse by drifting, though presumably not for long distances given the presence of direct developing larvae. Scant dispersal abilities may limit the potential poleward expansion projected for these annelid species, especially for populations found in semi-enclosed basins, such as the Mediterranean Sea. If tracking more favourable thermal conditions is not an option, climate-driven local declines or extinction events may be avoided through plastic or adaptive behavioural and/or physiological responses. Subtidal organisms are known to possess physiological adaptations to cope with thermally stressful events, either through acclimatisation or the evolution of higher thermal tolerance maxima and thermal performance breadth (Tomanek 2010;da Silva et al. 2019). In our simulations, the inclusion of species' physiological and reproductive thermal plasticity contributed to an expanded projected thermal habitat suitability for all species but did not prevent the occurrence of HSI decreases in some of the uncommon and rare species. This result suggests that the acclimation potential of a species does not guarantee its persistence in a climate-change context, as previously reported by other authors (e.g., Gunderson and Stillman 2015;Magozzi and Calosi 2015).
In conclusion, we provide general support for the possible use of the concepts of commonness and rarity to predict the susceptibility to rapid environmental changes in temperate benthic invertebrates, particularly in marine annelids, a ubiquitous and diverse group of metazoans that are unrepresented in macroecological and macrophysiological studies. By combining experimental observations and macroecological-physiological approaches, we found a greater thermal sensitivity in uncommon-rare species compared to their common relatives, although information on those species' geographic distribution was scant, patchy, or discontinuous. However, while the experimental observations and simulations based on the common species' thermal performance was in line with theoretical expectations, predictions on the fate of uncommon and rare temperate species under ocean warming scenarios were less accurate and more susceptible to methodological limitations. This highlights the need to develop an even more rigorous definition of the concept of rarity and its intermediate levels for integration with scenarios of climate change that are more relevant for the species under investigation if we want to more accurately predict the vulnerability of 'elusive' taxa to climate change.

Author contribution statement
The experimental design has been conceived and planned by GM-N. Experimental measurements were carried out by GM-N with the support of RS and DP. GM-N conducted statistical analyses with advice from RS and PC. Data modelling was conceived by GM-N, GR, WWLC, PC and performed by GR. GM-N wrote the first draft of this manuscript supported by PC. All authors contributed to the final version of the manuscript.
Funding This work was funded by the European Union through the Horizon-2020 Marie Curie Individual Fellowship awarded to , and by the PRIN grant of the Ministero Italiano dell'Universitá e della Ricerca awarded to RS (No. YS7AP3). PC was