Differential responses of trailing-edge populations of a foundation alga to thermal stress

ABSTRACT

Populations within the same species range edge may experience contrasting local conditions and exhibit diverse levels of environmental tolerance. This heterogeneity within a range boundary has seldom been considered in studies forecasting the impact of anthropogenic habitat alteration and climate change on species distributions. Moreover, any ecological prediction under changing environmental conditions requires a good understanding of the combined responses of organisms to multiple stressors, in particular the effects on key life cycle stages of species. The intertidal seaweed Fucus serratus is a dominant species on northern Atlantic shores and whose southern limit is in the NW Iberian Peninsula. We examined how early developmental stages of southern-edge populations of this foundation alga responded to the combined effects of environmental stressors, including salinity and aerial and seawater temperature. Four populations from two different areas of the NW Iberian coast were considered: two populations from open shores on the northern coast and two populations from the western rias on the Atlantic coast. The study findings revealed inter-population variability in the response to heat stress, with greater survival of germlings from the northern populations than those from western rias. Environmental conditions are still more benign within western rias, under the influence of strong summer upwelling events. The results also suggest the key role of aerial thermal stress in determining the southern limit of distribution of the target species. The future of these heterogeneous edge populations from NW Iberian Peninsula depends on trends in climate change and the ability of populations to cope with these. Environmental changes may already be occurring at rates that exceed the plastic and adaptive potential of edge populations in N Spain, while the future of western rias as climate refugia for this and other foundation alga is also uncertain.

KEYWORDS:

• Climate change
• climatic refugia
• foundation species
• intertidal
• phenotypic differentiation
• rocky shores

Introduction

It has long been recognized that species are not homogeneous entities, as conspecific populations exposed to different environmental conditions are often phenotypically and/or genetically distinct (Linhart & Grant, 1996). However, the inter-population variability in environmental affinities and tolerances has only recently been considered in studies forecasting the impact of anthropogenic habitat alteration and climate change on species distributions. Thus, these projections often consider the whole-species level, assumed to be a single entity with ecologically similar members (but see for instance Banta et al., 2012; Nicastro et al. 2013; Valladares et al., 2014; Bennett et al., 2015; Saada et al., 2016).

In this context, geographically peripheral populations have frequently been considered to live at the margin of their environmental tolerances and to be more susceptible to extinction than central populations, with which they share the same range of environmental tolerance (Sagarin & Gaines, 2002 and references therein). However, edge populations, in particular those at the low-latitude boundary of species distributions (rear-edge sensu Hampe & Petit, 2005) may have unique phenotypic and genetic traits and may be acclimatized or adapted to these stressful marginal environments (Kawecki, 2008; Bennett et al., 2015; Rehm et al., 2015). In addition, geographically peripheral populations are not always located in unfavourable, ecologically marginal sites (Pironon et al., 2017) as they can occupy ‘core-like’ locations, i.e. contemporary climatic refugia (Ashcroft, 2010; Keppel et al., 2012). Different populations within the same species range edge may experience contrasting local conditions and exhibit diverse levels of environmental tolerance. This edge-edge variability has seldom been considered, despite its importance in range limit dynamics and species responses to climate change.

A prerequisite to constructing any ecological prediction under changing environmental conditions is a good understanding of how organisms respond to interacting environmental factors (Gunderson et al., 2016). The reaction to combined stressors, which may be additive, synergistic or antagonistic, will in turn determine the impact on the performance of individuals and populations (Crain et al., 2008). The effects of different environmental components on key life cycle stages of species and how such effects vary among populations is a current priority in ecological research.

Marine organisms living in intertidal habitats are exposed to multiple environmental stresses, linked to both aquatic and aerial regimes, and they may thus be considered early tracers of the effects of changing climate conditions (Helmuth et al., 2006). Canopy-forming macroalgae play a key role in temperate coastal ecosystems by providing food and shelter for diverse species of fauna and flora, as well as sustaining complex food webs (Chapman, 1995; Schiel & Foster, 2006). Retraction of the low-latitude, rear edge of these dominant habitat-forming seaweeds is predicted in northern Atlantic rocky shores (Jueterbock et al., 2013) and many of these species are already in decline in this region (e.g. Lima et al., 2007; Díez et al., 2012; Nicastro et al. 2013; Fernández, 2016; Assis et al., 2017). However, the inter-population variability in environmental affinities and specifically the potential of southern edge populations to withstand environmental changes remain poorly investigated (Ferreira et al., 2014; Jueterbock et al., 2014; Bennett et al., 2015, Pearson et al., 2009). This particularly applies to how early developmental stages, which are usually bottlenecks in seaweed populations due to high mortality (Santelices, 1990; Vadas et al., 1992), respond to multiple stressors, as most of these inter-population studies focus on adult thalli.

The seaweed Fucus serratus L is a dominant canopy-forming species on northern Atlantic intertidal shores and its southern limit is in the NW Iberian Peninsula (Fischer-Piette, 1955; Lüning, 1990). The species has two range edges in this area, one in N Spain and the other in N Portugal (Lüning, 1990; Arrontes, 1993). In the last few decades, the distribution of F. serratus in N Spain has shifted towards the west, and its presence is now almost entirely limited to scattered populations in wave-sheltered Atlantic rias and a few semi-exposed populations in northern Spain, mostly, on the coast of the province of Lugo (Duarte et al., 2013; Araújo et al., 2014). Marginal populations from western Atlantic rias and those of N Spain are exposed to contrasting environmental conditions. In addition to differences in wave exposure, cooler water and greater nutrient supply in summer, and higher salinity variability over the year were detected inside western rias than in northern semi-exposed shores, due to the stronger influence of upwelling events and river flows (Martinez et al., 2012a, Abrantes et al., 2017; Duarte & Viejo, 2018 and references therein). Furthermore, environmental conditions are changing rapidly, particularly in northern coastal areas, where seawater and air temperatures are increasing (Gómez-Gesteira et al., 2008; Meinshausen et al., 2009; Abrantes et al., 2017) and the surface salinity is decreasing (González-Pola et al., 2005; Llope et al., 2006). The interactive effects of physical factors may affect these marginal populations in different ways, and the ability of the populations to persist will also vary. Only one previous study has investigated the combined effects of several physical factors in two marginal populations of F. serratus in experiments with vegetative fronds of the species (Martínez et al., 2012b).

