Four populations of E. singularis were sampled simultaneously in October 2008 at two depths, hereafter named either shallow (–10 to –15 m) or deep (–35 to –40 m), and at either the Riou Island (43°10′37′′N05°23′58′′E, Marseille, France) or the Medes Islands (42°02′49′′N03°13′21′′E, Estartit, Spain). Apical colony tips (10 cm long) were collected randomly by scuba divers from different healthy adults (colony height > 30 cm) in each population. All experiments were designed to limit the numbers and sizes of samples removed from the ecosystem. According to the geographic origin, –35 m does not correspond to the deepest distribution limit [27,28], but this however corresponds to a bathymetry at which specimens are living in contrasting thermal regimes by comparison to shallow ones [10]. At both geographic sites, samples were collected over an area of around 500 m². Gorgonian tips were glued onto holders with an Epoxy stick (HoldFast, Aquarium Systems) (Fig. 1) and equally distributed between two aquarium sets: one for the temperature treatments (see below) and the other for the control. Each set was composed of four 35-L tanks: two for the shallow-water populations and two for the deep water ones. Each tank was set up with tips from the Riou and the Medes Islands (see below for the number of replicates). All the tanks were supplied with oligotrophic Mediterranean seawater, which was changed weekly and homogenized thanks to submersible pump (Aquaball, EHEIM, Germany). Irradiance was set at 250 and 70 μmol photons m⁻² s⁻¹ for the shallow and deep populations respectively, to match the light intensity that the populations received in the field as closely as possible, with light/dark cycles of 12 h/12 h. Seawater was heated using submersible resistance heaters (Visi-Therm^® Deluxe) and the temperature recorded using Testo T175-T3 data loggers, at a frequency of four measurements per hour. All thermal indications in the text correspond to mean recorded values ± 0.3 °C. Salinity values were kept constant at 38‰ during all experiments. Specimens were fed twice a week with extracts of frozen nauplii of Artemia salina (4 g total quantity per tank).

Two heat stress experiments were carried out, hereafter named: HTD for “Heat Threshold Determination” experiment and TTN for “Time To Necrosis” experiment. For the HTD experiment, after a 2-week acclimation period at 20 °C, the temperature of the experimental set was gradually increased by 2 °C every week up to 30 °C (while control tanks were maintained at 20 °C) (Fig. 2). Each experimental tank was composed of eight colonies from Riou and eight from Medes Islands, which corresponds to n = 64 for control and n = 64 for treated samples. For the TTN experiment, the temperature was regularly increased by 2 °C every week from 20 °C to 26 °C, followed by a final step to 27 °C, which was maintained until total necrosis of all gorgonian tips had occurred (while the control tanks were maintained at 20 °C) (Fig. 3). Each experimental tank was composed of four colonies from Riou and four from Medes Islands, which corresponds to n = 32 for control and n = 32 for treated samples.

The morphometric analysis was conducted according to [29] (Supplementary data, Fig. S1). In the present work, gorgonian tip photosampling was performed by taking pictures of gorgonian fragments (Fig. 1) at exactly the same position and focal distance using a digital camera (Konica-Minolta Dimage G600). To obtain the best contrast between the healthy white tissues (when polyps are retracted) and the brown central axis (shown after tissue necrosis), a red screen was placed just behind the gorgonian fragments and a black square of 1 cm² was drawn on this frame to serve as an internal calibration area (Fig. 1). At the end of the photosampling, the pictures were analyzed using a custom macro program running under Image-J software. This macro detects semi-automatically and quantifies the area of healthy tissue from the transversal view of each gorgonian tip. At each time-point, the quantity of living tissue was expressed in cm² of planar white area. The source code of the Image-J macro, which can calculate the surface of various contrasting color specimens, and a complete set of instructions are available from the corresponding author on request. At the beginning of the TTN experiment, the transversal diameter of each colony was measured. After completing the morphometry analysis of the TTN experiment, the heat stress duration between the day when the first specimen reached 50% necrosis to the end of the experiment was subdivided into three groups of equal length, hereafter named sensitivity classes (Fig. 3). This allowed the gorgonian colonies to be qualitatively classed in sensitive specimens (noted hereafter “sensitive”), intermediate specimens (“intermediate”) and resistant ones (“resistant”) with respect to this specific stress. More precisely, “sensitive” class contains specimens that presented 50% necrosis between the 42nd and the 57th day of experimentation, “intermediate” class reached this level between the 57th and the 69th day, “resistant” class between the 69th and the 82nd day (Fig. 3). The final photosampling sizes of the HTD and the TTN experiments were respectively 23 and 21 pictures (Figs. 2 and 3).

