Understanding complex biogeographic responses to climate change

Predicting the extent and direction of species’ range shifts is a major priority for scientists and resource managers. Seminal studies have fostered the notion that biological systems responding to climate change-impacted variables (e.g., temperature, precipitation) should exhibit poleward range shifts but shifts contrary to that expectation have been frequently reported. Understanding whether those shifts are indeed contrary to climate change predictions involves understanding the most basic mechanisms determining the distribution of species. We assessed the patterns of ecologically relevant temperature metrics (e.g., daily range, min, max) along the European Atlantic coast. Temperature metrics have contrasting geographical patterns and latitude or the grand mean are poor predictors for many of them. Our data suggest that unless the appropriate metrics are analysed, the impact of climate change in even a single metric of a single stressor may lead to range shifts in directions that would otherwise be classified as “contrary to prediction”.

Changes in the distributional ranges of species are among the expected outcomes of climate change 1 . Several comprehensive studies report a broad prevalence of range shifts at poleward or upper range boundaries consistent with climate change predictions [2][3][4][5][6][7] . Still, shifts contrary to that expectation have been frequently reported [8][9][10][11] . Range shifts contrary to predictions may occur because organisms are responding to a different variable (related or unrelated to climate change), or because the predicted direction was wrongly established to begin with. In fact, by identifying appropriate controlling stressors, and refining how the predicted direction of change is established, recent analyses have shown that species may be tracking climate change even when distribution ranges are shifting in otherwise unexpected directions [12][13][14] . Making sense of shifts contrary to predictions is important as it may impact the confidence level of climate change attribution, and therefore influence public opinion and policy.
It is well recognized that the mean is a metric that oversimplifies much of the complexity of stressors [15][16][17][18] , and specific aspects of some stressors (e.g., minimum water temperature during winter) have been identified as playing key roles in determining biogeographic patterns 19,20 . However, since environmental data at the appropriate temporal and spatial scales are often lacking, the patterns of many aspects of stressors remain largely uncharacterized with appropriate detail. To fill this void, general perception seems to hold that whatever spatial gradient is detected in the mean must be reflected to a large extent by all other metrics, especially for cases where the mean neatly fits preconceived assumptions. There is, however, no statistical ground supporting that view, and we argue that this misunderstanding may lead to improper estimation of what exactly is the predicted direction of range shifts of a particular species as a response to changes in the patterns of a single stressor. This effect should be especially noticeable for systems following Liebig's law of the minimum to some extent, which emphasizes the role of the scarcest Scientific RepoRts | 5:12930 | DOi: 10.1038/srep12930 resource (or, in this case, the least favourable bioclimatic variable) in determining habitat suitability at a given location and time.
Using temperature extremes as an example, the present study aims at characterising the contrasting patterns encapsulated within a single stressor, as well as showing that climate change-induced alterations of aspects of that stressor can potentially lead to otherwise non-intuitive range shifts.