We evaluated how early developmental stages of F. serratus in populations in two marginal geographic areas, i.e. Atlantic rias and the northern coast of the NW Iberian Peninsula, respond to the combined effects of salinity, seawater and air temperature. We hypothesized that tolerance to thermal stress would be higher in populations of N Spain, presently exposed to harsh conditions (Duarte & Viejo 2018), while populations from rias, under the influence of freshwater inflow, would be more tolerant to low salinity levels. We also expected strong effects of air temperature in the performance of this intertidal species. In order to determine stage-specific variations in the vulnerability of the alga to environmental factors, the vital rates of two different developmental stages (germlings, juveniles) were estimated in mesocosm experiments under controlled conditions. Fine-scale empirical studies like these are necessary for more accurate prediction of species responses to climate change than obtained by broad-scale modelling approaches.

Material and methods

Collection of germlings and juvenile fronds

Germlings and juvenile stages of F. serratus were collected from four locations on the coast of Galicia (NW Iberian Peninsula): two locations in the province of Lugo, on the open northern semi-exposed Cantabrian coast (Peizás: 43°35′N, 7°16′W and San Pedro: 43°37′N, 7°20′W) and two locations within large embayments or rias on the western Atlantic coast (O Freixo, in Ría de Muros: 42°47′N, 8°56′W and Isla de Arousa, in Ría de Arousa: 42°33′N, 8°51′W; Fig. 1). Discs made from epoxy resin (4.5 cm diameter, Fetadit 55/63; Fetasa, Madrid, Spain), each with a rough surface and a central hole (see Johnson, 1994), were attached to polycarbonate plates (14 × 14 cm), which were then anchored to rocks with stainless steel bolts and placed under canopies of reproductive individuals of F. serratus. The discs were placed in the four sites between June and September 2016 to allow algal settlement and subsequent quantification and collection of germlings. A total of 80 discs on 20 plates were placed within each population in the rias, and 40 discs on 10 plates as well as 26 independent discs were placed under canopies at Peizás and San Pedro, respectively. Fewer discs were placed in the sites in the province of Lugo because of lower abundance of reproductive individuals of F. serratus at these sites on northern shores (Duarte & Viejo, 2018). The discs were exposed for three months in the field and were then recovered during the spring tides in September 2016. The number of germlings on each disc was counted under a stereomicroscope. In addition, 96 juvenile individuals of F. serratus were collected from each site (mean ± SE = 0.40 ± 0.005 g fresh weight, 4.37 ± 0.04 cm length, n = 384) and transported to the laboratory in cool boxes.

Fig. 1. Study locations on the NW Iberian Peninsula (Galicia). Pz, Peizás; SP, San Pedro; OF, O Freixo (Ría de Muros); IA, Isla de Arousa (Ría de Arousa). The black arrow shows the current rear-edge distribution limits along the Iberian Peninsula

Mesocosm experiments

The germlings that settled on the artificial discs and the juvenile fronds collected were used in mesocosm experiments, conducted at the Estación de Ciencias Mariñas de Toralla (ECIMAT, Vigo, Spain, https://cim.uvigo.gal/ecimat-cim/). Prior to the experiments, discs with germlings and the juvenile fronds were immersed in seawater of salinity ~35‰ and temperature ~16°C for 9 days with a water surface irradiance of ~ 450 μmol m⁻²s⁻¹. The experiments with germlings and juveniles were then started and lasted for 28 and 24 days, respectively (during September–October 2016).

The experimental set-up included ‘population’ (4 levels: Isla de Arousa, O Freixo, Peizás and San Pedro, hereinafter IA, OF, Pz and SP, respectively), ‘salinity’ (2 levels: low and high with mean values ~22 and 35‰), ‘seawater temperature’ (2 levels: cold and warm, with mean values ~16 and 20°C), and air temperature during emersion periods, hereafter ‘emersion’ (2 levels: low and high, with maximum temperatures of about 19 and 26°C, respectively) as fixed and orthogonal factors and ‘tank’ and/or ‘disc’ as random factors (see details of the experiments and the statistical analyses below). Given the smaller number of recruitment discs placed in the Lugo populations (Pz and SP), the complete orthogonal design in the case of germlings was only established for the ria populations (IA and OF). Germlings from Lugo were only exposed to different temperatures during emersion (at fixed conditions of 35‰ of salinity and 20°C of seawater temperature).

The salinity, seawater temperature and emersion levels were chosen to represent the different environmental conditions experienced by F. serratus individuals in the study areas (Duarte & Viejo, 2018). Salinity fluctuates within the rias and decreases from 35 to ~22–29‰ due to the river supply, while the values are less variable in the semi-exposed Lugo locations, ranging between 32 and 36.5‰ throughout the year (Intecmar, Xunta de Galicia, Spain, http://www.intecmar.gal/). The mean seawater temperature is lower within the rias than on the Lugo coast during summer (Duarte & Viejo, 2018). Finally, the air temperatures selected corresponded to daily mean values recorded by data loggers located underneath a canopy of F. serratus at O Freixo, during diurnal emersion periods, both under heatwave and non-heatwave conditions during summer of three consecutive years (2011–2013, authors unpublished data).

The experiments were set up in an isothermal walk-in chamber. The mesocosm system consisted of 32 transparent plastic tanks (3 l) supplied with 50 µm sand-filtered seawater from the Ría de Vigo (42°12′N, 8°48′W) adjusted to the two experimental salinities with fresh water and renewed every two days (mean values ± SE = 22 ± 0.04‰ and 35 ± 0.01‰, n = 16). Four tanks and three juveniles per tank from each population (i.e. 12 fronds per tank) were used for each combination of salinity × seawater temperature × emersion (i.e. a total of 32 tanks and 384 juvenile individuals). In the case of germlings, three tanks and two discs per tank and population were used for each combination of environmental conditions (in populations from rias, but only for emersion treatments at the fixed values of 35‰ of salinity and 20°C of seawater temperature for Lugo populations). A total of 120 discs (48 for each population from rias and 12 for each population from Lugo) were used. The experiments with germlings (settled on discs) and juveniles were conducted simultaneously in the same tanks (see Supplementary fig. S1 for an outline of the experimental design). The reverse side of each disc was numbered, and juveniles (hanging from ropes that ran through the tank) were also tagged with numbered plastic strips for subsequent identification.

Tanks were continuously aerated, and light was supplied from above by cool white fluorescent lamps (F18 W/840; water surface irradiance: 451.9 ± 6.4 μmol m⁻²s⁻¹ mean ± SE, n = 8) with a 12:12 h light: dark photoperiod. Tanks were randomly allocated to the salinity and seawater temperature treatments to minimize the effect of their position. The cold seawater treatment was determined by the chamber temperature (16°C mean value), while for the warm seawater treatment the water was heated using titanium aquarium heaters (100 W) regulated by a set of 3 AT Controller Twin (Aqua Medic®, Germany). Seawater temperature was recorded every 5 minutes by data loggers (IButton® Data Loggers, Germany). Mean daily values ranged between 15.67 (± 0.05 SE) and 16.84°C (± 0.01) in the cold seawater temperature treatment, and between 19.82°C (± 0.03) and 21.74°C (± 0.02) in the warm treatment, n = 288. To simulate emersion conditions during low tide, recruitment discs and juveniles were manually transferred once a day from the tanks to shelves located inside the same chamber and covered with a canopy of Fucus spp. (non-reproductive fronds). Juveniles and recruitment discs from the same tank were grouped on the shelves. The emersion time varied weekly throughout the experiment, alternating between 4 and 6 h, in order to simulate the monthly variations in the tidal range. Shelves for emersion periods were provided with light from above with cool light fluorescent lamps, at the same light intensity and photoperiod as in the tanks. The irradiance values under canopy were 328.50 ± 6.77 μmol m⁻²s⁻¹ (mean ± SE, n = 6). To produce temperature profiles similar to those experienced by macroalgae during emersion in the field, the air temperature was increased gradually using ceramic infrared heaters (150 W, Exo Terra, Mansfield, Ohio, USA) positioned over the shelves. The air temperature was regulated with digital temperature controllers (AT Control System controllers, Aqua Medic®, Germany) and recorded via individual temperature sensors inserted underneath the canopy. This system enabled continuous control and recording of air temperature with an error of 0.2°C. In the 4-hour emersion cycle, air temperature increased gradually, reaching maximal values of 18.77°C and 22.50°C in the low and high emersion treatments (Fig. 2a), while in the 6-hour cycle the maximal values were 18.71°C and 25.86°C, respectively (Fig. 2b). Care was then taken to simulate the emersion conditions experienced in the field by early developmental stages of intertidal algae, i.e. cover by adult canopies, the gradual increase in aerial temperatures during each emersion period, and the weekly variability in the length of the emersion periods.

Figure 2. Air temperatures in the two emersion treatments (high and low) in the first mesocosm experiment, a) 4-hour cycle; b) 6-hour cycle. Note the different scales on the axes. Average values ± SE are shown, n = 6