Final morphometric results are presented as means ± sem percentage of necrosis severity, obtained by dividing the surface of remaining living tissue by the initial tissue surface at each time-point. Since the calculation of the percentage of necrosis severity was repeated during the experiment for each individual, multivariate repeated measures General Linear Model ANOVA (GLM) was performed to detect the main and interactive effects of geographic and bathymetric origins respectively to the experimental treatment (STATISTICA 9, StatSoft Inc. 2009). Pairwise post-hoc comparisons were carried out using a Tukey’ post-hoc test. One-way ANOVA was used to determine if the initial diameter of gorgons depends on their geographical or bathymetry origins. Then, contingency test was carried out to determine statistic links between the initial diameter of specimens and their distribution within the sensitivity classes.

Genetic analyses were carried out at the beginning of both experiments by splitting a 1 cm long fragment from each of the sampled colonies. These fragments were conserved in 70% ethanol at –80 °C until DNA extraction, which was performed following a standard “salting out” procedure. Briefly, each fragment was incubated in TNES buffer (0.1 mM Tris, pH 7.5; 0.5 M NaCl; 50 mM EDTA; 0.6% SDS) plus 28 μg·ml⁻¹ proteinase K (Roche), for a minimum of 2 h at 56 °C and with gentle shaking. Proteins were precipitated with 5 M NaCl and removed by centrifugation at 13,000× g for 5 min at room temperature. The supernatant fraction was collected and DNA was precipitated overnight in 96% ethanol at –20 °C. After centrifugation at 13,000× g for 5 min, the DNA pellet was washed with cold 70% ethanol. The dried DNA pellet was then dissolved in TE buffer (1 mM TRISMA BAZE, pH 8.0; 1 mM EDTA) and stored at –80 °C until genetic analysis. All gorgonian tips were genotyped using 10 microsatellite loci, five from the host (named: C21, C30, C40, S10, and S14) and five from the symbiont (named: BVB2, CACV19, CACV30, CAES67, and S19) [30]. An additional table, available as supplementary data, summarizes the genetic diversity of E. singularis within host and symbiont samples. The study about the development of microsatellite loci from E. singularis and its symbiont is available under the following link: http://tomato.biol.trinity.edu/manuscripts/10–1/mer-09–0225.pdf. The PCR amplifications were carried out either locus by locus, or in multiplex using the QIAGEN^® Multiplex PCR Kit according to Molecular ecology resources primer development consortium et al. [30]. Genotype screenings were made on an ABI 3130 Genetic Analyser 53 (Applied Biosystems) and analyzed using STRAND software.

To explore the genetic differentiation between populations, we used the genetic differentiation index Jost D [31]. Host Jost D values were calculated using SPADE [32]. For the symbiont genetic differentiation analyses, since the genotyping method does not allow within-colony allele frequencies to be determined, we simply considered each allele as present or absent from each colony. Population allele frequencies over colonies were calculated on this basis, and from these values, we could calculate the symbiont Jost D values for each locus. Symbiont Jost D values over loci were computed as the mean of the values obtained from each locus. As SPADE only computes Jost D significance for a single locus (here from 500 permutations), significant departures from 0 over all loci were tested using Multitest V.1.2 [33], which allows multiple testing of the same hypothesis, without knowing the exact P value of each partial test. In our case, every test on a particular locus was considered as an independent test of the same null hypothesis, i.e., absence of genetic differentiation between populations.

To evaluate a putative linkage between the thermotolerance and the genotype, host and symbiont Jost D values were finally computed for the two geographical origins of the specimens (Medes and Riou Islands) and depending on the three sensitive classes (sensitive, intermediate and resistant).