Results and Discussion
We used a dataset comprised of 90 individual 4-year-long temperature time series (six microhabitats on 15 shores from 37 °N to 55 °N latitude, Fig. 1a) to evaluate to what extent the patterns of ecologically relevant metrics of a stressor are indeed captured by that stressor's mean. Extreme temperature is a major stressor in most ecosystems, and especially in the rocky intertidal 16 , where animals and plants a few centimetres apart can be experiencing dramatic differences in body temperature 21,22 . The overall mean temperature, which is often used as a reference for estimating species distributions 17 , was shown to match latitude surprisingly well (R = 0.98, Fig. 1b), thus reinforcing the idea of a relatively smooth, continuous gradient from warmer to colder temperatures with increasing latitudes along the European Atlantic coastline. However, with the exception of 'winter mean' and 'winter 95 th percentile' (Fig. 1h,i, blue lines), all other metrics ('summer mean' , 'summer 95 th percentile' , and winter and summer '7 day mean' , 'daily range' , 'microhabitat range' , 'minimum' , '5 th percentile' and 'maximum') exhibited patterns deviating substantially from that of the grand mean. These differences highlight the key role played is grand mean, calculated using all data from each shore. Red and blue lines (c-j) calculated using the warmest and coldest 30 days of each year (7 days for (c)), per shore. The shaded area is the pattern expected if each metric was perfectly correlated with latitude. Points in shaded area are "cooler than expected given latitude", and points outside shaded area are "hotter than expected". Correlation coefficients between each metric and latitude are depicted in the top right corner of each panel (blue for cold and red for warm periods). Map created in R 35 using Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) coastline data.
Scientific RepoRts | 5:12930 | DOi: 10.1038/srep12930 by climatic, geomorphologic and oceanographic factors at the local level, and more importantly, show that such factors can change skewness or kurtosis of the distribution of temperatures (or even causing multimodality) without affecting the mean. For example, seasonality, which is strongest within the Bay of Biscay -shores H and I -appears not to drive the grand mean for these shores too far away from the expected value given their latitude, but results in remarkably high summer temperatures (in fact the highest recorded in the study area; Fig. 1i,j, red lines) and equally remarkably low temperatures during winter (again, the lowest within the study area; Fig. 1f,g, blue lines). In another example, upwelling, which is typically stronger around shores F and L during summer 23 , can be seen driving 'summer daily range' and 'summer microhabitat range' (Fig. 1d,e, red lines), likely due to the co-occurrence of low water temperatures and high air temperatures. Again, this effect does not result in any important deviation of the grand mean from the latitudinal pattern for these shores, and would be largely missed if data at the appropriate spatial scale had not been collected. Additionally, the combination of regional factors and local context can result in surprising temperature distributions, such as seen at shore N. At this shore, all metrics were found to be lower or equal to the expected value based on latitude. However, the grand mean does not reflect the magnitude of this difference, especially considering that shore N is the coldest in the study area for some of the metrics calculated. The many patterns encapsulated within the distribution of values of a single stressor clearly indicates that the grand mean may largely misrepresent many other ecologically relevant aspects of that stressor. This is in accordance with previous studies 21,24 and reinforces the notion that a-priori knowledge of the physiological requirements of a species and a detailed characterisation of the thermal extremes at the study area are fundamental to ascertain the real stress landscape imposed on organisms.
Furthermore, using a theoretical example we show that complex biogeographic responses to climate change can be interpreted by using the appropriate metrics (Fig. 2). We assume that the distribution of a theoretical species is determined by a group of relevant metrics (in this case the thermal extremes measured as 'winter minimum' and 'summer 5 th percentile') and follows Liebig's law of the minimum (i.e., at each location density is dependent on the least favourable relevant metric). In the simplest form, the distribution pattern will be determined by the least favourable of a number of relevant metrics (Fig. 2a, light orange area). If climate change results in a favourable monotonic change of all metrics (Fig. 2b), the extent of suitable locations increases and a range expansion can be expected -the "general perception" poleward scenario. However, studies have highlighted that climate change not only can result in increased mean, minimum and maximum temperatures but also in increased variability -and that the exact signature of climate change varies regionally 18,[25][26][27] . In this case, if at least one metric becomes less favourable due to the increased variability, the whole distribution can be adversely affected, and an equatorward range contraction may occur (Fig. 2c). Using the metrics computed in this study it is possible to further expand the example. If the distribution of a species was found to be dependent on the interplay between extremes like 'winter minimum' and 'summer 5 th percentile' the initial distribution pattern should include a gap at shore H, and a polar range limit at shore B (Fig. 2d). If both 'winter minimum' and 'summer 5 th percentile' become warmer, a poleward range expansion can be expected (Fig. 2e), but if 'summer 5 th percentile' becomes warmer while 'winter minimum' becomes colder, the harshness of winter conditions prevail over the favourable summers and a equatorward range contraction should occur (Fig. 2f). Interestingly, in a few locations suitability would actually increase because the limiting factor was 'summer 5 th percentile' and not 'winter minimum' (shores L and N, Fig. 2h), highlighting the consequences of different mechanisms limiting species' densities across different locations 28 . The crucial point is that if field surveys were to reveal an equatorward range contraction for this species, this range shift would not be contrary to predictions, as general perception would suggest. Instead, it would be consistent with the predicted direction of change for this biological system's response to climate change, thus representing positive evidence towards the establishment of a link to climate change.
In addition, it is conceivable that some organisms' physiological requirements may include more complex interactions of aspects of a stressor than those depicted here 29,30 . For example, a mobile organism will be able to explore the various microhabitats available within a site for thermoregulation. If that organisms' physiology is found to be negatively impacted by 'summer maximum' , it is likely that the appropriate metric to study will instead be the difference between 'summer maximum' and 'summer microhabitat range' , as very high temperatures can be avoided if cooler microhabitats are available. These complex interactions between aspects of a stressor can generate new stress landscape patterns that do not match that of the grand mean, 'summer maximum' or 'summer microhabitat range' , further increasing the likelihood of erroneous expectations about the extent and direction of range shifts in face of climate change if the stress landscape is not properly characterized. The findings presented in this study reinforce the notion that using the appropriate metrics of a stressor and identifying the appropriate stressor 13 can provide decisive insights towards the detection, interpretation and prediction of complex distribution patterns, spatial and temporal variations of mechanisms controlling species distributions and the direction of range shifts. In addition, the conceptual framework here outlined emphasizes the paramount importance of coupling the collection of environmental data at the appropriate scales with a detailed characterization of species' physiological requirements (see Ashcroft et al. 31 or Greenberg et al. 32 for analogous approaches using modelled data). Although focused on the thermal regimes of the European Atlantic intertidal ecosystem, the concepts here outlined can be extended to other geographical regions, ecosystems, and stressors. Taken with necessary caution (see, for example, the cautionary advice by Scientific RepoRts | 5:12930 | DOi: 10.1038/srep12930 Austin 33 regarding assumptions of linear response to temperature), these results suggest that some of the cases where species have been shown to be shifting in directions contrary to expectations, the predicted direction of change may have been wrongly established.