Post-stress conditions of recruits

After the experiment, discs with surviving recruits (i.e. individuals visible to the naked eye) were maintained submerged in tanks under common ambient conditions of salinity (35.2‰ ± 0.08, n = 12) and seawater temperature (16°C ± 0.11, n = 12) for 32 days in the isothermal walk-in chamber. Only discs from the ria populations (IA and OF) previously exposed to low emersion conditions, were used at this stage, due to the elevated mortality experienced by germlings in some of the population and treatment combinations in the previous experiment (see Results).

Response variables

Survival of germlings, recruits and juveniles

The number of germlings/recruits per disc was calculated under a binocular microscope as the mean value from 3 random 1 cm² subsamples per disc. The proportion of live germlings/recruits was calculated at the end of the mesocosm experiment/start of the post-stress phase (after 28 days) and at the end of the post-stress phase (i.e. 60 days). The state (dead or alive) of each juvenile was determined 10 and 17 days after the start of the mesocosm experiment (22 and 29 September 2016, respectively). Dead individuals were identified by the presence of significant tissue depigmentation, broken apices and necrosis.

Growth of germlings, recruits and juveniles

Due to the small size of germlings at the beginning of the experiment, the first measurements, i.e. length measured from top of embryo body to tip of rhizoid, were taken one week after the experiment started (t = 7 days). Mean size per disc was calculated from 10 germlings, which were randomly selected and measured under a binocular microscope connected to an image analysis system (Nikon-SMZ1500, NIS Elements Basic Research 4.0). The mean size of germlings/recruits per disc was also estimated at the end of the mesocosm experiment and at the start/end of the post-stress phase.

In the case of juveniles, the fresh weight (± 0.001 g) of each individual was estimated weekly, after the specimens had been dried by gently blotting them with paper towels. The growth of each juvenile was calculated as the difference in the initial weight and the weight after 10 and 17 days. Final weights (24 days after the start of the experiment) were not considered for analysis due to the high mortality of juveniles at that time (see Results).

Data analyses

Potential differences between populations in the mean number of germlings/recruits and the size of juveniles/recruits were evaluated at the start of the mesocosm/post-stress experiments by using Generalized Mixed Models (GLMMs; Zuur et al., 2007) and Linear Mixed Models (LMMs) respectively. Population (fixed) and Tank (random) were included in the analyses as orthogonal factors. The initial size of germlings was not known, as the first measurements were taken one week after the start of the mesocosm experiment (see above).

GLMMs for count data (initial number of germlings/recruits) were fitted with a Poisson distribution of error terms and a log-link function. Data overdispersion was detected in the germling analysis (dispersion parameter, Φ: 7.42; Zuur et al., 2009) and was corrected by adding observation-level random effects, i.e. the disc was added as a random factor (OLRE; Harrison, 2014). As differences in the initial number of germlings between populations were detected in the mesocosm experiment (see Results), the possible influence of these initial densities on survival was subsequently determined by using 70th and 90th quantile regressions (Scharf et al., 1998; Cade & Noon, 2003). Goodness of fit (R¹) was calculated for each specific quantile regression according to Koenker & Machado (1999).

The effects of the physical factors and population of origin on the survival rate of germlings was analysed at the end of the mesocosm experiment by using GLMMs fitted with a binomial distribution of error terms and a logit link function. The average size of germlings per disc was analysed one week after the start and also at the end of the mesocosm experiment by using Linear Mixed Models (LMMs). The full or ‘beyond optimal’ models (sensu Zuur et al., 2009) included four orthogonal fixed factors (population, salinity, seawater temperature and emersion) and the tank × population interaction. In each case, the explanatory variable (either survival rate or size) was analysed separately for populations from rias (including all the treatment combinations) and for all populations (testing only the effects of emersion, see the mesocosm experimental set up). In the analysis of the survival of germlings for all populations, the random tank × population structure was not considered in the full model, due to the lack of a sufficient number of replicates to ensure adequate statistical power. In addition, data overdispersion was also observed in this analysis (dispersion parameter, Φ: 10.76) and was corrected by adding disc as a random factor.

Survival of juveniles 10 and 17 days after the start of the mesocosm experiment were likewise analysed by GLMMs with binomial distribution and a logit link function. The response variable was in this case State (dead vs alive; 0 vs 1), and the full model included the four previously mentioned orthogonal fixed factors, their interactions and the random factor tank. Growth of juveniles 10 days after the start of the experiment was analysed using an LMM of the same design. Three outliers were removed to achieve the data normality assumption. At t = 17 days, only juveniles exposed to the low emersion temperature treatment were considered for modelling, due to the high mortalities detected at the high emersion treatment.

