Geographical Variation in Phenotypic Plasticity of Intertidal Sister Limpet’s Species Under Ocean Acidification Scenarios

Introduction

Ocean acidification (OA) is a process where increasing amounts of anthropogenic CO₂ from the atmosphere are being diluted into the ocean, affecting its geochemical balance and decreasing its pH and carbonate saturation state (Ω) (Gattuso et al., 2015; IPCC, 2019). These changes affect the shell and skeleton formation of calcifying marine organisms (Riebesell et al., 2000; Orr et al., 2005; Ramajo et al., 2016). Moreover, as calcification has a core role controlling other processes such as growth, metabolism, and internal pH regulation, OA could have ubiquitous effects in these organisms (Pörtner and Farrell, 2008; Somero et al., 2016). In this context, to have a better understanding of how calcifying species are responding at the physiological and evolutionary level, it is crucial to model the impacts of this aspect of global change in marine ecosystems (Reusch, 2014; Sunday et al., 2014; Fox et al., 2019). It is increasingly documented that the marine animals’ responses to stressors such as ocean acidification can vary widely among and within species (Duarte et al., 2015, Shaw et al., 2016). Differential tolerance to environmental stressors is partly due to the conditions to which the individuals have been naturally exposed (Vargas et al., 2017), so it is important to consider those conditions before the experiments to understand the consequences of OA. Unfortunately, just a few investigations have considered the influence of geographical (e.g., Lardies et al., 2014; Gaitán-Espitia et al., 2017b; Broitman et al., 2018) and temporal environmental variation (Frieder et al., 2014; Eriander et al., 2015; Jarrold et al., 2017; Cornwall et al., 2018; Johnson et al., 2019) as drivers of differences in phenotypic plasticity or physiological responses among populations of marine organisms under OA.

The upwelling system along the southeast Pacific coast brings cold deep seawaters, rich in nutrients, and presents abnormalities in the surface temperature that produces oversaturation of CO₂ (Torres et al., 2011). The low pH of these waters shows a low degree of calcium carbonate saturation (Ω_{aragonite/calcite} < 1) (Broitman et al., 2018), which makes this region particularly prone to the effects of ocean acidification (OA) (Gruber et al., 2012; Vargas et al., 2017; Broitman et al., 2018). Waters in this naturally acidified coastal system can be corrosive and impact the calcification and physiology of benthic organisms, such as mollusks (Fabry et al., 2008; Ramajo et al., 2019). Thus, local populations of benthic organisms inhabiting this region can experience differences in upwelling regimens’ frequency and intensity (Ramajo et al., 2020). Specifically, Punta Lengua de Vaca (Talcaruca, Coquimbo region, 30°S), which is the most active wind-driven upwelling center along the Chilean coast (Rutllant and Montecino, 2002; Rahn and Garreaud, 2014), and also presents abnormalities in the surface temperature that produces oversaturation of CO₂ (Torres and Ampuero, 2009; Aravena et al., 2014). This natural scenario gives a unique opportunity to evaluate how the natural variability in the carbonate system affects the phenotypical and physiological characteristics of marine organisms (e.g., Lardies et al., 2017; Ramajo et al., 2019).