Methods
Microhabitat temperature. Intertidal microhabitat temperatures were recorded at 15 exposed to moderately exposed shores along the European Atlantic coast, spanning nearly 20° of latitude, from Southwest Scotland to South Portugal (Fig. 1a, A -South Cairn, B -Emlagh, C -Holyhead, D -Annascaul, E -Wembury, F -Landunvez, G -Batz-sur-Mer, H -Royan, I -Biarritz, J -Prellezo, K -La Caridad, L -Cabo Touriñan, M -Moledo, N -São Lourenço, O -Evaristo). Data were acquired using robolimpets (autonomous temperature sensing/logging devices mimicking the visual aspect and temperature trajectories of real limpets, see Lima and Wethey 34 for details). Loggers were deployed following Seabra et al. 22 . Temperatures were sampled from 6 distinct combinations of height above the low water mark (low, mid and high shore) and exposure to sun (shaded and sun-exposed), thus covering most of the spectrum of microhabitats occupied by intertidal species. Data were collected continuously between the summers of 2010 and 2014, at a sampling rate of 60 minutes and a resolution of 0.5 °C. For each microhabitat, logged temperatures were averaged whenever data from multiple sensors were available. All data manipulation and analyses were done using R 3.1.  include the mean using all available data for each shore (grand mean), the mean temperature during the hottest/coldest seven days of each year ('mean 7d'), and the mean daily range of temperatures ('daily range'), mean daily range of all microhabitats' maximum temperatures (i.e., microhabitat range, or 'micro range'), minimum, 5 th percentile, mean, 95 th percentile and maximum during the hottest/coldest 30 days of each year. For easier terminology, metrics computed during the hottest periods were prefixed "summer", and metrics computed during the coldest periods were prefixed "winter" (e.g., 'summer minimum' , 'winter mean 7d' , etc.). The correlation coefficient between each metric and latitude was also calculated.
Direction of range shifts. An example is presented to illustrate how range shifts driven by climate change can occur both towards the poles or the equator. We modelled the relative abundance (from 0 -absent, to 1 -highest abundance) of a theoretical species under two climate change scenarios. We modelled an equatorial species which is intolerant of cold stress and tolerant of heat stress. Abundance was determined as the lowest value of either 'winter minimum' or 'summer 5 th percentile' at each shore (following Liebig's law of the minimum). To allow the comparison of both metrics (which are not equivalent in absolute terms) we normalised each metric to vary from 0 to 1, reflecting the range observed within the study region. Zero abundance occurred at shores where at least one of the metrics had a value of zero, meaning that either 'winter minimum' or 'summer 5 th percentile' , or both, prevented the species from existing. Both metrics were used without any change for the initial conditions. The first climate change scenario considers climate change as a monotonic increase of both aspects of temperature, and the abundance pattern was computed using "hot" versions of both 'winter minimum' and 'summer 5 th percentile' (resulting in increased habitat suitability). The second scenario considers climate change as an increase of variability in which some metrics may actually become colder. In this case the abundance pattern was computed using the "cold" version of 'winter minimum' and the "hot" version 'summer 5 th percentile' . The initial range limit was identified as the most poleward shore with abundance greater than zero. The location of the range limit was re-evaluated for both climate change scenarios to determine the direction of change (poleward or equatorward).