At the end of the post-stress phase, the survival rate of recruits was analysed using a Generalized Linear Model (GLM), again fitted with a binomial distribution of error terms and a logit link function. As explained above, the analyses were only carried out at this phase with recruitment discs from rias populations previously exposed to the low emersion temperature. Three orthogonal fixed factors were considered (population, salinity, and seawater temperature in the previous mesocosm experiment). Tank was not included as a random factor in the full model because of the small number of tanks at this stage. Data overdispersion was detected (dispersion parameter, Φ: 10.30) and was corrected by adding observation-level random effects. The final mean size of recruits per disc was tested with a Linear Model (LM) using the same design as for survival in this post-stress phase. Prior to the analysis, data were log-transformed to comply with the assumption of data normality.

Model selection in mixed models, either GLMMs or LMMs, was performed in two steps, starting with the full (or beyond optimal) model. First, the optimal random structure was selected from nested models fitted with restricted maximum likelihood estimation by using the hypothesis testing approach. The p-values were corrected due to the problems of testing on the boundary (Zuur et al., 2009). The optimal fixed structure was then selected by comparing models fitted with maximum likelihood estimations. Given a selected random structure, models with different combination of fixed components were compared using Akaike Information Criteria for small samples (AICc; Burnham & Anderson, 2004). The differences in AICc values between each model and the model with the minimum value (Δi = AICc_i−AICc_min) were also calculated. The subset of candidate models was those with the lowest AICc values and Δi > 2 with the other models (Burnham & Anderson, 2004). The Akaike weights (ɷ_i) of the subset of candidate models was also determined, as was the relative importance of the factor in the candidate subset of models (Burnham, 2015). Akaike weights are equivalent to the probability of a given model being the best in the candidate subset of models (Burnham & Anderson, 2004). The relative importance of the factor is the sum of all the Akaike weights of the candidate subset of models containing each explanatory variable.

Data normality, homogeneity of variance and data dispersion assumptions were tested and where necessary were corrected in the selected model. Once the final model was selected, a Tukey’s post hoc test was applied to examine pairwise differences for significant interactions terms or main effects.

The complexity of the experimental design and the limited number of available populations per area (two within rias and two on northern shores), precluded the inclusion of Area as an additional main fixed factor and Population as a nested random factor in mixed models. A higher number of populations should be available (> 4 but preferable > 10; Zuur et al., 2007). We used a parsimonious approach by considering Population as a main fixed factor, determining the existence of divergent patterns between ria and northern populations by significant differences in pairwise comparisons.

Quantile regressions were conducted using the R package quantreg (Koenker et al., 2019), while LMs, LMMs GLMs and GLMMs were fitted using the lme4 package (Bates et al., 2015) and model selection with information criteria was conducted using the dredge function of MuMIn package (Barton, 2019) in the R 3.0.3 R software (R Development Core Team, 2011).

Results

Mesocosm experiment

Survival and size of germlings

Germling survival and growth showed inter-population variability in response to thermal conditions, with a relevant role of emersion temperature, either alone or in combination with seawater temperature. When all the treatment combinations were analysed for germlings from the ria populations, survival was only affected by the population of origin (Table 1a; relative importance of population factor = 1). Many more germlings from OF than from IA survived (mean percentage of survival ± SE = 48.78 ± 4% and 5.70 ± 0.06% for OF and IA, respectively, n = 48). However, analysis of all populations (at fixed values of 35‰ salinity and 20°C seawater temperature, see Methods) revealed that the effect of emersion temperature depended on the population of origin (see Table 1b, population × emersion interaction). At the high emersion temperature, more germlings from the Lugo populations than from ria populations survived (Fig. 3a). Initial differences between populations were found in the number of germlings that settled on each disc, with significantly higher densities in the rias than in Lugo (mean initial number of germlings per disc ± SE: IA = 88.40 ± 7.27; OF = 59.16 ± 4.12; Pz = 15.69 ± 2.53; SP = 13.08 ± 2.0, n = 48 and 12 for rias and Lugo populations respectively; Tukey post-hoc test, IA > OF > Pz = SP, for analysis see Supplementary table S2). Survival patterns did not seem to be significantly affected by these initial differences in density. Quantile regressions relating survival proportions to the initial density of germlings showed negative relationships between these variables, but with low goodness of fit (R¹) values and/or non-significant slopes (90th quantile regression: slope = −0.073, p = 0.025, R¹ = 0.041; 70th quantile regression: slope = −0.050, p = 0.156, R¹ = 0.016).

Table 1. Selection of the random and fixed structures of GLMMs on survival of germlings at the end of the first mesocosm experiment from a) rias populations and b) all populations. For the random structure, table shows likelihood ratio tests comparing different nested models. For the fixed structure, the LogLikelihood (logLik), AICc and Δi are shown for both, the subset of candidate models and the full model, and Weights (ɷ_i) only for the candidate models. Selected models are marked in bold. P: Population, S: Salinity, T: Seawater temperature, E: Emersion. In the fixed component, the multiplications signs (×) indicate the inclusion of all the interactions of lower order and the implicated main factors

a) Ria populations
Random structure	df	Model comparison	logLik	L. Ratio test	p value
MR1: P × S × T × E	16		−306.02
MR2: P × S × T × E + (1|Tank)	17	1 vs. 2	−305.49	1.0683	< 0.001
MR3: P × S × T × E + (Population|Tank)	19	2 vs. 3	−301.86	7.2617	< 0.001
Fixed structure	df	logLik	AICc	Δi	ɷi
MF1: P	6	−306.804	626.552	0	0.390
MF2: P + S	7	−305.993	627.258	0.705	0.274
MF3: P + E	7	−306.474	628.220	1.668	0.169
MF4: P + T	7	−306.488	628.248	1.696	0.167
Full model: P × S × T × E	19	−301.855	651.711	25.158
b) All populations
Random structure	df	Model comparison	logLik	L. Ratio test	p value
MR1: P × E	8		−386.8
MR2: P × E + (1|Tank)	9	1 vs. 2	−282.7	208.19	< 0.001
Fixed structure	df	logLik	AICc	Δi	ɷi
MF1: P × E	9	−282.703	588.1	0	1

Figure 3. Survival of germlings and final size (average values per disc) in the first mesocosm experiment. a) Survival of germlings for each combination of population and emersion (n = 6); b) final average size of germlings from the ria populations, for each combination of seawater and emersion temperature (n = 20–23). Mean values ± SE are shown. Different lower-case letters above bars indicate significant differences between means based on a posteriori Tuckey test. Population abbreviations as in Fig. 1