In the control set of tips maintained at 20 °C, neither signs of necrosis nor tissue thinning were observed in the colonies, whatever their geographical or depth origin (Fig. 2). The statistical tests further confirmed this, and revealed a very significant effect of temperature treatment (GLM repeated measures ANOVA; F19, 2059 = 305.93, P = 0.0000). The Tukey post-hoc test showed first signs of significant tissue thinning of heat stressed samples at the 29th day, corresponding to the 26 °C step, for all four populations (Tukey post-hoc test, P < 0.001) (Fig. 2). The first signs of severe necrosis, when the loss of tissue reveals the underlying brown axial skeleton, were thereafter measured on the 36th day, corresponding to the 28 °C step, at all sites and depths (Fig. 2). Total loss of tissue was finally observed from the 47th day for all specimens, corresponding to the 30 °C step. Consequently, 28 °C appeared to be the upper thermal threshold of rapid tissue necrosis in all the samples, collected either at Medes or Riou Islands and at either –15 or –30 m depth.

The TTN experiment was then carried out in order to highlight differences in thermotolerance capacity of E. singularis populations submitted to a sustained thermal stress of 27 °C (Fig. 3). In the control set of tips maintained at 20 °C, neither sign of necrosis nor tissue thinning were observed in the colonies, whatever their geographical or depth origin. The statistical tests further confirmed this statement and revealed a very significant effect of temperature treatment (GLM repeated measures ANOVA; F19, 1064 = 93.80, P = 0.0000). The first signs of significant tissue thinning of heat stressed samples appeared the 29th day at 26 °C for the four populations (Tukey post-hoc test, P < 0.001). The statistical tests also shown no statistical bathymetry effect, but a clear geographic effect (GLM repeated measures ANOVA; F19, 532 = 2.57, P = 0.0003). However, we could not further analyze this geographic effect, since assumptions of normality and homoscedasticity were not met at several time-points. Aiming to better analyze this dataset, we therefore considered individual sample distribution in three classes of thermal sensitivity (Fig. 3 Material and methods for details). Indeed, the overall analyze of the necrosis evolution for all samples displayed high within-population variability. For instance, the threshold of 50% necrosis severity was reached on the 42nd day for some specimens and only on the 78th days for others. Both from direct observations in aquarium, and from necrosis severity measurements, we clearly observed three groups of individuals presenting contrasting necrosis development. A first group of gorgonian tips presented early and rapid tissue degradation and were therefore placed in the sensitive class. On the other hand, some samples presented a very delayed onset of necrosis, and therefore constituted the resistant class. Between these two groups, other colonies displayed more progressive tissue damages, representing the intermediate sensitivity class. When considering the 50% necrosis threshold, a 5-day time lap was observed between the last individual of one class and the earliest of the other. Thus, no individual was at the border of two consecutive sensitivity classes, confirming the reliability of the class distribution. For the 50% necrosis severity threshold, the composition of these sensitivity classes revealed that within each geographical origin there was no depth-effect. However, Riou and Medes populations did not show the same trend of resistance (contingency tests, P < 0.05), and this observation holds true regardless the threshold considered: 50%, 75% or 90% necrosis severity. Indeed, when gathering shallow and deep populations, a first half of the Medes specimens suffered early necrosis, while the other half only displayed severe necrosis at the end of the TTN experiment (Fig. 4). In contrast, shallow and deep Riou populations mostly displayed necrosis in the intermediate sensitivity class (Fig. 4). All these statements explain why the mean necrosis evolution of the Medes samples fits a double-sigmoid function (Fig. 3), since it integrates both sensitive and resistant specimens. This also explains why at some time-points, assumptions of normality and homoscedasticity of these values were not met, which justifies the use of the sensitivity class approach.

Sample distribution among the three thermal sensitivity classes was further analyzed considering the initial diameter of the gorgonian samples. The diameter measurements showed that colonies from Riou Island were initially significantly larger than those of the Medes Islands, at 2.80 ± 0.013 mm and 2.33 ± 0.007 mm, respectively (ANOVA; F1, 09 = 14.85, P = 0.0002). However, no significant difference in size was detected between populations from shallow and deep sites within each location. Diameter was tested for correlation with sensitivity classes using contingency tests, but no correlation was found (results not shown).

