Warm and cold temperatures limit the maximum body length of teleost fishes across a latitudinal gradient in Norwegian waters

throughout Norwegian waters. Our results showed that the two largest study species, spotted wolffish ( Anarhichas minor ) and cusk ( Brosme brosme ), display the strongest negative temperature-size response. We also observed smaller maximum body lengths for multiple species within the coldest extent of their temperature range, as well as parabolic relationships between maximum length and temperature for Atlantic wolffish ( Anarhichas lupus ) and beaked redfish ( Sebastes mentella ). The smaller maximum body lengths for high latitude species at both warm and cold temperature extremes of species’ thermal ranges corroborate the temperature-size mechanisms of the gill-oxygen limitation theory (GOLT), whereby spontaneous protein denaturation limits growth at both warm and cold temperatures. and key that Here, we analysed the maximum body length of ten teleost fish species in line with temperature, dissolved oxygen concentration and geographic location (that encompasses multiple latent variables), across a broad (26°) latitudinal gradient


Introduction
Since the start of the twentieth century, mean global temperatures have increased approximately by 1.1 °C concurrently with increasing anthropogenic greenhouse gas emissions (IPCC 2021). In parallel, the geographic extent of low oxygen areas in global oceans has been expanding, whilst upper layer oxygen Abstract As the majority of marine organisms are water-breathing ectotherms, temperature and dissolved oxygen are key environmental variables that influence their fitness and geographic distribution. In line with the temperature-size rule (TSR), marine ectotherms in warmer temperatures will grow to a smaller maximum body size, yet the extent to which different species experience this temperature-size response varies. Here, we analysed the maximum body length of ten teleost fish species in line with temperature, dissolved oxygen concentration and geographic location (that encompasses multiple latent variables), across a broad (26°) latitudinal gradient 1 3 Vol:. (1234567890) levels have decreased, likely by 0.5-3.3% at surface and sub-surface layers throughout the past 50 years (Schmidtko et al. 2017;Breitburg et al. 2018;Bindoff et al. 2019). As the majority of marine species are water-breathing ectotherms, such environmental changes may present deleterious impacts on individual fitness.
In general, the response of marine water-breathing ectotherms to ocean warming and deoxygenation includes shifts in their geographic distribution to track their thermal affinities (Poloczanska et al. 2016). This may result from the complex interaction of temperature and dissolved oxygen concentration in sea water, as oxygen demands of water-breathing ectotherms increase, whilst the availability of dissolved oxygen may decrease with warming, enhanced ocean stratification and the accelerated respiration of dissolved organic matter . It is also suggested that a response to warming and deoxygenation includes the reduction of body size (Gardner et al. 2011).
Across an ectothermic species' geographical range, intraspecific clines in body size are apparent along latitudinal gradients, whereby individuals are generally smaller in warmer conditions (Horne et al. 2015). According to the temperature-size rule (TSR), those individuals reared in warmer temperatures grow faster and mature at a smaller size (Atkinson 1994). Whilst this 'temperature-size response' has been well documented in wild populations of aquatic ectotherms (Daufresne et al. 2009), the mechanisms underpinning this response have remained debated (Pauly 2010(Pauly , 2021Lefevre et al. 2017) and several explanations have been offered to resolve this interaction between environment and growth (Deutsch et al. 2015;Pörtner et al. 2017;Clarke et al. 2021;Pauly 2021;Verberk et al. 2021).
The strength of temperature-size responses will vary between species, likely based on factors including typical body size (Rubalcaba et al. 2020), sensitivity to environmental changes (Forster et al. 2011;Hoefnagel and Verberk 2015), scope of thermal niche (García Molinos et al. 2016;Burrows et al. 2019;Brito-Morales et al. 2020), life history (Weber et al. 2015;Audzijonyte et al. 2020), geography (Deutsch et al. 2020;Clarke et al. 2021), feeding strategy and behaviour (Block et al. 2011) and physiology (Atkinson et al. 2006). Therefore, rather than testing, if an ectothermic species displays a temperature-size response at all, it may be useful to compare the degree of responses between species, in order to reveal factors associated with strong temperature-size responses (Verberk et al. 2021).
In the present study, we investigated the relation between the maximum body length (hereafter simply referred to as 'maximum length') of ten teleost fish species and ocean temperature and dissolved oxygen concentration. This was done by fitting generalised additive models (GAMs) to open-source long-term bottom trawl survey data collected in Norwegian shelf and offshore waters (Djupevåg 2021 Whilst much of the temperature-size research on marine ectotherms has focused on warming and maturation, we investigated the influence of temperature, dissolved oxygen concentrations and geographic location (representing other latent environmental variables) on the maximum length of species, as the mechanisms between these responses may vary (Hoefnagel et al. 2018). In turn, we hypothesise that the largest study species will exhibit the strongest negative temperature-size response (i.e. smaller in warmer temperatures), as larger ectotherm species have been found to be more sensitive to increased temperatures and deoxygenation (Rodnick et al. 2004;Messmer et al. 2017;Rubalcaba et al. 2020).

Study area
The study area includes coastal and offshore waters of Norway, extending from the temperate North Sea, Skagerrak and the Norwegian Sea to 1 3 Vol.: (0123456789) the Barents Sea and southern portion of the Arctic Ocean along the northern shelf of Svalbard ( Fig. 1). Spanning a latitudinal range from 56°N to 82°N, the study area is marked by a temperature gradient from an annual mean sea surface temperature (SST) of 10.7 °C in the North Sea, to 4.8 °C in the northern Barents Sea and southern Arctic Ocean (Gonzalez-Pola et al. 2019).

Trawl survey data
The survey data employed in the present study were collected during the Norwegian Institute of Marine Research's long-term bottom trawl surveys between 1980 and 2020 (Djupevåg 2021). These long-term trawl surveys provide data on the abundance of targeted species to inform stock assessments across To account for variation in sampling methods over the period 1980-2020, we restricted the data utilised in the present study to a single gear type and sampling method ('shrimp trawl', 15-35 mm mesh size) between 1990 and 2020. This provided us with a total of 256,788 records of 198 unique species caught in 22,262 research trawls (Fig. 2). The mean length of trawling events was 1.2 nautical miles. We selected those species that had more than 2000 records. This reduced our dataset to 239,902 records for 36 species. The information gathered from each record included the maximum length (total length (TL), cm) of fishes collected in a trawling event, the bottom depth of the trawl sample and geographic location of the start of the trawl station (latitude, longitude).

Environmental data and exploratory analyses
Mean temperature (°C) and dissolved oxygen concentration (mols/m −3 ) at the maximum bottom depth of each raster cell (5 arcmin) were extracted from Bio-Oracle (Assis et al. 2018). Bio-Oracle environmental layers were produced via monthly averages of climate data between 2000 and 2014; thus, they represent the prevailing, long-term environmental conditions experienced by the study species. Environmental measurements were associated with maximum length records based on the latitude and longitude of sample location. For data analyses, longitude and latitude were expressed in Universal Transverse Mercator (UTM) coordinates, i.e. eastings and northings. For reporting purposes, the dissolved oxygen concentration values from Bio-Oracle were converted from mols/m −3 to mg/L using the unit conversions provided by the International Council for the Exploration of the Seas (ICES 2022).
We tested for multicollinearity between environmental covariates by calculating the Pearson's correlation coefficient (r) between temperature, dissolved oxygen concentration, eastings and northings for the 36 selected species. Such multicollinearity analyses are necessary as regression models, including GAMs, are sensitive to highly correlated continuous variables (Guisan et al. 2002;Dormann et al. 2013). For the majority of the 36 species, the Pearson's r between temperature and dissolved oxygen concentration was greater than 0.7. In order to avoid problems with model fit via multicollinearity, the relationship between temperature and dissolved oxygen concentration across the study extent was examined using a GAM, and the residuals from this model (henceforth 'OxyResid') were extracted and subsequently used as a covariate within our maximum length GAMs (Supplementary Fig. 1) (Leathwick et al. 2006;García et al. 2020). The OxyResid covariate indicates, at each sample site, the deviation in mean dissolved oxygen concentration that is expected at its given temperature, i.e. accounts for regional-scale variation in dissolved oxygen concentration (Leathwick et al. 2006).
The Pearson's correlation coefficients were also calculated between temperature, OxyResid, eastings and northings for the environmental data pertaining to where each species was sampled. Pearson's r's greater than 0.7 in absolute value were considered indicative of multicollinearity, and species where this was the case were excluded from subsequent analyses. A total of ten species (67,029 records) had a Pearson's r less than 0.7 in absolute value and were, therefore, included in maximum length GAM analyses (Table 1). Before fitting GAMs, we explored the distribution of the covariates included to visually identify outliers and subsequently removed them (Zuur et al. 2010). After this data treatment, the total number of records amongst the ten study species was reduced to 66,875 .

Study species
The ten species selected for maximum length GAM analyses include two pelagic species, one bathypelagic species, one benthopelagic species and six demersal species. Capelin is a pelagic, pan-Arctic planktivore species that represents an important prey item for piscivorous fish (Fall et al. 2018). Golden redfish is a pelagic, long-lived and late-to-mature Atlantic species that is found across shelf waters and the upper bathyal zone (Bakay 2017). Beaked redfish is a bathypelagic, long-lived species that inhabits deep-water habitats, as well as pelagic areas of continental slopes of the Atlantic Ocean (Cadrin et al. 2010). Greenland halibut is a benthopelagic flatfish species with circumpolar distribution that inhabits coastal fjords and bays, as well as slope areas of continental shelves (Giraldo et al. 2018). Norway redfish is a demersal species that occurs in Atlantic coastal and shelf waters and prefers structurally complex habitat, including boulder fields or sponge and coral habitat (Kutti and Fosså 2015). Cusk is a demersal Atlantic species that inhabits the upper slope to deep continental shelf waters (Knutsen et al. 2009). Atlantic wolffish is a demersal, denning species that inhabit high-rugosity habitat in coastal waters, as well as gravelly and boulder field habitats in deeper waters (Novaczek et al. 2017). Spotted wolffish is demersal, and similarly inhabits complex and high-rugosity habitat around or within sheltered dens across the Atlantic Ocean (Baker et al. 2012). Daubed shanny is a demersal, circumpolar, cold-adapted species that inhabits soft bottoms (Meyer Ottesen et al. 2011). Lastly, polar cod is a circumpolar demersal species that is the most abundant Arctic fish to inhabit areas under and around pack ice, although it is also found in ice-free areas (David et al. 2016).

Model fitting
GAMs were employed to investigate the influence of temperature and dissolved oxygen concentration (more precisely, OxyResid) on the maximum length of the ten study species. In addition to these environmental covariates, a tensor product smooth between eastings and northings was included to represent the effects of geographic location. This term represents a proxy for latent environmental or ecological covariates that were otherwise not included in the model (e.g. variation in depth). This tensor product smooth term also accounts for the broad-scale spatial autocorrelation that may be present within the maximum length data (Wood 2006;Grüss et al. 2016Grüss et al. , 2018aGrüss et al. , 2021. For all ten study species, GAMs were fitted to maximum length data using the 'gam' function from the mgcv package in the R environment (Wood 2006; Team R Development Core 2021) as follows, where L is the maximum length; g is the log-link function between maximum length and each term on the right side of the equation; te(X,Y) is the tensor product smooth between eastings and northings; and s is the thin plate regression spline with shrinkage (bs = 'ts' specified in the smooth function of 'gam' within the mgcv library) fitted to temperature and OxyResid (Grüss et al. 2018a).
The Gamma distribution was identified as the most appropriate distribution for the maximum length data employed in the present study, using diagnostics from the fitdistrplus package (Delignette-Muller and Dutang 2015) in R (Team R Development Core 2021). The thin plate regression splines employed in GAMs were limited to four degrees of freedom (k = 4) in order to prevent overfitting and to preserve the interpretability of model results (Mannocci et al. 2017;Grüss et al. 2018b;Weijerman et al. 2019). In addition, an extra penalty was applied to model covariates as their smoothing parameter approached zero, which allowed for the complete removal of a covariate when its smoothing parameter equalled zero (Marra and Wood 2011;Grüss et al. 2014).

Model validation
To evaluate the GAMs of the ten study species, we relied on two performance metrics: (1) the adjusted coefficient of determination (adjusted-R 2 ), which indicates the percentage of the variance explained by the GAMs; and (2) Spearman's rank correlation coefficients (Spearman's rho's) between predicted maximum lengths and observed maximum lengths. Our validation procedure consisted, for each GAM, of the following: (1) the dataset of interest was split into a training dataset which was randomly assigned 70% of the data and a test dataset which was assigned the remaining 30% of the data; (2) this was repeated 10 times, such that 10 training datasets and 10 corresponding test datasets were produced; (3) the GAM was fit to each of the 10 training datasets; and (4) the 10 GAMs fitted using the 10 training datasets were evaluated using the corresponding test datasets and the performance metrics (adjusted-R 2 and Spearman's rho's) (Grüss et al. 2016;Egerton et al. 2021

Analyses with the fitted and validated models
To determine the relative importance of temperature, OxyResid and geographic location in explaining the maximum length of the ten study species, we implemented the approach of Grüss et al. (2016). This approach consists, for each study species, in comparing the predictions of the final GAM with those of GAMs where the values of a given predictor (temperature, OxyResid, eastings or northings) are randomly permuted within the dataset fitted to the model ('random GAMs') (Thuiller et al. 2012;Grüss et al. 2016). An index of relative importance is obtained for each predictor (temperature, OxyResid, eastings and northings) by computing one minus the Pearson's correlation coefficient between the predictions of the final GAMs and the predictions from the random GAMs (Grüss et al. 2016;Dove et al. 2020;Bolser et al. 2020). We also used the results from the GAMs to predict the maximum length of study species as a function of temperature using the 'predict.gam' function from the mgcv package (Wood 2006; Team R Development Core 2021). Specifically, we predicted maximum lengths over a vector of values ranging between the minimum and maximum values of temperature encountered by the species of interest where it was sampled, whilst (1) setting OxyResid to its mean value from the modelled dataset and (2) setting easting and northings to their values at the barycentre of the study area (Grüss et al. 2018c.
Examining the potential effects of fishing pressure on the maximum length of study species In order to acknowledge the potential effects of fishing pressure on the maximum length of our study species, we calculated a 'pseudo-Z', or pseudo-total mortality (Z) value, for the 9 species included in the present study that are targeted by fisheries (all species except for daubed shanny). This was done by employing our maximum length records to construct a length-converted catch curve (Pauly 1990), using the FiSAT II Fish Stock Assessment Tool (Gayanilo et al. 2005). Catch curves were calculated using two growth parameters included in von Bertalanffy growth functions: asymptotic length (L ∞ ) and growth coefficient (K, see Supplementary Table 1 and Appendix for methods). For each species, pseudo-Z values were calculated across 5-year periods. We then compared these values to determine whether trends in fishing mortality (F) changed across the 30-year study period. Here, we assume that natural mortality (M) for species is constant across time, so that changes in pseudo-Z reflect changes in F (Z = M + F; Beverton and Holt 1956).

Results
The length data included in analyses were generally representative of adult size classes for study species (Table 1). Across the study area, mean sea surface temperature ranged from − 1 in the northern extent to 11.5 °C in the southern extent, whilst mean temperature at bottom depth ranged from − 1 to 8.9 °C ( Supplementary Fig. 12). The gradient of dissolved oxygen concentration ranged from 8.45 in the southern extent to 10.9 mg/L in the northern extent of the study area ( Supplementary Fig. 13).
For all ten study species, temperature, Oxy-Resid and the interaction between eastings and northings were found to have a significant 1 3 Vol:. (1234567890) effect (p < 0.005) on their maximum length (Table 2). In estimating the relative importance of model covariates in explaining maximum length, eastings and northings were found to be the most important in explaining maximum length for all ten study species (Supplementary Fig. 14). All GAMs passed the model evaluation test, as (1) median adjusted-R 2 values ranged between 0.124 (spotted wolffish) and 0.331 (Greenland halibut); and (2) Spearman's rho values ranged between 0.361 (spotted wolffish) and 0.573 (beaked redfish) and were all significantly different from zero at the α = 0.05 level (Supplementary Fig. 15).
Predictions of the maximum length of the study species across their temperature range displayed three patterns (Fig. 3). The two species with the largest predicted maximum lengths, spotted wolffish and cusk, were predicted to decrease in maximum length as temperature increased-this was also the case for Norway redfish (Fig. 3). In contrast, capelin, Greenland halibut, golden redfish, daubed shanny and polar cod were predicted to increase in maximum length as temperature increased in the study area. Atlantic wolffish and beaked redfish were predicted to increase in maximum length until a given temperature, after which maximum length decreased (Fig. 3). Lastly, the fishing mortality proxy pseudo-Z values suggested that no observable positive or negative trends occurred across the 30-year study period (Supplementary Table 1).

Discussion
Overall, our results confirm our original hypothesis that large species will exhibit the strongest negative temperature-size response. The two largest species, spotted wolffish and cusk, were predicted to experience the largest relative decrease in maximum body length across their observed temperature range (Fig. 3). These results are consistent with the body of literature that suggests larger species are more likely to experience stronger negative temperature-size responses (Forster et al. 2012;van Rijn et al. 2017). In contrast, the maximum lengths of capelin, Greenland halibut, golden redfish, daubed shanny and polar cod were shown to increase with temperature. Interestingly, all of the aforementioned species display smaller maximum lengths within the coldest extent of their temperature range.
Our results corroborate the temperature-size mechanism that is part of the gill-oxygen limitation theory (GOLT, Pauly 2019Pauly , 2021. According to the GOLT, as the 3-dimensional bodies of ectotherms grow towards their asymptotic sizes, their 2-dimensional gills become gradually unable to supply sufficient oxygen (per unit weight) required for the synthesis of body proteins to exceed the rate of spontaneous protein denaturation (Pauly 1979(Pauly , 1981. As a result, growth ceases when synthesis equals denaturation and, thus, relatively high temperatures for the thermal niche of a species, which promote denaturation, cause ectotherms to remain smaller (Pauly and Cheung 2018;Pauly 2019Pauly , 2021, as exhibited here by spotted wolffish, cusk and Norway redfish. Table 2 Number of data points fitted to the generalised additive models (GAMs) of the ten study species, and the adjusted-R 2 of these GAMs. All GAMs express maximum body length as a function of temperature, the residuals of dissolved oxygen and temperature (OxyResid) and an interaction term (tensor product smooth) between eastings and northings (Eq. 1)

Species
Number of data points  Fig. 3 The maximum length (TL, cm) of the ten species as a function of temperature (°C), as predicted by the generalised additive models (GAMs) developed in the present study, with 95% confidence intervals. The dashed line indicates the reported temperature range of species, including the mean (dot) as reported by species-specific mapping parameters in AquaMaps (Kaschner et al. 2019). This figure was created using R Studio (Team R Development Core 2021) Environ Biol Fish (2022) 105:1415-1429 1423 Spontaneous protein denaturation, i.e. the loss of the quaternary structure which give proteins their shapes, and thus enable them to function (e.g. as enzyme), however, is not only accelerated by increasing, but also by decreasing temperature, a process known as 'cold denaturation' (Privalov 1990;Graziano 2014;Sanfelice et al. 2015;Yan et al. 2018). The result has been called 'metabolic cold adaptation', whereby some species have evolved to elevate their metabolic rate to live at cold temperatures (Wohlschlag 1962(Wohlschlag , 1964Pauly 1979). Thus, because of the need to replace proteins denatured by the cold temperature they inhabit, some polar fishes display a metabolic rate that is higher than expected given a extrapolation of standard metabolism-temperature relationships (Wohlschlag 1962(Wohlschlag , 1964Pauly 1979). Also, Pauly (1980) reported that fishes in water < 4 °C exhibited higher natural mortality rates than expected, other things (i.e. growth parameters) being equal, and body size was significantly negatively correlated with mortality. These results suggest that cold-water fishes are relatively smaller in body size (Pauly 1980). Cold denaturation has received relatively little attention vis-à-vis the temperature-size relationship for water-breathing ectotherms (Todgham et al. 2007(Todgham et al. , 2017Fraser et al. 2007;Peck 2016Peck , 2018. Although it was a constitutive element of the first presentation of the GOLT (Pauly 1979) and was corroborated in Pauly (1980), it was not incorporated into the further development of this theory. In fact, Pauly (2019, pg. 169) sent it into 'profound hibernation'. However, the humped shapes of the temperature-size response for Atlantic wolffish and beaked redfish present evidence for cold denaturation causing a response similar to elevated temperatures, thus justifying reviving the cold denaturation hypothesis as an explanatory mechanism for reduced maximum sizes of fish at very low temperatures. Recent findings from the Antarctic also support the reviving of the cold denaturation hypothesis. Fraser et al. (2022) found that the growth of Antarctic fishes was limited by high 'degradation' rates of proteins relative to temperate species. These and our results add to the growing body of literature highlighting the limitations of growth for marine ectotherms at low temperatures.
As previously mentioned, capelin, Greenland halibut, golden redfish, daubed shanny and polar cod are predicted to increase in maximum length with temperature (Fig. 3). These findings suggest that these species have optimal temperatures which exceed the temperatures at which they were sampled, and beyond, where their maximum length would start to decline, as the data employed in the current study only covers a part of the temperatures our study species experience (Fig. 3). We interpret the parabolic temperature-size response of Atlantic wolffish and beaked redfish, to result from the study area spanning a range of temperatures which include the temperature to which they are best adapted, and thus their maximum length can increase or decrease depending on the local conditions. Overall, we expect that the variable temperature-size responses at both warm and cold extents reflect species-specific growth sensitivities to either warming or cooling temperatures across their geographic range, as well as species traits. We, therefore, recommend applications of our modelling approach that would incorporate further species' traits, including their sensitivities to temperature and dissolved oxygen concentration, to elucidate temperaturesize responses at both warm and cold temperatures between species of variable size, milieu, activity and trophic level (van Rijn et al. 2017).
Regarding the relative importance of model covariates, for all study species, geographic location was found to be most important in explaining maximum length. This is to be expected, as these variables were chosen to represent latent environmental and ecological variables in the study region (Grüss et al. 2016(Grüss et al. , 2018aBolser et al. 2020), and may reflect variation in depth, the availability of resources, ecological competition, size selective predation or mortality, as well as fishing effects (Connolly and Roughgarden 1999;Pauly 2010;Verberk and Bilton 2011;Cheung et al. 2013;Tu et al. 2018;Grüss et al. 2021). Regarding temperature and dissolved oxygen concentration, we recommend that future research compares the temperature-size response between taxonomic groups, as well as the interactive effect between temperature and dissolved oxygen concentration on determining the maximum length of marine ectotherms. For example, those species that are more vulnerable to low oxygen conditions may experience a stronger negative temperature-size response (Rubalcaba et al. 2020).
The maximum body lengths analysed in the current study approximate adult body sizes for all study species (Table 1). Overall, fishing pressure may truncate the size and/or age distribution of fished populations, and conversely, reductions in fishing pressure relative to historical levels could see the increase of maximum body lengths of targeted populations (Worm et al. 2009). Yet, considering the lack of any trend observed in the pseudo-Z values, we assume that fishing pressure was unlikely to confound the predictions of maximum length for the study species that are fished within the area (Supplementary Table 1). Nonetheless, the interaction of further ocean warming and fishing pressure may affect the maximum body length of targeted populations (Tu et al. 2018). Therefore, we recommend that the future studies that will utilise our modelling framework to investigate the relationship between maximum length, temperature and dissolved oxygen concentration also consider fishing pressure (e.g. trawl footprint).
In conclusion, by predicting the maximum length of various species across Norwegian waters, our results confirm our hypothesis that the largest species experience the strongest negative temperature-size responses. Furthermore, we have observed limitations to the maximum body length of several species within the coldest extent of their observed temperature range. These results support the temperature-size mechanism as described by the GOLT, and offer evidence to revive the cold denaturation hypothesis as a mechanism that limits the maximum body size of marine ectotherms at very low temperatures.
Author contribution All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Charles P. Lavin, Cesc Gordó-Vilaseca, Fabrice Stephenson and Arnaud Grüss. The first draft of the manuscript was written by Charles P. Lavin and all authors contributed to writing of the manuscript. All authors read and approved the final manuscript.

Funding Open access funding provided by Nord University
Data availability Species' length data analysed during the current study are available on GitHub at https:// github. com/ charl es-patri ck-lavin/ Norway-maxim um-length-EBFI-2022-. git.

Code availability
The code used for analysis in the current study is available on GitHub at https:// github. com/ charl es-patri ck-lavin/ Norway-maxim um-length-EBFI-2022-. git. Code for specific functions is available upon request.

Declarations
Ethics approval Not applicable.

Consent to participate Not applicable.
Consent for publication Not applicable.

Conflict of interest
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.