The combined effect of seawater and emersion temperature affected the final size of germlings from the ria populations (Table 2a; relative importance of the seawater temperature × emersion interaction = 0.42). Germlings exposed to cold seawater and low emersion temperatures reached the largest sizes (Fig. 3b). Salinity or population of origin, on the other hand, appeared in a smaller proportion of the subset of candidate models, suggesting limited effects of these factors on the final size of germlings (Table 2a; relative importance of salinity and population, 0.31 and 0.18, respectively). One week after the start of the experiment, salinity affected the size of germlings from rias, which were smaller in the high salinity (35‰) treatment (mean size of germlings ± SE = 647.42 ± 35.79 µm and 444.70 ± 18.11 µm for low and high salinity treatments respectively, n = 33; for analysis see Supplementary table S2). However, analysis of all four populations under fixed conditions of seawater temperature and salinity revealed that temperature during emersion had a negative effect on the final size of survivors (Table 2b, mean final size of germlings ± SE = 385.69 ± 35.40 µm and 276.74 ± 28.27 µm for low and high emersion temperatures respectively, n = 21). No differences in responses were detected between populations. However, one week after the start of the experiment, germlings from the Lugo populations (Pz and SP) growing under the high emersion temperature were smaller than those from the ria populations (IA and OF; the best model included the interaction term population × emersion, see Supplementary table S3 and Fig. 2).

Table 2. Model selection of the random and fixed structures of LMs models on size of germlings at the end of the first mesocosm experiment from a) rias populations and b) all populations. For the random structure, table shows likelihood ratio test comparing different nested models. For the fixed structure, the LogLikelihood (logLik), AICc and Δi are shown for both, the subset of candidate models and the full model, and weights (ɷ_i), only for the candidate models. Selected models are marked in bold. Abbreviations as in Table 1. Null = model with no fixed structure. In the fixed component, the multiplications signs (×) indicate the inclusion of all the interactions of lower order and the implicated main factors

a) Ria populations
Random structure	df	Model comparison	logLik	L. Ratio test	p value
MR1: P × S × T × E	17		−88.064
MR2: P × S × T × E + (1|Tank)	18	1 vs. 2	−88.063	0.003	0.481
MR3: P × S × T × E + (Population|Tank)	20	2 vs. 3	−88.009	0.107	0.496
Fixed structure	df	logLik	AICc	Δi	ɷi
MF1: T × E	5	−60.141	131.031	0	0.172
MF2: E	3	−62.412	131.117	0.085	0.165
MF3: P + E	4	−61.913	132.319	1.287	0.090
MF4: S + T × E	6	−59.633	132.330	1.299	0.090
MF5: P + T × E	6	−59.649	132.361	1.330	0.089
MF6: E + T	4	−61.940	132.373	1.342	0.088
MF7: Null	2	−64.128	132.401	1.369	0.087
MF8: S + E	4	−62.023	132.539	1.508	0.081
MF9: S × E	5	−60.989	132.729	1.697	0.074
MF10: S × E + T × E	7	−58.791	133.000	1.990	0.064
Full model: P × S × T × E	17	−55.805	154.611	23.580
b) All populations
Random structure	df	Model comparison	logLik	L. Ratio test	p value
MR1: P × E	9		−29.237
MR2: P × E + (1|Tank)	10	1 vs. 2	−29.226	0.022	0.439
MR3: P × E + (Population|Tank)	19	2 vs. 3	−28.451	1.549	0.604
Fixed structure	df	logLik	AICc	Δi	ɷi
MF1: E	3	−18.936	44.503	0	1
Full model: P × E	9	−16.832	57.289	12.786

Survival and growth of juveniles

Temperature during emersion was relevant in determining juvenile performance, sometimes interacting with seawater temperature. Thus, juvenile survival decreased in the second week after the start of the experiment, especially for the high emersion temperature; this pattern became more evident over time (Fig. 4a). Ten days after the start of the experiment (22 September), the effect of emersion was dependent on seawater temperature (Fig. 4b and Table 3a; relative importance of seawater temperature × emersion interaction = 0.43). At this time period, survival of juveniles was lower for warm seawater and high emersion temperatures, regardless of the population of origin (Fig. 4b). However, 17 days after the start of the experiment (29 September), the effect of emersion appeared to depend on the population of origin (Table 3b; relative importance of population × emersion = 0.21). At the high emersion temperature, juvenile survival was much lower, irrespective of the population origin. By contrast, at low emersion temperature, survival of the juveniles from OF (ria population) was highest and survival of the juveniles from Pz (northern population) was lowest. Differences in survival between emersion treatments were also more evident in individuals from OF and less obvious in those from Pz (Fig. 4c).

Table 3. Selection of the random and fixed structures of GLMMs on survival of juveniles: a) 10 days (22 September) and b) 17 days (29 September) after the start of the mesocosm experiment. For the random structure table shows likelihood ratio tests comparing different nested models. For the fixed structure, the LogLikelihood (logLik), AICc and Δi are shown for both, the subset of candidate models and the full model, and weights (ɷ_i) only for the subset of candidate models. Null= model with no fixed structure. Selected models are marked in bold. Abbreviations as in Table 1. In the fixed component, the multiplications signs (×) indicate the inclusion of all the interactions of lower order and the implicated main factors

a) 10 days
Random structure	df	Model comparison	logLik	L. Ratio test	p value
MR1: P × S × T × E	32		−165.16
MR2: P × S × T × E + (1|Tank)	33	1 vs. 2	−146.43	37.459	< 0.001
Fixed structure	df	logLik	AICc	Δi	ɷi
MF1: T × E	6	−51.986	133.165	0	0.301
MF2: T + E	5	−49.513	133.311	0.146	0.280
MF3: Null	2	−53.834	134.381	1.216	0.164
MF4: S + T	5	−51.558	134.833	1.668	0.131
MF5: S + T × E	7	−49.018	134.936	1.771	0.124
Full model: P × S × T × E	33	−38.693	208.434	75.269
b) 17 days
Random structure	df	Model comparison	logLik	L. Ratio test	p value
MR1: P × S × T × E	32		−199.10
MR2: P × S × T × E + (1|Tank)	33	1 vs 2	−189.97	18.269	<0.001
Fixed structure	df	logLik	AICc	Δi	ɷi
MF1: P × E	7	−62.327	139.587	0	0.213
MF2: P + T + E	8	−61.504	140.217	0.630	0.156
MF3: E	4	−66.021	140.368	0.780	0.144
MF4: T + E	5	−65.233	140.958	1.371	0.107
MF5: P + T × E	9	−60.730	140.985	1.397	0.106
MF6: P + S × E + T × E	13	−55.967	141.128	1.540	0.099
MF7: P + S + E	8	−61.996	141.203	1.615	0.095
MF8: T × E	6	−64.434	141.562	1.975	0.079
Full model: P × S × T × E	33	−49.628	192.848	53.261