Jost D genetic differentiation indexes were computed for both host and symbiont, based on all the specimens used in the TTN experiment (n = 128). Considering the bathymetric distribution, no significant difference was found between populations from shallow and deep environments at either site for the host or the symbiont (Table 1). However, considering the geographic origin of the host, all Jost D values between populations of different sites were significant, with values ranging from 0.070 (Riou Island deep vs Medes Islands shallow) to 0.123 (Riou Island shallow vs Medes Islands deep). Host populations from different depths are therefore genetically homogenous within a geographic site, although a genetic differentiation was measured between the Riou and the Medes Islands. The same pattern was found for the symbiont (Table 1). Values of zooxanthellae genetic differentiation between geographic sites were three times higher than between the different depths (except for Riou Island shallow vs Medes Islands deep). In conclusion, neither host, nor symbiont was differentiated according to depth, although they were genetically different between the two sites.

Jost D genetic differentiation indexes were also computed for specimens belonging to the sensitive classes. Concerning the symbionts, no significant genetic differentiation was observed between the sensitive classes whatever the geographic origin (Table 2, upper right). Concerning the host (Table 2, lower left), since specimens from the Medes Islands are mostly distributed in both sensitive and resistant classes, this could suggest two different subpopulations. However, further Jost D analysis indicated that the sensitive specimens were not genetically different from the resistant ones. No statistical difference was observed among the three classes of Riou Island tips. However, a significant difference was revealed between the intermediate and resistant classes of Medes Islands’ samples. These results revealed several correlations between genotypic and heat resistance differences (Table 2, lower left). For example, individuals within the sensitive class from Medes islands were genotypically different to those within the intermediate and resistant classes from the Riou Island. Moreover, individuals within the resistant class from Medes islands were genotypically different to those within the intermediate and resistant classes from the Riou Island.

In the temperate Mediterranean Sea, experiments on gorgonian thermal threshold are still scarce [8,27,34] or focused on the non-lethal temperature stress response [25,26,35]. This study is the first to investigate the lethal thermal threshold of the symbiotic gorgonian E. singularis in experiments with a controlled gradual increase of seawater temperature. In the HTD experiment, we showed that the same temperature threshold induced abrupt tissue necrosis in all populations, whatever their depths or geographical origin.

Neither necrosis nor tissue thinning were observed for seawater temperatures below 26 °C. This observation is in agreement with previous data obtained from E. singularis hyperthermia experiments [25,26,35]. At 26 °C, we observed a thinning of tissue, but without denudation of the axial skeleton for all populations. Accordingly, previous studies [25,26] have shown that a seawater temperature of 26 °C, applied for more than 1 week, induced significant physiological and molecular modifications, but without apparition of complete tissue necrosis. In the present work, coenenchyme loss leading to partial axial denudation was only observed at 28 °C, followed by rapid and complete tissue destruction at 30 °C. It is clear that the length of the applied heat stress is a crucial parameter defining any lethal temperature threshold. Moreover, in conservation biology context, critical heat thresholds should be considered in conjunction with other many environmental factors, such as light intensity (including UV), water salinity, food availability (as further discussed), water turbidity, seasonality and the possible presence of heat-activated pathogens [13,14]. In any case, the present observations, made in controlled conditions and for a given stress length, which is perfectly corroborate the Garrabou et al. study [6] on the 2003 mass mortality events. Indeed, field investigations of several gorgonian species of the NW Mediterranean Sea revealed that for E. singularis the degree of necrosis was:

• • maximal (> 85% necrosis) when the summer temperature maxima (Tmax2003) were around 28 °C;
• • intermediate (> 15% necrosis) for a Tmax2003 of 26.5 °C;
• • minimal (< 5% necrosis) for a Tmax2003 of 25.4 °C [6].

Accordingly, the higher thermal limit driving E. singularis from a heat stress response to lethality is between 26 and 28 °C whatever the site of origin or bathymetry.

