Local scale pH conditions, genetic diversity and copepod fecundity-related traits were characterized in a temperate-estuarine system in Southern Chile (Valdivia River Estuary, 39.8°S 73.2°W) and a coastal upwelling system off northern Chile (Antofagasta, -23.4°S -70.6°W). Strong riverine discharges characterize the estuarine system, with maximum runoff during the austral autumn-winter period (130 ± 93 mm y^-1 of precipitation), and occasional rainy events occurring throughout the year (Pérez et al., 2016). The influence of freshwater runoff in this estuarine system is widely affected by tidal cycles that determine regular events of low temperature (< 14°C), low pH (pH_T < 7.7) and low salinity conditions (< 33 psu) in the adjacent coastal zone (Aguilera et al., 2013; Garcés-Vargas et al., 2020; Osma et al., 2020). In contrast, the Antofagasta upwelling system is located off the most arid global desert, the Atacama Desert, which lacks seasonality and freshwater discharges (Hartley et al., 2005). Permanent equatorward winds promote year-round uplifts of subsurface cold and low pH (pH_T < 7.7) water (Torres et al., 2002; Aguilera et al., 2020a) over an extremely narrow (< 5 km) continental shelf. The estuarine and upwelling habitats are embedded in the Humboldt Current System and located in distant coastal provinces (> 1700 km, > 15° latitude distant), separated by oceanographic clines efficiently structuring the latitudinal distribution and genetic diversity of benthic invertebrates (Haye et al., 2014).