Fig. 4. Survival of juveniles in the first mesocosm experiment. a) Survival over time for each combination of population and emersion treatment (n = 16); b) survival of juveniles 10 days after the start of the experiment for each combination of seawater and emersion temperature (n = 16); c) survival of juveniles 17 days after the start of the experiment for each combination of population and emersion treatment (n = 4). In a) mean values and the smallest and largest SE bars are shown; in b) and c) mean values ± SE are shown. Different lower-case letters above bars indicate significant differences between means based on Tukey post hoc test. Population abbreviations as in Fig. 1

The growth of juveniles 10 days after the start of the experiment (22 September) was influenced by additive effects of population of origin and emersion temperature (Table 4a; relative importance of population = 1, emersion = 1). The weight of all juvenile survivors (particularly those from IA) decreased, with the exception of individuals from OF, which increased slightly in weight (Fig. 5a). Likewise, at high emersion temperature the weight of the juveniles decreased, while in the low temperature treatment the weight increased (Fig. 5b). The initial weight of juveniles differed among populations, with juveniles from IA being heavier than those from the other populations (mean initial weight ± SE, IA = 0.51 ± 0.010 g; OF = 0.35 ± 0.006 g; Pz = 0.39 ± 0.008 g; SP = 0.36 ± 0.008 g, n = 96, Tukey post-hoc test, IA > OF = Pz = SP; see Supplementary table S4 for model selection). Similar results were obtained by considering relative growth rates rather than absolute growth values (results not shown) despite the differences in initial sizes. The growth of juveniles 17 days after the start of the mesocosm experiment (29 September) was only analysed for low emersion conditions (see Methods). At that time, juvenile growth depended on the additive effects of population of origin and seawater temperature (Table 4b; relative importance of population: 1, temperature: 0.61), and growth tended to increase in juveniles from rias (particularly those from OF) and in juveniles exposed to the seawater temperature (Figs 5c–d).

Table 4. Selection of the random and the fixed structures of LMMs on growth of juveniles: a) 10 days (22 September) and b) 17 days (29 September) after the start of the mesocosm experiment. For the random structure table shows likelihood ratios test comparing different nested models. For the fixed structure, the LogLikelihood (logLik), AICc and Δi are shown for both, the subset of candidate models and the full model, and weights (ɷ_i) only for the candidate models. Selected models are marked in bold. Abbreviations as in Table 1. In the fixed component, the multiplications signs (×) indicate the inclusion of all the interactions of lower order and the implicated main factors

a) 10 days
Random structure	df	Model comparison	logLik	L. Ratio test	p value
MR1: P × S × T × E	33		450.06
MR2: P × S × T × E + (1|Tank)	34	1 vs. 2	457.99	15.876	<0.001
MR3: P × S × T × E + (Population|Tank)	43	2 vs. 3	466.36	16.736	0.003
Fixed structure	df	logLik	AICc	Δi	ɷi
MF1: P + E	19	613.933	−1187.123	0	0.367
MF2: S + P × E	20	615.001	−1186.959	0.163	0.339
MF3: T + P × E	20	614.233	−1185.422	1.701	0.157
MF4: S + T + P × E	21	615.255	−1185.150	1.973	0.137
Full model: P × S × T × E	43	623.584	−1146.211	40.911
b) 17 days
Random structure	df	Model comparison	logLik	L. Ratio test	p value
MR1: P × S × T	17		147.30
MR2: P × S × T + (1|Tank)	18	1 vs. 2	146.18	2.244	0.226
MR3: P × S × T + (Population|Tank)	27	2 vs. 3	143.18	5.986	0.524
Fixed structure	df	logLik	AICc	Δi	ɷi
MF1: P + T	9	210.042	−400.547	0	0.392
MF2: P	8	208.879	−400.538	0.008	0.390
MF3: P × T	12	213.059	−399.382	1.165	0.218
Full model: P × S × T	20	215.000	−382.074	18.473

Figure 5. Growth of juveniles in the first mesocosm experiment. a) Growth (difference in final and initial fresh weights) 10 days after the start of the experiment for each population (n = 67–79); b) growth 10 days after the start of the experiment for high and low emersion treatments (n = 140–157); c) growth 17 days after the start of the experiment for each population (n = 37–23); d) growth 17 days after the start of the experiment for cold and warm seawater temperatures (n = 63–64). Mean values ± SE are shown. Different lower-case letters above bars indicate significant differences between means based on a Tukey post hoc test. Population abbreviations as in Fig. 1

Post-stress experiment

The population of origin and the effect of previous treatments influenced the performance of recruits maintained under common conditions during the post-stress phase. Thus, survival of recruits from OF was much higher than that of recruits from IA and of those previously exposed to cold seawater temperature, i.e. 16°C, that also was the set temperature in the post-stress phase (Table 5a, mean percentages of survival ± SE, population effect = 39.16 ± 8.74% and 63.02 ± 8.13% for IA and OF, n = 16 and 15; seawater temperature effect = 63.16 ± 7.32% and 33.46 ± 9.32%, at cold and warm seawater temperatures, n = 18 and 13). Discs from OF presented higher initial recruit densities than those from IA (the selected model was the full model which included population, df = 4, AICc=505.95; mean initial number of recruits per disc ± SE, IA = 33.04 ± 8.92; OF = 40.8 ± 13.05; n = 16–15).

Table 5. Selection of the best linear model on survival and growth of recruits from rias populations in the post-stress experiment for: a) survival of recruits, b) initial size of recruits and c) final size of recruits. The table shows the LogLikelihood (logLik), AICc, Δi and AICc weights (ɷ_i). Selected models are marked in bold. Abbreviations as in Table 1

	df	logLik	AICc	Δi	ɷi
a) Survival of recruits
MF1: P + T	4	−95.156	199.852	0	0.715
MF2: P × T	5	−94.648	201.695	1.844	0.285
Full model: P × S × T	9	−93.793	214.157	14.305
b) Initial size of recruits
MF1: T	3	−110.524	227.936	0	0.318
MF2: Null	2	−111.908	228.245	0.309	0.273
MF3: P + T	4	−109.474	228.486	0.550	0.242
MF4: P	3	−111.164	229.218	1.281	0.167
Full model: P × T × S	9	−108.283	243.137	15.201
c) Final size of recruits
MF1: Null	2	−111.908	228.245	0	0.520
MF2: T	3	−111.338	229.566	1.320	0.269
MF3: P	3	−111.582	230.053	1.808	0.211
Full model: P × S × T	9	−108.430	243.432	15.186

At the start of the post-stress experiment, the recruits from ria populations previously exposed to cold seawater temperature were larger than those exposed to warm seawater (Table 5b, mean initial size ± SE = 588.17 ± 57.39 µm and 512.33 ± 84.11 µm at cold and warm seawater temperatures respectively, n = 18 and 13). However, at the end of the experiment, the size of recruits was not affected by prior exposure to different seawater temperature or salinity or by the population of origin (Table 5c).