On the basis of the HTD experiment results, the TTN experiment was then performed with a gradual increase of seawater temperature from 20 °C to 27 °C, which was further kept constant to observe the gradual evolution of gorgonian necrosis. This heat stress experiment, whose design corresponds to a classical test of resistance, allowed us to discriminate differential necrosis responses, and also highlighted strong intra-population variability in both Riou and Medes populations (Figs. 3 and 4). This variability, observed using random sampling of wild individual colonies, is in agreement with the diversity classically observed in the field [6,8,12]. The distribution of samples among distinct groups accordingly to their necrosis evolution drove us to better analyze the TTN dataset by splitting individual heat-responses in three sensitivity classes. Thanks to this approach, we highlighted a statistical difference according to the geographic origin but not according to the bathymetry.

Deep and shallow specimens were not statistically differentially distributed either using the GLM repeated measures ANOVA or the contingency test on the sensitivity class distribution. These observations are peculiar when considering the original living temperature distribution with respect to the bathymetry. Indeed, mean thermal records (2004 to 2006) confirm a 3 °C decrease between –10 and –35 m depths at both sites [10]. Accordingly, the sustained step at 27 °C can be considered as 3 °C more intense for specimens from deeper environments than for specimens from shallower ones. Nevertheless, specimens from deeper origin did not undergo necrosis at earlier stage or at an accelerated rate in comparison to those from shallower ones. This result, which may seem “illogical”, is in fact reminiscent to our previous findings showing that E. singularis fragments collected at –35 m depth were more resistant than specimens from shallower depths when submitted to a steady stress at 26 °C [25,26]. In these previous studies, we proposed that shallow E. singularis populations are living close to their current thermal maxima and that conversely, the deeper populations display a better physiological plasticity to hyperthermic condition. Remarkably, at some locations Gori et al. [27,28] observed high densities of E. singularis populations until 70 m depth. Since, in the present study, the depth-dependant thermotolerance of the white gorgonian was only investigated on populations living between –15 and –35 m, it would be very interesting to further consider the thermal sensitivity of specimens living deeper, at the lowest limit of distribution.

The geographic origin effect was first detected using the GLM repeated measures ANOVA, and further analyzed with the contingency test on the sensitivity class distribution. This revealed that Medes Island specimens presented a bimodal distribution pattern while Riou samples were mostly distributed in the intermediate class (Fig. 4).

There could be several explanations for the differences in sensitivity class distribution according to geographic origin. The two sites were indeed chosen because they have contrasting thermal regimes [8,10], which could therefore have driven the gorgonian to develop contrasting thermal resistance capacities. Mean temperatures recorded from 2004 to 2006 indicate that Medes Islands waters are about 2 °C warmer than Riou Island waters at both depths [8,10]. Furthermore, the mean associated variance of seawater temperature is very different between sites. Indeed, the Riou Island is associated with high temperature variability at all depths, linked to wind regimes that are responsible for frequent upwelling of cold deep water. These water movements also affect the Medes Islands, where they instead induce downwelling of warm surface water. Accordingly, at the Medes Islands, shallow temperatures are fairly stable (almost three times less variable than at Riou), while deep temperatures are very variable (as much as at Riou). Thus, field observations of thermal conditions could explain, at least in part, the observed contrasting distribution of the colonies among the sensitivity classes. Another important point of comparison between these sites is the difference in impact of the summer mass mortality events of 1999 and 2003 (for the shallow populations). Despite higher overall temperatures, several field observations reported that Medes gorgonian populations were less impacted than Riou ones [5,8,10]. This implies that the past heat waves have exerted different levels of selective pressures. Of course, retrospectively, it is not possible to analyze the consequences that these past selective pressures could have had on unknown previous sensitivity class distributions. However, the present Riou population seems to be more homogenous in its response to stress, which could either result from the highly variable thermal regime, or from selection imposed by past mortality, or both.