Environmental characterization in the estuarine system was assessed during the seasonal (austral spring-summer) occurrence and reproduction of the local A. tonsa population, between the years 2010 and 2012, under the influence of flood tides (Aguilera et al., 2013). The environmental variability was characterized at 7 m, likely reflecting the most prevalent environmental niche of adult females in which fecundity-related traits were assessed. Indeed, the spatial niche of A. tonsa might be smaller in the estuarine system due to hydrographic conditions potentially reducing the habitat suitability on the surface (less dense and low pH freshwater lens) and near the bottom (i.e., turbulent and sediment enriched tidal currents) (Pino et al., 1994; Vargas et al., 2003; Osma et al., 2020). Environmental drivers such as temperature, salinity, pH_T, and food concentration were characterized on an inter-daily scale (i.e., each 4 – 6 days) by means of 18 oceanographic surveys ( Table 1 ). A small CTD (Ocean Seven 305 Plus, www.idronaut.it) deployed from above the bottom to the surface provided continuous temperature (°C) and salinity (psu) measurements. In addition, water samples were collected at 7 m depth with a 10-L Niskin bottle for determinations of pH_T, total alkalinity (A_T), and copepod food. Within 2 h of collection, pH_T was measured in a closed 25-mL cell thermostated at 25.0°C using a Metrohm 713 pH meter and a glass combined double junction Ag/AgCl electrode (Metrohm model 6.0219.100) calibrated with Tris buffer at 25°C. Samples for A_T were collected in 250 mL high-density polypropylene bottles, poisoned with 50 µL of a saturated HgCl₂ solution and then stored in darkness at room temperature until analysis no longer than 1 month. A_T was analyzed using the automated potentiometric titration method (Haraldsson et al., 1997) controlling the accuracy against certified reference material (CRM, Batch #101, Scripps Institution of Oceanography, San Diego, USA). With pH, A_T, temperature, and salinity data, pH_T (total scale) and other carbonate system parameters were estimated using the CO2sys_v3.0 software (Pierrot et al., 2021), set with Mehrbach solubility constants (Mehrbach et al., 1973) and refitted after Dickson and Millero (1987). Food concentration was determined as the biomass of nanoflagellates and phytoplankton cells in carbon units (μg C L^-1) (e.g., Vargas and González, 2004).

Environmental characterization in the upwelling system was carried out at variable intervals (4 ± 3 d) during 2015 ( Table 1 ). In this system, A. tonsa reproduces throughout the year (Ruz et al., 2015) and its vertical niche seems to be constrained to the upper 20 m of the water column due to a shallow oxygen minimum zone (Escribano et al., 2009). Thus, environmental characterization in the upwelling system was accomplished at 10 m depth by means of 28 oceanographic surveys, all during 2015 ( Table 1 ). Temperature and salinity vertical profiles were recorded from above the bottom (∼ 30 m) to the surface, through vertical deployments of a calibrated SeaBird SBE19 Plus CTD. Seawater was collected through oceanographic sample collections at 10 m depth for determinations of pH_T, A_T and copepod food. Seawater pH_T was measured in triplicate as described above, within 1 h from the time of collection in a closed 25 mL cell thermostated (25°C). A_T samples were handled and determined as described above by means of an automated Alkalinity Titrator AS-ALK2 Apollo SciTech. The accuracy for A_T determinations was controlled in combination with certified reference material (A. Dickson, Batch #140). The pH_T was calculated as described above. Food concentration was determined as chlorophyll-a (Chl, µg Chl L^-1) concentration, and then Chl concentration was converted to carbon units by using a Chl:C ratio=120, according to Vargas and González (2004), for the same study area. Triplicate samples (200 mL) of sieved (200 μm mesh) seawater were filtered onto a GF/F filter (nominal pore size=0.7 µM). Chl was extracted for 24 h in 90% acetone v/v and measured in a TD Turner fluorometer (Strickland and Parsons, 1972).

Three different cohorts of both the estuarine and upwelling A. tonsa populations were assessed. Three different cohorts were assessed in the estuarine habitat over the period of 2010–2012, while in the northern-upwelling habitat, copepods were sampled throughout 2015. Cohorts were identified because they were either sampled in subsequent seasons (estuary) or by accumulated changes in the body size of adult females across an annual period (upwelling) (Aguilera and Bednaršek, 2022). Copepod samples were collected by gentle oblique hauls with a 200-µm mesh size WP2 net equipped with a non-filtering 1 L cod-end, from 12 to 7 m depth (estuary) and 15 to 10 m depth (upwelling) strata. Samples were gently transferred to a temperature-controlled and well oxygenated container over the duration of the transport to the laboratory.

Genetic diversity was assessed using females randomly sorted from samples collected during surveys encompassing months (upwelling, n=29) to years (estuary, n=32) and analyzed individually. For molecular analyses, partial sequences of the mitochondrial gene cytochrome c oxidase subunit I (mtCOI, c. 700bp; Forward LCO1490: 5’-GGTCAACAAATCATAAAGATATTGG-3’ and Reverse HCO2198: 5’-TAAACTTCAGGGTGACCAAAAAATCA-3’; Folmer et al., 1994) and nuclear gene for 18S ribosomal RNA (c.1600 bp; Forward 18A1mod: 5’-CTGGTTGATCCTGCCAGTCATATGC-3’) and Reverse 1800mod: 5’-GATCCTTCCGCAGGTTCACCTACG-3’; Raupach et al., 2009) were amplified through the polymerase chain reaction (PCR). The mtCOI gene is related to the maternal mode of inheritance with sharp intraspecific (up to 4%) and interspecific (∼9 to >25%) genetic divergence (Bucklin et al., 1999; Bucklin et al., 2003; Hebert et al., 2003b; Bucklin et al., 2010). The nuclear 18S gene evolves at a lower rate than mitochondrial genes (Moriyama and Powell, 1997), which has been a valuable complementary tool to disentangle complex phylogenetic relationships and reconstruction in copepod species (Blanco-Bercial et al., 2011; Cornils and Blanco-Bercial, 2013). All PCR amplifications were performed in 30 µl volume reaction with 3µl PCR buffer (Promega), 3µl MgCl₂ (Promega), 0.5 µl dNTPs (10 mmol µl^-1), 0.5 µl of each primer (100 pmol µl^-1), 0.2 µl Taq polymerase (Promega), 17.3 µl of sterile water and 5 µl of DNA template. For mtCOI amplification, a temperature profile with an initial denaturation at 94°C for 1 min was applied, followed by 32 cycles at 94°C for 40 s, annealing at 48°C for 40 s, extension at 72°C for 1 min, and finishing with another extension at 72°C for 1 min. Amplification of the 18S fragment consisted in an initial denaturation at 94°C for 1 min, then 37 cycles at 94°C for 40 s, annealing at 54°C for 1 min, extension at 72°C for 1 min and a final extension at 72°C for 90 s. Successful amplification of fragments was confirmed by 1.2% agarose gel electrophoresis and staining with GelRedTM fluorescent DNA dye. Both DNA strands were directly sequenced (Macrogen, Seoul, Korea; http://www.macrogen.com). Mitochondrial and nuclear sequences were edited and aligned in Geneious R11.1.5 (2018) using a MUSCLE alignment algorithm (Edgar, 2004) with default settings. Population diversity indexes, as the number of polymorphic sites (p), number of mutations (nm), nucleotide diversity (π), mean number of nucleotide differences (k), haplotype diversity(h), number of haplotypes (nh) and number of private haplotypes (np), were estimated in DnaSp software V.6. (Rozas et al., 2017). Genetic divergences between populations (Φst) were estimated in the software Arlequin V.3.5 with 10000 permutations (Excoffier and Lischer, 2010). Finally, genealogical relations between haplotypes were represented with a minimum spanning haplotype network (Bandelt et al., 1999) in PopART software (Leigh and Bryant, 2015).

Assessed fecundity-related traits were female body size, egg production rate and egg size. Within 2 h of collection, only mature and visibly healthy A. tonsa females were sorted out under a stereomicroscope ( Supplementary Figure 1 ). Up to 40 copepod females were immediately preserved in 90% ethanol for further measurements of copepod body length (mm) within two weeks of preservation to diminish ethanol effects on body size (Moksness and Fossum, 1992). Additional copepods were incubated in seawater collected during sampling with the aim of estimating phenotypic plasticity in fecundity-related traits as related to the habitat-specific variations in environmental conditions. Mean incubation temperature was 14 ± 1°C, which corresponded to the mean temperature conditions experienced by the various cohorts in both habitats. Mean egg production rates (EPR) were estimated based on batches of 30–40 copepod females, individually incubated for 17–20 h in 200 mL clean crystallizing dishes filled with sieved (<200 μm) natural seawater. The concave walls of the dishes that converge towards the floor, allow the settlement and grouping of the relatively dense copepod eggs (Tang et al., 1998) onto the floors of the dishes (although the walls of the dishes were surveyed as well). Eggs produced over this period by each group were counted under a stereomicroscope, standardized to daily duration (24 h), with the mean EPR expressed as the egg fem^-1 d^-1 (± SD) (Aguilera et al., 2011). After being counted, some of the produced eggs (20–30) were preserved (90% ethanol) and their size (i.e., diameter in µm) measured under an inverted microscope within two weeks of preservation to diminish ethanol effects on egg size (Moksness and Fossum, 1992).

Environmental conditions in acclimation experiments should be consistent with the habitat conditions of the respective populations (Stillman, 2003), providing the opportunity of examining phenotypic outcomes in response to specific attributes of environmental variability, such as extreme conditions. This contrasts with common garden experiments in which the mean habitat conditions of one (or more) population is displaced to a common nominal level to allow comparison of phenotypic responses among different populations (Thorpe et al., 2005). Three acclimation experiments with three different cohorts were carried out within each population immediately after estimating phenotypic responses to field conditions. Females were acclimated during 96-h under mean field temperature and salinity, controlled food supply, and high (i.e., control) and extreme low pH levels ( Supplementary Figure 3 ). Events with pH < 7.7 were considered extreme low pH events (Bednaršek et al., 2021). This threshold of pH 7.7 corresponds to the 5th percentile of measurements at the estuarine and upwelling habitats. These conditions are associated either to the predominance of freshwater in the estuary (Aguilera et al., 2013; Osma et al., 2020) or very recently upwelled waters (Torres et al., 2002). In the estuarine system, seawater was collected 8 km northward from the estuarine system and transferred into a micro-mesocosm laboratory at the Calfuco Marine Station, Austral University of Chile. In the upwelling system, seawater was collected during three different opportunities from the Antofagasta Bay and transferred to laboratory facilities at the Marine Sciences Faculty of the Antofagasta University. Mean seawater temperature was maintained relatively constant (± 0.2°C) in a free circulating water system (estuary) or cold room (upwelling) during the incubations. In both systems, the target pH_T levels were achieved by mixing seawater with air containing different pCO₂ concentrations (Vargas et al., 2017). Mean control and extreme pH values reproduced during acclimation experiments were 7.73 ± 0.05 and 8.0 ± 0.06 in both the estuarine and upwelling system. Copepod females (4 – 6) were pipetted into three 660-mL borosilicate acid-washed bottles filled with pH equilibrated water, with daily exchange of incubation water and food, the latter composed of Isochrysis sp. (estuary) and Isochrysis sp. + Tetraselmis suecica (upwelling) in concentrations above the saturation level for this species (>350 µg C L^-1 and > 4 µg Chl L^-1; Thompson et al., 1994). After the last 24 h incubation within the 96 h of acclimation, fecundity-related traits (i.e., EPR) were estimated as described above.

To test differences in environmental conditions at the inter-population level (i.e., between A. tonsa populations) we used a one-way ANOVA, while body-size dependent fecundity was compared with an ANCOVA analysis. Parametric tests were carried out after successfully fulfilling requirements of normal distribution (Lilliefors test, p<0.01) and homogeneity of variance (Levene’s tests). Inter-population comparisons with significant differences (p-value< 0.005) were a posteriori compared with Tukey’s HSD test (Sokal and Rohlf, 1995). The numerical relationships of fecundity with environmental pH conditions were explored with Pearson’s correlation tests, such as plastic responses which were expressed as a mean reaction norm that had a significantly non-zero slope (Stearns, 1992). Such a mechanistic approach to detect the relationship between phenotypic plasticity and specific environmental features can yield insights into ecological speciation among populations (Chen and Hare, 2008). Acclimation experiments were carried out with three different cohorts of each population. Although these cohorts were sampled in subsequent chronologies, such sampling was not continuous and thus, assessed cohorts correspond to discrete subgroups within each population (group). To evaluate inter-population differences in the EPR during acclimation experiments, we conducted a linear mixed effects model. Extensions of simple linear models and mixed effects models allow the simultaneous assessment of fixed and random effects on non-independent data sets, such as phenotypic/genetic variations across cohorts of a given population (i.e., hierarchical structure). Experimental conditions (temperature, salinity, and body size) and populations (estuarine and upwelling) were considered as fixed effects, while experiments (I, II, III) as a random effect. Differences in fecundity plasticity were indicated by a significant effect of population. Statistical analyses were performed in PRIMER6+ and STATISTICA package10.

A total of 61 mtCOI sequences, 29 from the upwelling and 32 from the estuarine population were aligned in a region of 508 bp (Genbank access code OQ877133-OQ877193). COI sequences showed a moderate to high diversity in both localities with 13 haplotypes in total. Particularly, a relatively higher diversity was revealed in the estuarine with respect to the upwelling population (details in Table 2 ), with the presence of 6 private haplotypes. In the upwelling population, we found 3 private haplotypes ( Figure 1 ). The most abundant haplotypes were present in both populations, while the private haplotypes were less frequent ( Figure 1A ). However, the pairwise Φst comparison showed a considerable and significant genetic divergence between populations (Φst=0,414; p-value < 0.001). An alignment of 1460 bp to nuclear 18S was constructed with 19 sequences, 12 from the estuarine and 7 from the upwelling population (Genbank access code OQ875872-OQ875890). These sequences showed a high diversity in both populations with the presence of 12 haplotypes ( Figure 1B ). As with mtCOI, estuarine individuals showed a greater diversity than upwelling copepods ( Table 2 ), with the presence of 8 haplotypes. No shared haplotypes were found, all being private haplotypes. However, a low and not significant genetic divergence was found between populations (Φst =-0.008; p-value= 0.374), likely related to high nucleotide differentiation (diversity) among estuarine individuals (π=0.133 and see detailed pairwise comparisons in Supplementary Table 1 .

There was a marked divergence between inter-daily temperature variation assessed for the estuarine (11.9–16.9°C) and upwelling (14.0–17.3°C) populations ( Figures 2A, B ). Temperatures below 13°C were significantly more prevalent (>70%) in the estuarine system ( Table 3 ). Unimodal (13.1°C) distribution of temperature at the estuarine system contrasted with warmer (>15°C), unimodal distribution observed at the upwelling system. A clear divergence in salinity ranges and variability was evidenced between populations ( Figures 2C, D ). The highest salinity values (33.7 psu) experienced by copepods in the estuarine system were below the lower threshold (34.7 psu) of salinity variations experienced by upwelling copepods. Mean ( ± SD) pH_T values, 7.94 ± 0.14 (estuary) and 7.92 ± 0.12 (upwelling), overlap between both coastal habitats ( Figures 2E, F ). Differences in mean pH conditions were not significant between A. tonsa populations ( Table 3 ). The magnitude of extreme low pH events was similar for both populations (pH_T <7.7), although the frequency of such events was three times higher (9.4%) for the estuarine than upwelling (3.6%) population. Food in the estuarine system was higher and less variable (190 ± 76 µg C L^-1) than in the upwelling (175 ± 200 µg C L^-1) site ( Figures 2G, H ), such that, there was a significant difference in the food conditions between systems. Mean (± SD) values of environmental conditions experienced by each cohort of both populations are in Supplementary Figure 2 .

The distribution of female body size tended to overlap to a certain extent between populations ( Figures 3A, B ). The prevalence (>20%) of individuals larger than the mode (>1.1 mm) in estuarine females were significantly higher than upwelling congeners (p-value; Table 3 ). Within EPR ranges extended up to 70 egg ind^-1 d^-1, the mode in both systems were below 20 egg ind^-1 d^-1. Relatively smaller copepod females were found in the upwelling habitat, which produced larger brood-sizes (13 egg fem^-1 d^-1) than those produced by the largest females collected in the estuarine system (11 egg fem^-1 d^-1) (p-value; Figures 3C, D ). There was a relationship between body size and fecundity (ANCOVA, p<0.05) in the estuarine population. Similar ranges of egg size (13–15 µm) were found in both populations ( Figures 3E, F ), though lower egg size values (<82 mm) were significantly more prevalent (>70%) in the estuarine population ( Table 3 ).

There were population-specific patterns of phenotypic plasticity in fecundity to field pH conditions ( Figure 4 ). Such pH variations explained 58% (estuary, y =-543 + 70*x; p<0.0001; r^{2 =} 0.58) and 47% (upwelling, y=-259 + 34*x; p<0.0001; r^{2 =} 0.47) of the EPR variability according to Pearson’s correlation coefficient, while student’s t-test_0.05 indicate significant differences (p=0.002) between population EPR slopes. Accordingly, a positive non-zero EPR slope was observed in both populations with pH increasing ( Figure 4 ), although the sensitivity (i.e., EPR change rate) of estuarine copepods was twice as high as that of the upwelling counterparts.

Relatively similar in magnitude, extreme low pH levels were three times more frequent in the estuarine than in the upwelling system. Natural extreme (pH_T 7.68–7.78) and high (pH_T 7.94–8.06) pH levels were replicated in three 96-h acclimation experiments involving females from both populations. There were significant differences in temperature (F_{1,24 =} 1488, p=0.0001) and salinity (F_{1,24 =} 249, p=0.0001) levels between experiments conducted with both populations ( Supplementary Table 2 ). In addition, the temperature (F_{2, 24 =} 11, p=0.0004) of experiment I and salinity (F_{2, 24 =} 6, p=0.006) of experiment III with the estuarine population were significantly lower than other treatments. Both control and extreme pH scenarios were similar among the experiments for both populations (F_{1,24 =} 1.1, p=0.3). The phenotypic plasticity in EPR to extreme pH conditions varied among the three different cohorts assessed for each A. tonsa population ( Figure 5 ). The full linear mixed effects model indicates reproductive outcomes were significantly influenced by salinity, body size, pH-condition, population, and the interaction experiment*pH-condition*population (adj. R^{2 =} 0.78) ( Table 4 ). The effect of body size did not appear in a best-fitting model, while its adjusted R² was 0.68.

The analysis of molecular diversity based on the mtCOI and nuclear 18S genes including the number of haplotypes and polymorphic sites, haplotype, and nucleotide diversity ( Table 2 ), indicates the estuarine population was genetically more diverse and possibly older than that present in the upwelling habitat. The considerable level of genetic variation also suggests that both populations have not faced recent events of drastic population reductions (bottlenecks), which could be caused by strong selection pressure (Lloyd et al., 2016; Gurgel et al., 2020). Furthermore, the observed molecular divergence based on the mtCOI gene, and the absence of shared haplotypes to 18S between estuarine and upwelling populations, suggest genetic structuration and limited genetic connectivity between examined A. tonsa populations. Cryptic lineages have been observed in A. tonsa across latitudinal and environmental gradients (Caudill and Bucklin, 2004; Chen and Hare, 2008), sometimes coexisting in sympatry and manifesting variable phenotypic responses (i.e., reaction norms) to specific habitat drivers (Chen and Hare, 2011).

Environmental data ( Figure 2 ) well represented those environmental drivers, for example, temperature, salinity, and pH conditions, prevailing either seasonally during the reproductive period of the estuarine population (Aguilera et al., 2013; Garcés-Vargas et al., 2020; Osma et al., 2020) or daily and synoptic variations in the upwelling system due to a combination of solar heating and surface mixing (Kaplan et al., 2003; Piñones et al., 2007). Fluctuations in pH conditions, including extreme low events, in the assessed coastal habitats are related to river discharges (Aguilera et al., 2013; Pérez et al., 2016) or wind-driven upwelling (Torres et al., 2002; Vargas et al., 2017; Aguilera et al., 2020a), both characterized by cold water. The significant temperature-pH_T relationship ( Supplementary Figure 3 ) provides confidence that these events were well represented in the current study. However, we also acknowledge some potential caveats. At a mean temperature between 15–16°C, observed in both examined systems, A. tonsa individuals can complete its ontogenetic development within 17–18 days (Miller et al., 1977), with a mean life span around 35 days (Ceballos and Kiørboe, 2011; Kiørboe et al., 2015). This suggests our sampling approach might have largely underestimated high frequency (hours-day) variability affecting environmental pH in both coastal systems (see for example Hofmann et al., 2011; Carstensen and Duarte, 2019; Kapsenberg and Cyronak, 2019). Such environmental variability might be critical in shaping the phenotypic plasticity of short life cycle copepods (Gaitán-Espitia et al., 2017; Vargas et al., 2017). For example, day-night pH cycles have been observed in estuaries (Baumann and Smith, 2018) and upwelling (Saderne et al., 2013) ecosystems associated with photosynthesis/respiration balance. In this sense, assessed ranges of pH_T in both coastal systems overlapped to a certain extent despite contrasting oceanographic, climatic, and geographic forcing operating in distant eco-geographic provinces (Vargas et al., 2017; Vargas et al., 2022). However, Pearson’s analysis indicated population-specific phenotypic responses to local pH variations ( Figure 4 ) that might reflect physiological adequations to local scale pH conditions (Fitzer et al., 2012; Aguilera and Bednaršek, 2022). Contradictory results have been found with regards to pH effects (null, positive, and negative) on Acartia species through in situ (Hansen et al., 2019; Aguilera et al., 2020b) or laboratory observations by means of short (day-week) to long term (transgenerational) experiments (Langer et al., 2019; Dam et al., 2021). This emphasized that both approaches are needed to draw more comprehensive patterns in the response of copepod species and related biological communities to climate change (Reum et al., 2016). In the case of copepod physiological responses to low pH, such interaction can be modulated if ample and nutritious food is available (Aguilera et al., 2020b; Cominassi et al., 2020), such as that observed in both coastal habitats. Significant correlations between EPR and field pH conditions might also correspond to the effect of other simultaneous environmental drivers affecting copepods physiology, including changes in temperature (Sasaki and Dam, 2019; Dam et al., 2021) and food spectra (Kleppel and Burkart, 1995; Jónasdóttir et al., 2009). For example, river-discharges can introduce organic carbon and nutrients from different land use and the resultant physical-chemical conditions might influence the tolerance of marine organisms to environmental perturbations associated with climate change (Pérez et al., 2016). However, the concentration of nutrients and carbon in the Valdivia estuary are relatively low due to its low human intervention and greater vegetation coverage (Pérez et al., 2015).

An additional characterization of EPR phenotypic variations was accomplished in this study through sequential acclimation experiments (e.g., Stillman, 2003), in which mean habitat conditions (temperature, salinity, pH and food condition) of each population were replicated. Temperature and salinity levels varied among experiments conducted with the estuarine population ( Supplementary Table 2 ), likely affecting phenotypic outcomes of the local population. Since experimental temperature deviations (0.5°C) were still within the most prevalent thermal conditions in that habitat ( Figure 2A ), and salinity differences (0.2 psu) were well below changes able to impact the metabolic scope of these copepods (Calliari et al., 2008), observed experimental variations might have a relatively low influence on observed phenotypic responses. This is supported by results of the linear mixed effects model that indicated no significant effect of temperature on EPR variation ( Table 4 ). However, body size, which is linked to habitat temperature (Miller et al., 1977), did influence EPR during acclimation experiments according to the full linear mixed effect model ( Table 4 ). Multiple discrete cohorts inherently contain genetic information and manifest phenotypic plasticity in response to specific local environmental conditions (Sasaki and Dam, 2020). Differences in female age and related effects on fecundity (Rodríguez-Graña et al., 2010), not visually detectable during 96-h acclimation experiments, could have reduced the influence of body size on experimental EPR since this effect did not appear on a second best-fitting model.

During acclimation experiments, copepods from both populations were fed at similar food concentration (400 µg C L^-1) but different composition. Estuarine A. tonsa females were fed unialgal Isochrysis sp. Diet, whereas upwelling copepods a mixed diet composed by Isochrysis sp. And Tetraselmis suecica (3:1 proportion). Subtle (< 30%) chemical (lipids, carbohydrates, and proteins) differences have been found between different Isochrysis and Tetraselmis strains (da Silva Gorgonio et al., 2013; Custódio et al., 2014). However, such differences seem to not lead to relevant changes in copepods reproduction (Lee et al., 2006), likely due to Isochrysis being the main source of nutritional compounds (Knuckey et al., 2005), which are highly required during invertebrates’ reproduction (Müller-Navarra et al., 2000). Acclimation experiments considered 96-h to acclimate copepod EPR to experimental conditions. Copepod ingestion might require at least 48 h to acclimate to laboratory non-limiting food conditions while digestive enzymes dampened short-term food variations (Mayzaud et al., 1992). Thus, acclimation to food (ingestion, digestion) might have occurred within the first 48–72 h of our experiment. Considering A. tonsa can convert ingested food into egg production within less than 10 h (Tester and Turner, 1990), it is highly possible our acclimation period might have been enough to evaluate copepods physiological outcomes (acclimation) during our experiment.

According to our GLM analysis, the population*pH-condition did not significantly affect EPR during acclimation experiments ( Figure 5 ). Although observed phenotypic variation tended to be higher in the estuarine population in agreement with genetic diversity and plastic variations to field pH conditions, only a single cohort (CIII) of this population exhibited a significantly different performance during acclimation experiments. Environmental conditions observed during the sampling of the CIII in the estuarine habitat were significantly less variable than that affecting other cohorts of the local population ( Supplementary Figure 2 ). The environmental conditions, involving mean values, fluctuations and extreme events, experienced seasonally across the development of the cohort can induce changes in the phenotypic plasticity and tolerance to environmental perturbations in Acartia species (Sasaki and Dam, 2019). Such a plasticity does not necessarily imply adaptive plasticity nor local adaptation. In both cases, transgenerational experiments are required to remove environmental influence on and elucidate adaptive plasticity and genetic components of phenotypic plasticity (Brennan et al., 2022). In order to assess the extent to which this phenotypic variation and plasticity are influenced by genetic components, potentially linked to local adaptation, a different experimental approach will be required using common garden conditions.