Discussion

The study findings revealed inter-population variability in the resistance to heat stress at the southern rear-edge of the distribution range of Fucus serratus. A growing body of literature shows greater thermal tolerance in rear-edge than in central populations, particularly in several habitat-forming algae including fucoids and kelps (e.g. Pearson et al. 2009; Ferreira et al., 2014; Bennett et al., 2015; Saada et al., 2016; King et al., 2018). However, edge-edge comparisons as shown in this study are scarce, despite their importance in unravelling range dynamics under present and future climate change scenarios. Our results also highlighted the strong influence of aerial thermal stress, either alone or in synergy with seawater temperature, on the performance of early stages of F. serratus, and thus its potential role in establishing the southern boundary of the distribution. Finally, the effect of salinity was subtle, initially only affecting the size of germlings in the ria populations.

Germlings from the northern semi-exposed populations on the Lugo coast were more resistant to thermal stress during emersion (greater survival) than those from western rias on the Atlantic coast. In juvenile stages, the effect of emersion temperature was predominant, without any evident inter-population differences between the northern and the ria populations in vital rates. However, the greater resistance of Lugo populations to heat stress was also detected in experimental studies with adult plants. The survival threshold of adult plants to seawater temperature was 3°C higher in one population from Lugo coast than the thresholds of other populations along the species distributional range, including populations of western rias of Iberian Peninsula (A. García, unpublished data).

The higher thermal tolerance of germling and adult F. serratus occurs under current harsh conditions on open northern shores. Habitat quality recently declined in the western part of northern Spanish coasts, as shown by inhibition of the growth of juvenile F. serratus transplanted from rias to the Lugo coast, and the decrease in plant size, reproduction and recruitment of the extant populations (Duarte & Viejo, 2018). Significant warming and limited nutrient supply in summer due to increased water stratification have been detected on the northern coast of Spain in the last few decades (Llope et al. 2006, 2007; Gómez-Gesteria et al., 2008). Smaller germlings and reduced growth of juveniles from northern populations were also recorded in our mesocosm experiments. This may reflect physiological costs associated with stress tolerance. Metabolic costs are involved in the synthesis of antioxidants, phenolic compounds or heat shock proteins, as part of the ubiquitous cellular responses to heat stress (Feder & Hofmann, 1999; Contreras-Porcia et al., 2017). High levels of expression of constitutive heat shock proteins have been detected in Spanish populations of F. serratus in relation to core populations of N Europe (Jueterbock et al., 2014).

The underlying mechanisms of the greater thermal tolerance in northern Spanish populations than in western rias may be due to local adaptation and/or acclimatization (Kawecki, 2008). Experiments often do not enable genetic responses to be disentangled from plastic responses, which may also be transgenerational (Munday et al., 2013). Our findings indicated that germlings/recruits may rapidly acclimate to specific environmental conditions, as survival was higher during the post-stress phase in those specimens previously exposed in the mesocosm experiment to the same temperature. Southern populations of F. serratus are genetically unique and exhibit pronounced inter-population differences relative to core populations in central and northern Europe, although with lower genetic diversity (Coyer et al., 2003). The limited dispersal capacity of the species and the high genetic differentiation at a local scale may favour ecotypic differentiation at the southern edge of distribution (Coyer et al., 2003).

Whether adapted or acclimatized to the present conditions, the greater resistance of germlings and adult plants to thermal stress may offset poor recruitment, favouring persistence of the Lugo populations. Nevertheless, the long-term viability of these peripheral northern populations is uncertain and will depend on the rate of future climate change and the adaptive potential of these populations to face the forthcoming conditions at a fast enough speed to avoid local extinction. Selection may be strong under increasing environmental stress, as observed at small scales on the steep intertidal gradient (e.g. Hays, 2007). Spanish populations of F. serratus exhibit signs of increased selection pressure for heat tolerance over time (Jueterbock et al., 2018). However, the lower genetic diversity detected in the southern populations of this species, in relation to core locations in northern European shores, may limit their evolvability (Coyer et al. 2003; Pearson et al., 2009). Future climate change scenarios may exceed the potential for edge populations from northern Spain to persist. In fact, north-eastern populations located in the province of Asturias, which were abundant in the 1990s and constituted then the boundary of F. serratus distribution, are currently almost extinct (Arrontes, 2002; Duarte et al., 2013).

While conditions are becoming harsher for F. serratus and other habitat-forming algae in northern Spain, such as F. vesiculosus, Himanthalia elongata or kelps (Nicastro et al., 2013; Casado-Amezúa et al., 2019 and references therein), the western rias still provide benign conditions for several macroalgae at the southern edge of their distributions (Duarte & Viejo, 2018). Rias can act as contemporary refugia, providing favourable environmental conditions related to the presence of intense summer upwelling events, shelter from wave-action and inputs from rivers (Fraga, 1981; Álvarez et al., 2005). The role of upwelling areas as refugia for ecosystem engineers, including canopy-forming fucoids, has previously been highlighted (Chollett et al., 2010; Lourenço et al., 2016).

Higher heterogeneity in the performance of early stages between ria populations than between northern populations was also detected in the mesocosm experiments, with more survivors and faster growth rates in the O Freixo population than in the Isla de Arousa population. These sites are located in the inner part of two rias characterized by different flow rates and morphological and topographical conditions (Rosón et al., 1995; Carballo et al., 2009). Cooler waters in summer and higher inorganic nutrient supply were recorded in O Freixo than in the Isla de Arousa (Duarte & Viejo, 2018). The geographic isolation and higher environmental heterogeneity may promote larger phenotypic and genetic divergence between these two populations than between populations from northern open shores.

The study findings also revealed the important role played by thermal stress during emersion, either alone or in synergy with seawater temperature, on the vital rates of early stages of F. serratus from edge populations. Although variable, the responses of several species of genus Fucus to combined abiotic stresses were mostly additive (Wahl et al., 2011). In particular, F. serratus adults were affected additively by solar radiation, seawater and air temperature (Martínez et al. 2012b; Fernández et al., 2015). Nevertheless, our contrasting results are not surprising, as responses often differ among life stages, with early stages being more vulnerable to stressors than adults, particularly to aerial exposure (Vadas et al., 1992; Davison & Pearson, 1996; Nielsen et al., 2014). Germling size and juvenile survival were both synergistically affected by air and seawater thermal stress. When multiple stressors occur simultaneously or in rapid succession, synergistic effects are likely to occur, because the second stressor will increase the intensity or duration of the first stressor (Gunderson et al., 2016). Repeated exposure to high air and water temperatures may lead to higher physiological costs.