Since we observed coenenchyme thinning prior to axial denudation, we suspected that heat resistance capacities could be directly correlated with initial sample diameter, which could be a proxy of the initial biomass. This would have introduced a bias in the subsequent analyses, since Riou samples were found to be statistically larger than Medes ones, suggesting potentially higher metabolic reserves. However, no correlation was observed between initial diameter and sensitivity class distribution, suggesting that energetic stores are not the main factor affecting differential thermal resistance capacities in the TTN experiment. In this, another important parameter to consider would be the role of the food availability in the thermal resistance. Such factor was previously observed for the temperate non-symbiotic gorgonian Paramuricea clavata, since a moderate heat stress (+3 °C above ambient temperatures) exerted less damaging effects when applied under high-food conditions, than in ambient-food ones [29]. These authors also proposed the hypothesis that the energetic constraint during the summer dormancy is the main mechanism triggering the heat wave mortality events [36,37]. It is important to remember that, in our experimental design, all symbiotic E. singularis tips were both fed and exposed to light/dark cycles at levels ensuring good photosynthetic activities of the zooxanthellae [25]. In conclusion, complementary experiments using different feeding regimes are needed to precisely define the respective contribution of heterotrophic and autotrophic metabolisms to the thermal sensitivity. A very recent work confirms that the nutritional status and the heat stress resistance are not directly correlated for this symbiotic sea fan [38].

Since both host and symbiont genotypes could contribute to the thermal resistance of the cnidarian holobiont [15,16], we genotyped both partners of the E. singularis symbiosis. Previous work has confirmed there is no within-colony genetic variability for both the host and the symbiont (Forcioli, pers. comm.). Consequently, the genetic differences analyzed here are true differences among individual colonies. In our work, we showed in Table 1, a lack of depth-dependent thermotolerance associated with a lack of depth differentiation for the symbiotic algae as previously demonstrated by Forcioli et al. [24]. The absence of genetic differentiation according to bathymetry is an interesting result since the analyzed populations are living under contrasting environmental conditions [10] and since E. singularis is suspected to present a low dispersal potential [39]. A recent study has been carried out on the genetic structure along a depth gradient (from 20 to 70 m) of Corallium rubrum, another Mediterranean gorgonian [40]. The authors have pointed out both a reduction of the genetic variability along the depth gradient and the presence of a threshold between –40 and –50 m depth in the connectivity among the samples, supporting the hypothesis of a local adaptation. Accordingly, in order to assess the depth-related genetic variability of E. singularis, it would be interesting to consider locations, such as Cape de Creus in the northwestern Mediterranean Sea, which harbors populations from 5 to 70 m depth [27,28]. In addition, Gori et al. [41] revealed that differences in reproductive cycle were observed between populations living at –20 m and–60 m depth, and this could suggest local adaptation depending on the bathymetry, which can therefore drive to a significant genetic differentiation. In tropical symbiotic corals, although all zooxanthellae belong to Symbiodinium sp., differences in clade composition are believed to confer a certain sensitivity or resistance to heat [20]. However, previous studies showed that nearly all Mediterranean symbiotic cnidarians harbour the A’ temperate Symbiodinium clade [21,22]. Nevertheless, we were very interested to test whether the intraclade diversity detected with symbiont microsatellite analyses could be linked to resistance phenotype. To date, our results did not indicate such a correlation (Table 2, upper high). On the other hand, we also inspected whether host genetic differentiation could explain the contrasting thermal sensitivity. In particular, the bimodal distribution of the Medes colonies could be due to the sampling of two different subpopulations. However, the Jost D value analyzes refuted this hypothesis. Several genotypic differences were correlated with individuals presenting differences in heat resistance, suggesting a potential link of the gorgon genotype to their phenotypic resistance. It must be remembered that genotypic analysis revealed a differentiation according to the geographic origin, which could explain, in part, these correlations between Riou and Medes Islands. However, among Medes Islands’ samples, the significant difference between the intermediate and resistant classes suggests that such a link exists. Since the microsatellite loci used in this study were neutral (and non-coding) and also considering our limited sample size, we could not formally identify the genetic markers correlated to the thermal sensitivity. It is clear that more samples and replications in each population are required to correlate a specific genotype with a thermo-sensitive or heat-resistant phenotype. In future studies, it would require the combination of a better genome coverage (as offered by next-generation sequencing technologies) and cell biomarker approaches, to answer this fundamental question: what genetic characteristics explain the phenotypic differences observed in the present study?