Phenotypic plasticity is defined as the capacity of a single genotype to exhibit different phenotypes depending on the environment (Schlichting and Pigliucci, 1998; Gianoli and Valladares, 2012), and in a variable environment, the level of plasticity can define the survival of a species or population (Schlichting and Pigliucci, 1998; Sultan, 2004; Lardies et al., 2014; Osores et al., 2017). Thus, in this study, it is hypothesized that a population thriving in a more variable environment will present higher phenotypic plasticity in physiological and morphological traits in response to different pCO₂, compared to a population of the same species from a more stable environment. To achieve this, two intertidal mollusk sister species were used as study models, the limpets Scurria zebrina and Scurria viridula. These species are common calcifying marine organisms and two of the most abundant herbivores in the mid rocky intertidal level, along the transitional biogeographic break situated at the southeast Pacific coast (30–32°S) (Espoz et al., 2004; Aguilera et al., 2013, 2020). They are sister species (Espoz et al., 2004) that share morphological traits (e.g., maximum size) and habitat but have contrasting behavioral responses to heat stress and predators (Espoz and Castilla, 2000, Broitman et al., 2018; Aguilera et al., 2020). Furthermore, it has been recently reported that these species show metabolic and growth plasticity in response to environmental changes (Aguilera et al., 2018; Broitman et al., 2018). We compared the plastic phenotypic responses between populations of these two species, inhabiting the aforementioned zone of Punta Lengua de Vaca, Talcaruca (where the rear and leading edges of the two sister species overlap; Aguilera et al., 2020), and an adjacent area, with more stable environmental conditions (Los Molles, northern Chile). Firstly, the relationship between shell length and calcium carbonate content was calculated in field specimens. Then, individuals from both populations were exposed to mesocosms representing the current pCO₂ in the Chilean coast (500 μatm) and the levels predicted for the year 2100 in upwelling zones (1500 (μatm) (Vargas et al., 2017). In these experimental conditions, the carbonate content, metabolism, and physiological traits such as hemocyanin content and growth rate were recorded in both populations.

Materials and Methods

Study Sites

Individuals of Scurria zebrina and Scurria viridula were collected from different populations in the upwelling center (i.e., biogeographic break) and in a near zone without the influence of semi-permanent upwelling in the northern-central Chilean coast (Figure 1): Talcaruca (30° 29′ S, 71°41′ W) and Los Molles (32° 24′ S, 71° 50′W), respectively. These sites are characterized by different variability in upwelling dynamics and the physical-chemical characteristics of coastal waters (Torres et al., 2011; Aravena et al., 2014; Vargas et al., 2017; Broitman et al., 2018). The environmental variables (i.e., temperature and salinity) and seawater carbonate system parameters were measured in the intertidal zone (see below and Table 1). For the temperature measurements, submersible temperature data loggers (HOBO^®, Onset Computer Corp., MA, United States) were installed, housed inside PVC pipes embedded in concrete blocks, and deployed at the mid-intertidal level. The loggers recorded the ocean temperature data every 30 min, which was downloaded on a monthly or seasonal schedule; see details in Gaitán-Espitia et al. (2017a). Using the same methodology and time-frequency described above, we anchored conductivity dataloggers (HOBO U24-001 Conductivity Data Logger, ONSET^®) to determine salinity in the intertidal zone (∼1 m below the seawater surface at the mid-intertidal). The pH in the field was measured biweekly for 2 years using a Metrohm 713 Meter (Metrohm^®) connected to a combined electrode (Metrohm^® model 6.0219.100) previously calibrated with the Metrohm^® pH 4 (6.2307.200), pH 7 (6.2307.210), and pH 9 (6.2307.220) standard buffers. The pH values are reported on the NBS scale. Furthermore, seawater samples were taken with a biweekly frequency (for 2 years) to measure Total Alkalinity (TA) as described later. Given this information, we proceed to estimate carbonate system parameters (see below for details of estimation) in both localities.

Figure 1. Location of the study sites localized near the upwelling center Punta Lengua de Vaca around 30°S on the Chilean coast. Map of the study sites (Talcaruca and Los Molles) in north-central Chile and their relation to the South America map. The distribution of both species is that reported in Aguilera et al. (2020), and it is indicated here with a red (Scurria viridula) or blue (S. zebrina) arrow. The distribution overlap zone for both species is marked with dashed lines. A representative photograph of both species also shows S. viridula in a red square and S. zebrina in a blue one.

Table 1. Summary of salinity and seawater carbonate chemistry variables measured in two study locations.

Animal Collection and Maintenance

Firstly, 24–30 individuals were collected during low tide from each species and population. These individuals were transported in coolers to the laboratory and then measured to determine the relationship between shell length and calcium carbonate content.

Subsequently, 40 individuals from each species were collected during low tide from each locality in the spring-summer season and were subjected to pCO₂ manipulation in experimental mesocosms as described below. The experimental design is detailed in Supplementary Figure 1. In the laboratory, the limpets were maintained in aquariums with seawater and each individual was marked and identified with bee tags (Beeworks^©) glued to the rear zone of the shell. The seawater for the aquariums was obtained from Quintay (33°11′S; 71°1′W, Valparaíso Region), and the animals were maintained with a salinity of 33 ppm and 14 ± 1°C (using a water bath SunSun) for 14 days of acclimation. Limpets were assigned to two different pH treatments: pCO₂; 500 μatm (7.9 pH), and 1500 μatm (7.5 pH), having four aquariums per treatment for each population in both species (16 aquariums in each limpet species). Each aquarium had four to five individuals from each species. The phenotypic traits were measured after 15 days in each treatment, and the limpets were moved to the other treatment for the next 15 days. Limpets were assigned randomly to different treatments and fed with Ulva spp. three times a week. The frame time of 15 days was chosen because the ecosystem influenced by upwelling systems, specifically in Punta Lengua de Vaca, shows that natural scales of variation in environmental variables (i.e., pH, pCO₂, and temperature) have an average duration of 13 days during upwelling activation (see Ramajo et al., 2019).

Experimental Mesocosms and Carbonate System Parameters

The pCO₂ was achieved in experimental mesocosms as described in Torres et al. (2013) and Benítez et al. (2018). In brief, 500 μatm and 1500 μatm were obtained by pumping dry air and pure CO₂ until the target was achieved using a mass flow controller (MFC, Aalborg, model GFC) as described in Benítez et al.(2017; 2018). The air injected into the air MFC controller was previously dried using a drierite desiccant column (W. H. Hammond Drierite Co), filtered, and the CO₂ was extracted using a soda-lime column. The pH, salinity, and temperature were controlled daily. Treatment conditions in the aquariums were maintained bubbling with the corresponding CO₂ concentration. Total or partial water changes were performed according to the registered pH. The pH samples were collected in 50 mL syringes and immediately transferred to a 25 mL temperature-controlled cell maintained at 25.0 ± 0.1°C for standardization. The pH was measured each day using a Metrohm 713 Meter (Metrohm)^® connected to a combined electrode (Metrohm^® model 6.0219.100) previously calibrated with the Metrohm^® pH 4 (6.2307.200), pH 7 (6.2307.210), and pH 9 (6.2307.220) standard buffers. The pH values are reported on the NBS scale.

Samples for Total Alkalinity (TA) were poisoned using HgCl₂50 μL, and the bottles (Pyrex^®, Corning) were sealed using parafilm. Then, TA was determined using the open-cell titration method (Dickson et al., 2007) using an automatic alkalinity titrator (Model AS-ALK2, Apollo SciTech) equipped with a combination pH electrode (8102BNUWP, Thermo Fisher Scientific) and temperature probe (Star ATC, Thermo Fisher Scientific) connected to a pH meter (Orion Star A211, Thermo Fisher Scientific). All samples were analyzed at 25°C (± 0.1°C), and the temperature was regulated using a water bath (Lab Companion CW-05G). Accuracy was controlled using certified reference material (CRM, supplied by A. Dickson, University California San Diego), and the TA repeatability was 2 to 3 μmol kg^–1 on average. Temperature, salinity, pH, and TA data were used to calculate the carbonate system parameters (pCO₂, CO₃^2–). To achieve this, analyses were performed using CO2SYS software in MS Excel (Pierrot et al., 2006) set with Mehrbach solubility constants (Mehrbach et al., 1973) refitted by Dickson and Millero (1987). The carbonate systems in the mesocosms for both localities are listed in Supplementary Tables 1, 2.

Physiological and Metabolic Measurements

Calcium carbonate content (net calcification) of the limpets was estimated from changes in their buoyant weight (i.e., underwater weight) and verified with the dry weight measurements of the shells (Palmer, 1982). All the weight measurements were conducted using an analytical balance (Shimadzu AUX220). For the field samples, the calcium carbonate content and shell length were analyzed immediately once they arrived in the laboratory, under controlled conditions in seawater parameters (see above). For the individuals in the mesocosm, estimations were made on day 0 (at the beginning of the experiment) and at the end of a 15-days period exposition to each treatment. Shell dissolution rates, in turn, were estimated by weighing 2–3 empty shells of limpets (per treatment) at the beginning and the end of the experiment (see Duarte et al., 2015). The inorganic material content was calculated by putting the shells in a furnace (Vulcan model A-550) for 4 h at 450–550°C and weighing the remainder with an analytical balance.

Growth rates were estimated from changes in the total body weight (to the nearest 0.0001 g) of the limpets and expressed as a daily rate of increase. Before oxygen uptake measurements, all individuals were kept under starvation for 24 h in their specific experimental conditions (i.e., pCO₂ concentration, temperature, salinity, and photoperiod). Then, limpets were placed individually into a closed respirometric chamber (113 mL), filled with filtered seawater from the corresponding pCO₂ treatment. Once sealed, the chambers were placed into a tank with seawater at 14°C controlled by a chiller. Oxygen consumption rates were measured using a fiber-optic oxygen optode connected to a PreSens Microx TX3 temperature compensated oxygen meter (Precision Sensing, GmbH, Regensburg, Germany) with a tip diameter of 140 mm. Oxygen partial pressure measurements were run for at least 60 min and were never allowed to decrease below 85% O₂ saturation to avoid animals experiencing hypoxia. The first 10 min and the last 5 min of determinations were eliminated to avoid possible disturbances caused by the stress of animal manipulation. Oxygen consumption values were expressed as O₂ mg h¹g^–1. Metabolic and growth rates were calculated for at least 8 individuals from each population and species in each treatment.

Hemocyanin content was analyzed as described in Pascual et al. (2003), extracting 10 μL of hemolymph from the foot using a hypodermic syringe with EDTA anticoagulant solution and quantifying in a spectrophotometer (GeneQuant 1300 Healthcare BioSciences) at 325 nm and using the equation:

Hemocyanin concentration (Hc) = (Absorbance/extinction coefficient) x dilution factor.

The phenotypic plasticity index (PI_md) was calculated as indicated in Valladares et al. (2006) for the metabolic and growth rate traits. Briefly, it was estimated comparing each species and population separately, considering the maximum medians (absolute values) as follows:

PI_md = [(| Median_max|−| Median_min|)/| Median_max|] x 100.

Statistical Analyses

All data analyses were performed using GraphPad Prism version 8.4.3, GraphPad Software, San Diego, CA, United States. A simple linear regression was used to describe the relationship between the length (L) of the individuals and the carbonate content (CC) in each species from each population using data from at least 24 individuals per treatment. The slopes were compared between each population. Metabolism, growth rate, and hemocyanin content were compared using a two-way ANOVA with pCO₂ and population as fixed factors, and a posteriori Sidak’s multiple comparisons test. Data were previously evaluated for normality using the Kolmogorov-Smirnov test and heteroscedasticity with the Spearman’s test.

Results

Environmental Conditions in Localities and Their Relationship With the Carbonate Content of Both Scurria Species

Sea Surface Temperature (SST) in Los Molles presents an annual media of 14.24°C, with an annual maximum of 16.40°C and a minimum of 10.80°C (see Table 1). On the other hand, Talcaruca has an annual media of 13.20°C, with an annual maximum of 16.8 °C and a minimum of 8.40°C; the coefficient of variation was 38% higher in Talcaruca than in Los Molles (Table 1). Seawater carbonate chemistry parameters in samples collected from the field varied highly among study sites (Table 1). The highest and more variable levels of pCO₂ were observed in Talcaruca (870 ± 312.9 μmol/Kg), while a pCO2 of 458 ± 141.4 (μmol/Kg) was observed in Los Molles. Talcaruca also exhibited lower pH values than Los Molles (7.93± 0.43 vs. 8.08 ± 0.11, respectively). The saturation states for calcite (_Calcite) and aragonite (Ω_Aragonite) were 27.14 and 22.5% lower in Talcaruca than in Los Molles, respectively. Salinity and TA (Total Alkalinity) presented similar values in both sites (see Table 1).

In order to analyze how the locality affects the ratio carbonate content/shell length, a simple linear regression was performed comparing these parameters in each species (Figure 2). In the case of S. zebrina, slopes were significantly different from zero (Talcaruca: CC = 0.1006^∗L - 1.330; F₍₁,₂₂₎ = 6.69, p = 0,0168 and Los Molles: CC = 0,3819^∗L - 8,354; F₁,₂₈ = 147.8, p < 0.0001). Furthermore, slopes were significantly different between both populations (F₁,₅₀ = 29.39, p < 0.00001) (Figure 2A). In Talcaruca, individuals presented more variability, and the maximum values of calcium carbonate were lower than in Los Molles (Figure 2A). For S. viridula, slopes were also significantly different from zero (Talcaruca: CC = 0.2418^∗L -- 4.990; F₁,₂₅ = 160.3, p < 0.0001 and Los Molles: CC = 0,2811^∗L - 5,566; F₁,₂₄ = 90.8, p < 0.0001). In this case, slopes from both populations were not significantly different (F₁,₅₀ = 1.26, p = 0.2671), and a similar variability was observed in both populations (Figure 2B).

Figure 2. Carbonate content in populations of two congeneric Scurria species. Simple linear regressions (lines) of length and calcium carbonate (Log₁₀ transformed) in S. zebrina (A) and S. viridula (B) from two different localities, Talcaruca (empty circle), and Los Molles (black square). Each dot represents one organism.

pCO₂ Effects in Growth and Physiological Traits of Two Different Scurria Species

To understand how contrasting pCO₂ levels could induce different effects in the metabolism of the populations of Talcaruca and Los Molles, two mesocosms were installed to reach 500 and 1500 μatm of CO₂. The carbonate systems in the experimental mesocosm for Talcaruca and Los Molles are detailed in Supplementary Tables 1, 2, respectively. Interestingly, in both species, individuals of Talcaruca showed more variability and a higher mean in metabolic rate in both conditions than those from Los Molles (Figure 3). Two-way ANOVA analyses were performed considering the pCO₂ and population as factors. For S. zebrina, metabolic rates were not significantly different between pCO₂ treatments (Figure 3A; F₁,₃₆ = 0.34, p = 0.56), however, the mean metabolic rates and variation of the populations from Talcaruca were higher than those from Los Molles, regardless of pCO₂ treatment (Figure 3A; population factor F₁,₃₆ = 11.97, p = 0.0014). Additionally, a significant effect was observed at 1500 μatm of CO₂ between the two populations (p = 0.023). A similar pattern was observed in S. viridula, where both the pCO₂ and population factors had significant effects on the metabolic rate (F₁,₃₀ = 4.7, p = 0.037; and F₁,₃₀ = 11.36, p = 0.0021, respectively). Again, a significant effect was observed at 1500 μatm of CO₂ between the two populations (p = 0.0229) (Figure 3B). The phenotypic plasticity index (PI_md) was compared between the pCO₂ treatments in each species and population. For S. zebrina, this index for metabolism was higher (12.14%) in Talcaruca than in Los Molles (5.9%). The opposite pattern was observed in S. viridula, where the PI_md was higher in Los Molles (89%) than in Talcaruca (43%).

Figure 3. pCO₂ effects on the metabolic rate of two different Scurria species. Box and whiskers plot representing the metabolic rate of S. zebrina (A) and S. viridula (B) under 500 or 1500 μatm of CO₂. Data from Talcaruca is presented in white and from Los Molles in black. The boxes represent the interquartile range, while the inner line indicates the median value, and the vertical lines represent the maximum and minimum values. Asterisks represent statistical differences (Two-way ANOVA and Sidak multiple comparison test, p < 0.01).

Growth rates were also analyzed with two-way ANOVA (pCO₂ and population as factors), finding no significant effects of the interaction, both in S. zebrina (Figure 4A) and S. viridula (Figure 4B) (F₁,₃₃ = 1.28, p = 0.26, and F₁,₂₈ = 0.35, p = 0.55, respectively). Alike, no significant effects were detected for the pCO₂ or population factors in both species (Figures 4A,B). PI_md for growth rates was higher for both species in Talcaruca than in Los Molles, 39.9% versus 12.9% in S. Zebrina, respectively; and 91% versus 75.6% for S. viridula, respectively.

Figure 4. pCO₂ effects on the growth rate of two different Scurria species. Box and whiskers plot representing the growth rate of S. zebrina (A) and S. viridula (B) under 500 or 1500 μatm of CO₂. Data from Talcaruca is presented in white and from Los Molles in black. The boxes represent the interquartile range, while the inner line indicates the median value, and the vertical lines represent the maximum and minimum values.

The differences in the dissolution rate of shells belonging to the two different populations were measured, finding that the highest dissolution rate was observed at 1500 μatm for the shells from Los Molles in both species, S. zebrina (Figure 5A) and S. viridula (Figure 5B).

Figure 5. pCO₂ effects on the dissolution rates of two different Scurria species. Bars represent the mean and the standard error mean in the dissolution rates of S. zebrina (A) and S. viridula (B) under 500 or 1500 μatm of CO₂. Data from Talcaruca is presented in white and from Los Molles in black.

Finally, the hemocyanin content did not show significant differences between populations and treatments, but a higher concentration of hemocyanin in S. viridula (Figure 6B) than in S. zebrina was observed (Figure 6). Additionally, both species presented a different pattern regarding the population, where S. zebrina showed a higher value when coming from Los Molles (Figure 6A), and S. viridula presented higher values when coming from Talcaruca (Figure 6B).

Figure 6. pCO₂ effects on the hemocyanin content of two different Scurria species. Bars represent the mean and the standard error mean in the hemocyanin content of S. zebrina (A) and S. viridula (B) under 500 or 1500 μatm of CO₂. Data from Talcaruca is presented in white and from Los Molles in black.

Discussion

Phenotypic plasticity is a primary way by which organisms respond to environmental variability, including those imposed by global change (Bonamour et al., 2019), and phenotypic changes registered in natural populations have been correlated with climate change (Charmantier et al., 2008; Reusch, 2014). Still, the distinction between genetically or plasticity-based changes is not always provided (Merilä and Hendry, 2014). Here we combined field characterizations with experimental approaches to analyze the metabolic and physiological responses of two populations from contrasting environments, using two congeneric species of limpets.

Both study sites, despite being separated on a small spatial scale (i.e., approximately 200 km), differ abruptly in their environmental conditions. These differences were mainly due to the semi-permanent upwelling center at Punta Lengua de Vaca, which is localized on the coast of Talcaruca (Torres et al., 2011; Lardies et al., 2017; Broitman et al., 2018). Temperatures near the coast in Talcaruca are influenced by prolonged cold events associated with upwelling; these factors are important since such changes could affect the dynamics and structure of ecological communities (Aravena et al., 2014; Broitman et al., 2018). Our results indicate that the salinity levels in both locations were similar, whereas the temperature, pH, and carbonate system parameters (i.e., states of saturation of aragonite and calcite) showed abrupt spatial variation, characterizing this coastal zone as highly patchy in terms of temperature and carbonate system parameters.

Most experimental studies reveal that growth rates and calcification are negatively affected under OA conditions in a wide range of taxa (Hofmann et al., 2010; Kroeker et al., 2013; Duarte et al., 2015; Bednarsek et al., 2019). Although, a few studies show positive effects on growth (e.g., Wood et al., 2008; Findlay et al., 2009; Gutowska et al., 2010; Lardies et al., 2017) and calcification rates (Wood et al., 2008; Ries et al., 2009; Lagos et al., 2016). This suggests that the effects of OA on shell growth and carbonate precipitation are species-specific (Ries et al., 2009; Kroeker et al., 2013), and that carbonate layers can develop at low pH and under subsaturation conditions, namely below the threshold of Ω_aragonite = 1 proposed in other studies (Findlay et al., 2011). OA affects the pH in body fluids, such as the extrapallial fluid, affecting calcification and shell structures (see Cohen and Holcomb, 2009; Ries, 2011).

In general, we observed that the changes in phenotypic traits in scenarios of ocean acidification correlated with the environmental variability experienced in the native habitat of the populations. The shells’ carbonate content showed higher variability in S. zebrina individuals from Talcaruca than those from Los Molles, and some individuals from Los Molles reached larger body sizes. Interestingly, this correlates with the fact that S. zebrina is on the edge of its range in Talcaruca, where the species has a north endpoint of distribution and is not able to overstep this biogeographic break of the Chilean coast (Espoz et al., 2004; Rivadeneira and Fernaìndez, 2005; Aguilera et al., 2020). On the other hand, S. viridula showed a similar carbonate content/size relation in the two analyzed populations (see Figure 2B). In smaller individuals, higher quantities of carbonate content were observed in those individuals from Talcaruca, but the individuals of greater sizes showed similar carbon content ratios, and interestingly, similar maximum sizes were detected in both localities. This indicates that S. viridula can compensate the effects of this variable environment along its ontogeny, in terms of growth rate and calcification, and consequently surpass the biogeographic break of the Talcaruca site, extending its distributional range to the northern coast (see Aguilera et al., 2013). Growth rates estimated in experimental mesocosm did not significantly differ among populations. Still, data were less variable in S. viridula than in S. zebrina (see Figures 4A,B, respectively) and S. viridula seems to be less affected in the acidified conditions than the control, both in the population of Talcaruca and in Los Molles. Greater differences in phenotypic plasticity between populations have been implicated as key contributors to persistence under challenging or new environmental conditions, especially under ocean acidification (Murren et al., 2014). The degree of phenotypic plasticity can vary greatly between related species (Gibbin et al., 2017; Hattich et al., 2017), as was registered in S. viridula and S. zebrina, and can be influenced by the spatial differentiation of the niche or even by competition (Broitman et al., 2018; Aguilera et al., 2020). This is even more relevant for these species, which thrive in areas with semi-permanent upwelling, generating higher physiological demand, particularly in larval stages (Doney et al., 2012). On the other hand, the intertidal zone already exhibits considerable variability in pH (Hofmann et al., 2011; Jellison et al., 2016) and temperature (Aravena et al., 2014; Helmuth et al., 2016), but ocean acidification and ocean warming are likely to result in greater intertidal zone extremes in temperature and pH in the future (Solomon et al., 2007; Paganini et al., 2014), which could add significant physiological stress on intertidal marine fauna. In this sense, our experimental setup where limpets were submerged during the period of acclimation and experimental could be underestimating the effects of OA in physiological process compared with variability experienced in their natural habitat. In this sense, in our experimental setup, limpets were submerged and not exposed to intertidal cycles, which could cause an underestimation of the effects of OA in the physiological processes in their natural habitat. Despite this, as both populations were exposed to the same conditions, our experimental setting allowed us to detect differences in their plastic responses (see also Broitman et al., 2018).

The pattern of energy usage of an organism is reflected in energy expenditure measurements, the most common being the rate of metabolism. Furthermore, physiological variation within an individual’s life history can have profound implications for fitness (Lardies and Bozinovic, 2006; Koop et al., 2011; Ramajo et al., 2020). Increasing evidence suggests that organisms associated with environments that naturally present high pCO₂ levels may have physiological and metabolic adaptations and consequently be better acclimatized to ocean acidification (e.g., Duarte et al., 2015; Calosi et al., 2017; Vargas et al., 2017; Ramajo et al., 2019). Biomineralization in calcifying organisms is a highly regulated and potentially energy-demanding biological process (Sokolova et al., 2012); therefore, organisms exposed to hypercapnic conditions face higher energy demands to supply processes that are vital for the maintenance of cellular and extracellular homeostasis and/or maintenance of calcareous integrity (Stumpp et al., 2011). It has been shown that low pH and the associated reduction in carbonate ion availability can affect multiple physiological pathways in bivalves, beginning with calcification (Gobler et al., 2014). Concordantly, experiments carried out with Mytilus edulis exposed to acidic conditions responded with higher metabolic rates at the expense of a decrease in calcification and growth rates (Thomsen et al., 2010). Nevertheless, Milano et al. (2016) suggest that organisms react differently to acidified waters and that different responses may depend on the specific capacities of organisms to regulate their internal pH (Ries et al., 2009; Gutowska et al., 2010). Here, we found that the metabolic rate was highly dependent on the population habitat, but in general, was lower in Los Molles than in Talcaruca (Figure 3). When individuals were exposed to acidification, those of S. zebrina maintained or even increased their metabolic rate, whereas the individuals of S. viridula were more likely to decrease their metabolic rate at higher pCO₂ levels.

Besides, the shell dissolution rate was lower for both species in limpets from the Talcaruca population, and especially in S. viridula, the shells were more resistant to dissolution. Comparable results have been reported in other studies of mollusks from estuarine and upwelling zones (e.g. Duarte et al., 2015; Lagos et al., 2016; Lardies et al., 2017; Ramajo et al., 2020), and also in other mussel species that appear resilient to elevated pCO₂ (e.g., Thomsen and Melzner, 2010; Mackenzie et al., 2014; Clements et al., 2018). However, the changes in empty shell weight were used as a proxy of shell dissolution but must be interpreted with caution, given the role of organic layers such as the shell periostracum protecting live animals from dissolution (Tunnicliffe et al., 2009; Nienhuis et al., 2010).

Hemocyanins in the hemolymph of marine invertebrates are functional analogs of the hemoglobin present in vertebrates and exhibit cooperativity in their oxygen-binding properties and allosteric regulation to efficiently uptake and deliver oxygen under physiological conditions (Taylor et al., 1995; Birk et al., 2018). The results obtained showed no significant differences in the hemocyanin content among treatments. Despite this, it is important to note that S. viridula individuals from Talcaruca could produce a higher concentration of hemocyanin at higher levels of pCO₂, compared to those from Los Molles (see Figure 6B), which is consistent with the levels of oxygen consumption recorded during the experiment (Figure 3B).

We observed contrasting phenotypic responses between the populations of S. zebrina and S. viridula, across the biogeographic boundary where their geographic distributions overlap; therefore, local environmental conditions are key in the maintenance of coastal marine patterns of distributional ranges. Nevertheless, it is worth mentioning that there is an inherent limitation in scope between the time-scale of the incubation period considered in our study (i.e., 15 days) with the scales at which natural variation occurs, and the projected event of OA will take place in nature. Nevertheless, an ecosystem influenced by an upwelling system, specifically on Punta Lengua de Vaca, shows that natural scales of variation in environmental variables (i.e., pH, pCO₂, and temperature) have an average duration of 13.8 days during upwelling activation (see Ramajo et al., 2019). Despite these discrepancies, such experiments provide insights into stressor tolerance and adaptation levels, mainly because the rates of acclimation to hypercapnia have been reported to occur especially fast in some marine invertebrates (Widdicombe and Spicer, 2008; Calosi et al., 2013; Ramajo et al., 2019). Then, it is probable that the observed differences in studied physiological traits are the end-product of local adaptation.

Taken together, our results from the field and laboratory analyses show that responses to OA are species-specific, even in this related species. Interestingly, S. viridula showed better performance under OA than S. zebrina (i.e., similar sizes and carbonate content in individuals from both populations; lower effects of acidification on the growth rate combined with a reduction of metabolism at higher pCO₂). Interestingly, these characteristics could explain this species’ success in overstepping the biogeographical break of Talcaruca, which S. zebrina cannot achieve. Additionally, the results show that the habitat factor has a strong influence on some traits. For instance, individuals from Talcaruca presented a higher growth rate plasticity index than those from Los Molles. Although we cannot conclude (within the experimental frame time) that this higher plasticity is associated with a higher fitness or performance, the literature indicates that this could be related to the maintenance of the physiological performance in individuals exposed to acidic waters (Osores et al., 2017; Ramajo et al., 2020). According to the literature, differences in populations’ metabolic responses might be related to their aerobic capacity as an adaptation to different habitats (Clarke and Fraser, 2004; Watson et al., 2014). Therefore, compared populations that occupy more stable environments (Los Molles) with populations that inhabit more variable environments appear to be more tolerant to extreme acute acidification events because of their higher metabolic limits (Osores et al., 2017; Jahnsen-Guzmán et al., 2020; Kurihara et al., 2020). However, many studies indicate that these populations may be at higher risk because of the adverse effects of chronic exposure to acidification, where higher metabolic costs may imply trade-offs (i.e., a decrease in calcification) (Lannig et al., 2010; Lagos et al., 2016; Osores et al., 2017). Thereby, this work highlights the relevance of incorporating the biogeographical dimension to have a better understanding of the environmental challenges that the marine organisms are already facing and brings insights to model the future effects of OA in different populations and species.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

ML acquired funding for the study. ML and CD conceived the idea and designed the study. MP, ML, and PC carried out experiments and analyzed the data. All authors drafted and commented on the manuscript.

Funding

This work was funded by ANID FONDECYT No. 1190444 and FONDECYT No. 1190710. ML and CD appreciate the support given by the Millennium Nucleus Center for the Study of Multiple drivers on Marine Socio-Ecological Systems (MUSELS, MINECON NC120086) and Programa Investigación Asociativa (ANID ACT 172037). Further, support by the Millenium Institute SECOS to ML. is greatly appreciated. MP also thanks ANID PIA/BASAL FB0002 and ANID-FONDECYT No. 1190634.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank Carolina Fernández for field and laboratory assistance and Nelson Lagos for providing equipment for alkalinity measurements. We would also like to thank Laura Ramajo and Roberto García-Huidobro for their help with carbonate system parameter estimations and Bernardo Broitman for providing the maps in Figure 1.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.647087/full#supplementary-material

References