In our experiments, the synergies occurred in moist conditions during emersion, under adult canopies, simulating field conditions. A previous study with adult fronds of F. serratus indicated that humidity ameliorated aerial thermal stress (Martínez et al., 2012b; Fernández et al., 2015). However, desiccation may confer protection during aerial exposure and increase thermo-tolerance in intertidal algae by maintaining fronds in an inactive state (Davison & Pearson, 1996; Mota et al., 2015). Further work is necessary to elucidate whether moisture exacerbates or alleviates thermal stress during exposure of early developmental stages.

The southern limits of the geographic distribution of various species, particularly of intertidal seaweeds, have frequently been attributed to the major influence of seawater temperature (e.g. Lüning, 1990; Eggert, 2012; Jueterbock et al., 2013; Saada et al., 2016). Our results indicate that thermal stress during aerial exposure, and the synergistic effects of air/water warming, probably play a key role in determining the southern limits of distribution of F. serratus. Previous field observations also pointed out the importance of stress at low tide in the decline and subsequent extinction of eastern populations of F. serratus in northern Spain (Viejo et al., 2011; Martínez et al., 2012b). Indeed, molecular thermal damage of intertidal organisms in temperate regions occurs at temperatures almost exclusively experienced during emersion periods (Helmuth et al., 2002 and references therein).

In summary, populations from rias and northern shores in NW Iberian Peninsula exhibited distinct resistance to heat stress and experienced contrasting local conditions. The future of these two distinct groups of marginal populations remains an open question and depends on trends in climate change and the ability of populations to cope with these. Environmental changes may already be occurring at rates that exceed the plastic and adaptive potential of edge populations in northern Spain, despite the detected greater resistance to heat stress of these populations than those from rias. Global warming and specifically heat waves are increasing in frequency and intensity (Smale et al., 2019). Warming is highly heterogeneous in coastal areas and is much less pronounced in upwelling regions (Varela et al., 2018). However, several studies have highlighted the weak upwelling events both in N and W Iberia during the last few decades (Llope et al., 2006; Pérez et al., 2010). Furthermore, a recent study using high-spatial resolution climate models predicts the future weakening of the NW Iberia upwelling (Sousa et al., 2020) under emissions scenario RCP 8.5, which is characterized by large increases in greenhouse gas emissions (Riahi et al., 2011). The future of western rias as climate refugia for F. serratus and several foundation species will thus be threatened, and the stressful conditions on the northern coast will also be exacerbated. Whether these heterogeneous rear-edge NW Iberian populations persist or become extinct is a key question with important ecological and evolutionary implications. The disappearance of habitat-forming species would probably cause bottom-up cascading effects, leading to functional simplification and a decrease in the productivity potential of coastal systems (Duarte et al., 2015; Filbee-Dexter & Wernberg, 2018). Furthermore, F. serratus and other fucoids and kelps harbour unique genetic and phenotypic variation at the warmer rear edge, in NW Iberia (Coyer et al., 2003; Provan, 2013; Neiva et al., 2015 and references therein). The loss of this unique and heterogeneous phenotypic and genetic component at the rear edge of the distribution may compromise the adaptive potential of these key species as a whole to increasing warming conditions. Global actions that slow the future rate of greenhouse gas emissions could reduce the impacts of climate change, by favouring the potential of marine species to acclimate or adapt to the new conditions.

Supplementary figure S1. Diagram of the mesocosm experimental design showing the different treatment combinations and populations per tank. S: Salinity, T: Seawater temperature, E: Emersion (HE: high emersion temperature; LE: low emersion temperature; S22 and S35: treatments with salinity of 22 and 35‰; T16 and T20: treatments with seawater temperatures of 16 and 20°C). See Materials and methods for details. For juveniles, 4 tanks and 12 individuals (3 for each population) per tank were included in each S × T × E combination. For germlings, 3 tanks and 8 (or 4) discs were used per tank (2 discs per population). IA, Isla de Arousa; OF, O Freixo; Pz, Peizás; SP, San Pedro.

Supplementary figure S2. Size of germlings one week after the start of the mesocosm experiment for each combination of Population and Emersion (n = 36). Mean values ± SE are shown. Different lower-case letters above bars indicate significant differences between means based in a posteriori Tukey test. Population abbreviations, as in Appendix 1.

Supplementary table S1. Selection of random and fixed structure in the GLMM for the initial number of germlings. For the random structure the LogLikelihood (logLik) and the results of the likelihood ratio test comparing nested models are shown. For the fixed structure, logLik, AICc, Δi and AICc weight (ɷ_i) are shown for the candidate (and full) model. Selected models are marked in bold. P: Population; E: Emersion

Supplementary table S2. Selection of the random and fixed structure in LMMs for size of germlings one week after the start of the experiment for the ria populations. For the random structure, the LogLikelihood (logLik) and the results of the likelihood ratio tests comparing different nested models are shown. For the fixed structure, logLik, AICc and Δi are shown for both the subset of candidate models and the full model, and the weights (ɷ_i) are only shown for the subset of candidate models. Selected models are marked in bold. Null: model with no fixed components. P: Population; S: Salinity; T: Temperature; E: Emersion.

Supplementary table S3. Selection of random and fixed structure in LMMs for size of germlings from all populations one week after the start of the mesocosm experiment. For the random structure, the LogLikelihood (logLik) and the results of the likelihood ratio tests comparing different nested models are shown. For the fixed structure, logLik, AICc, Δi and weight (ɷ_i) are shown for the candidate (and full) model. Selected models are marked in bold. P: Population; E: Emersion.

Supplementary table S4. Selection of the random and fixed structure in LMMs for initial weight of juveniles. For the random structure, the LogLikelihood (logLik) and the results of the likelihood ratio tests comparing different nested models are shown. For the fixed structure, logLik, AICc, Δi and AICc weight (ɷi) are shown for the candidate (and full) model. Selected models are marked in bold. P: Population.

Author contributions

A.G. Garcia: fieldwork, experimental set-up, execution of the experiment, statistical analysis, writing and editing the manuscript; C. Olabarria: fieldwork, experimental set-up, execution of the experiment, writing and editing the manuscript; O. Álvarez-Losada: fieldwork, experimental set-up, execution of the experiment, statistical analysis and editing the manuscript. R.M. Viejo: fieldwork, experimental set-up, execution of the experiment, statistical analysis, writing and editing the manuscript.

Acknowledgements

We are grateful to Christine Francis for improving the English. We wish to thank ECIMAT for providing the facilities to develop the experiments and allowing us to use their equipment. Constructive comments of editors and two anonymous reviewers greatly improved the manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary information

The following supplementary material is accessible via the Supplementary Content tab on the article’s online page at https://doi.org/10.1080/09670262.2020.1842909

Additional information

Funding

This research was supported by project CGL2014-60193-P of the Spanish Ministry of Economy and Competitiveness.