Isolation of strains from coastal Chile, the Juan Fernandez Archipelago, and oceanic waters west of Peru between 2011 and 2013 was previously described (Bendif et al., 2016; von Dassow et al., 2018). In November-December 2015, strains were isolated from 3 stations during the LowpHOx 1 cruise aboard the R/V Cabo de Hornos. To be able to compare to previous strains collected from near the Punta Lengua de Vaca upwelling at 30°S, water was collected from the onboard continuous seawater system as the ship passed this latitude, and temperature and salinity were recorded from the ships thermosalinometer. At two other stations, water was sampled from 5 m from the CTD rosette equipped with 10L Niskin bottles. In August 2016 a one-day field campaign was conducted on the R/V Stella Maris II from Coquimbo to south of Punta Lengua de Vaca, and water was collected from beside the boat using a torpedo system connected to a CTD pump as previously described. Following previously described protocols, water samples (100 ml) were concentrated to 1–2 ml by gentle centrifugation and coccolith-bearing cells were individually isolated based on depolarization of forward scatter light by an InFlux Mariner Cell Sorter flow cytometer (von Dassow et al., 2012, 2018; Bendif et al., 2016). In the LowpHOx 1 cruise, the flow cytometer was onboard and samples were processed immediately. In 2016, samples were carried in coolers to the lab in Santiago for processing the next day. Details of strain origins are provided in Supplementary File 1.

Once established, strains were identified by scanning electron microscopy (SEM) to species-level (E. huxleyi or Gephyrocapsa species) and morphotype (E. huxleyi) and classified (von Dassow et al., 2018). Strains were then maintained in L1 medium in an environmental chamber on a 12:12 light-dark cycle where temperature fluctuated from 14 (lights off) to 15°C (lights on). Cultures were maintained by transfer to fresh medium every 3–4 weeks.

Cultures were collected by filtration or centrifugation and DNA was extracted using either a DNeasy Plant Mini kit (Qiagen) or following a modified CTAB protocol: Cells pelleted by centrifugation were extracted in 700 μl of 100 mM Tris-HCl (pH 7.5), 25 mM EDTA, 1.5 M NaCl, 2% (w/v) CTAB, and 0.3% (v/v) β-mercaptoethanol, ground with acid-washed glass beads, and 3 μl of proteinase K (20 mg ml^–1) was added. Subsequent steps followed a protocol for recalcitrant plants (Healey et al., 2014). PCR and partial sequencing of cox1, cox3, and tufA were performed as described previously (Bendif et al., 2016). Additionally, partial cox1 and cox3 sequences were downloaded from Genbank and those that could be aligned over sections of 1375 (cox1) or 810 bp (cox3) were kept for analysis. Sequences were aligned with ClustalW. Haplotype networks were constructed individually for all three markers, as well as on the concatenated cox1–cox3 alignment, by median joining (Bandelt et al., 1999) using PopART (Leigh and Bryant, 2015). Genetic groups were determined for all three alignments using Automatic Barcode Gap Discovery (ABGD) (Puillandre et al., 2012). Gene trees were constructed using Maximum Likelihood using RaxMLGUI 2.0 (Bootstrap = 1000, GTR+I+G) on CIPRES, and in parallel with Bayesian Inference on Mr. Bayes (HKY+I+G, 50M of generations, sample frequency every 1000, temperature 0.30, burning 25%).

For most strains for which cox1 and/or cox3 sequences are available in Genbank, location and month of isolation, but not associated environmental data, is available as meta-data either from the source literature or from public data provided by culture collections (see Supplementary Files 1, 2). Reconstructed ocean environmental parameters sea surface temperature (SST), sea surface salinity (SSS) and mixed layer depth (MLD) was assigned to all strains isolated after 1993 (from which cox1 and/or cox3 sequences were used here) using the Multi Observation Global Ocean ARMOR3D L4 analysis and multi-year reprocessing (Guinehut et al., 2012; Mulet et al., 2012) actualized from the Copernicus Marine Services website¹ on 20 Feb 2020. Minimum distance to mainland coast of strain origins was estimated using Google Earth. In that analysis, the mainland coast was the nearest continent or coast of large islands (>100,000 km², e.g., Iceland or Japan). For each strain origin location and date, the monthly mean SST, SSS, and MLD was extracted. Multinomial logistic regression was performed using the nnet package (Rippley and Venebles, 2021) in R (R Core Team, 2020). For strains isolated from the Southeast Pacific by the lead author or lab members for which direct measurements of SST and SSS are available (Supplementary File 2), an additional analysis was performed using these data.

Thermal reaction norms here are based on fitting Thermal Performance Curves (TPCs) for growth. Two parallel TPC analyses were performed. Recent comparative work analyzing or re-analyzing published TPC data in diverse phytoplankton (Thomas et al., 2012; Boyd et al., 2013; Anderson and Rynearson, 2020) have fit such data to the Norberg equation (Norberg, 2004), a general function for a temperature response for growth based on the Eppley curve, an envelope function proposed to describe an evolutionary limit to maximum growth rate vs. temperature (Eppley, 1972). Such an approach has also been used for investigating asexual evolution of thermal reaction norms in clonal populations of E. huxleyi subject to selection in the lab at different temperatures (Listmann et al., 2016). For comparison with this literature base, it is our primary analysis. The general Norberg equation gives growth rate μ_i as a function of temperature T for individual or species i as:

(1) μ i ⁢ ( T ) = a i ⁢ e b ⁢ T ⁢ ⌊ 1 - ( T - Z i W i / 2 ) 2 ⌋

where the trait Z_i is the temperature at which growth approaches the Eppley curve optimum, and W_i is the thermal niche width (in °C). Thus Z_i is not the temperature of maximal growth but the temperature at which maximal growth is closest to the maximum exhibited by phytoplankton generally at that temperature, a represents a competitive optimum. Norberg’s derivation fixed the Eppley curve coefficient a at 0.59 and the exponent b at 0.0633 (Norberg, 2004). Both a and b have been considered by other authors to be free parameters to fit in published analysis of experimental TPC curves (Thomas et al., 2012; Boyd et al., 2013; Listmann et al., 2016). We returned to the original interpretations of Norberg and Eppley (Eppley, 1972; Norberg, 2004), that the envelope represents a constant shape defining the observed upper limits of phytoplankton performance with temperature across all species, and kept b fixed at 0.0633 while letting a_i vary among strains. In this way, a_i has a simple biological interpretation: a_i/0.59 is the degree to which strain i reaches the envelope maximum at T = Z_i. In fact, for many strains it was not possible to obtain adequate fits, or reasonable values of Z_i, when b was free. Also, Thomas et al. (2012) used b = 0.0631 as a theoretical value, but we found slightly better fits keeping b = 0.0633 (not shown). Curve fitting for the Norberg equation was performed in Prism 9 (GraphPad Software). We note that previous studies have derived T_opt,i, the temperature of maximum growth, from the Norberg equation and reported this. We were interested in comparing the both thermal niche width W_i and T_opt,i. T_opt,i was derived numerically in R from both both Z_i and W_i, so is not independent, and only Z_i and W_i were used for statistical testing.

A key previous study on E. huxleyi thermal reaction norms (Zhang et al., 2014) instead used the Template Mode of Variation (TMV) approach, so we also performed TMV analysis in order to compare differences among groups most directly to what was observed in that study. The TMV approach also has the advantage of accommodating more complex TPC shapes (Izem and Kingsolver, 2005). This involves finding a common curve shape f(x) to the individuals to be compared, and then for each individual i an individual curve is fitted:

(2) μ i ⁢ ( T ) = 1 w i ⁢ f ⁢ ( T - m i w i ) + h i

Here, m_i is the temperature of maximal growth (from here on referred to as T_opt), h_i is the height parameter (average growth rate), and w_i is the non-dimensional width parameter. Here a polynomial of degree 4 was fit for f(x), following Zhang et al. (2014). The fit was performed in Matlab using the scripts provided by Izem and Kingsolver (2005). A large number of strains (most E. huxleyi, all G. oceanica, and most G. parvula) grew positively at 27°C and it was necessary to include this point to ensure a value of μ_i(T) at a temperature higher than the temperature of maximum growth, while no G. muellerae grew at 24°C and so could not be tested at 27°C using the progressive acclimation approach. The TMV scripts require at least one data point at all temperatures tested, so G. muellerae had to be analyzed separately. By visual inspection, we detected isolated cases where the TMV analysis appeared to strongly over-estimate max growth rate (and therefore the height parameter h_i) compared to the Norberg curve fits. This might relate in part to the fact that we could not measure negative growth rates. The empirical minimum growth rate was constrained at 0 but the fit polynomial in TMV analysis could go below 0, which would affect the weighted sum-of-squares error in the fitting procedure. Therefore, we used the TMV analysis only to compare differences in T_opt observed here to those reported by Zhang et al. (2014).

For E. huxleyi haplogroups I and II in the first experiment and all groups of E. huxleyi strains in the second experiment, the number of strains was sufficient to allow normality testing. The distributions of a, Z and W from the Norberg fits and the w, h, and T_opt from TMV fit were all consistent with normality according to the Kolmogorov-Smirnov and Shapiro-Wilkox tests. For these parameters, significance of differences among groups was evaluated using 1-way ANOVA tests with Tukey’s multiple comparison test to identify significant pairwise differences. The distribution of w from TMV fits did not pass normality tests for E. huxleyi haplogroup I from either experiment, so this parameter was compared using the Kruskal-Wallis test with Dunn’s test for determining significances of pairwise differences. These tests were performed in Prism 9 (GraphPad Software).

A haplotype network was constructed for the concatenated cox1–cox3 alignments of 186 strains for which both mitochondrial markers were available (Figure 1), a number which includes 163 strains newly sequenced here (119 new strains from 2015 and 2016 campaigns, 44 strains isolated by the same methods between 2011 and 2013), and 23 sequences from Genbank. ABGD analysis supported three or four separate haplogroups (see Supplementary Figures 1–3). The three group separations corresponded to the clades I, II and III previously defined (Hagino et al., 2011; Bendif et al., 2014, 2016; although termed α, β and γ, respectively, in the latter two studies, we revert to the roman numeral designation to avoid confusion with morphotype classifications A, B, and C, and adopt the term haplogroup). No haplogroup III (G. oceanica) strains were isolated by this study in the Southeast Pacific, but both haplogroups previously documented for E. huxleyi were present. However, with the four group separation proposed by ABGD, clade I was separated into two haplogroups, with haplogroup Ib composed only of strains isolated from the Southeast Pacific. Both Maximum Likelihood and Bayesian trees (Figure 2) agreed with the placement of haplogroup Ib as basal to haplogroup I, but with only moderate support (77% bootstrap support and 75% posterior probability).

Sequences were available only for the shorter cox3 from 324 strains, while only cox1 was available for 17 strains in Genbank and 3 strains sequenced here. To extend the study to include these 344 strains, haplotype networks for cox1 and cox3 were analyzed separately (Supplementary Figures 4–9). ABGD analysis assigned the four G. oceanica strains for which cox1 sequences were available (clade III) to two distinct haplogroups. More importantly, cox1 and cox3 haplotype networks were congruent in the assignation into haplogroup I composed of G. parvula, and G. ericsonii, and some strains of E. huxleyi, and haplogroup II composed of G. muellerae and the other strains of E. huxleyi. Additionally, although ABGD analysis did not provide support for separating haplogroup Ib from I based on the shorter cox1 or cox3 sequences alone, strains could be assigned to cox1–cox3 haplogroup Ib based on the haplotypes assigned from networks built with either cox1 or cox3 separately, so there was no evidence for haplogroup 1b outside of the 412 strains from the Southeast Pacific.

A tufA haplotype network was constructed from the alignment of 172 total sequences, including 161 sequences from our Southeast Pacific isolates and 11 sequences that could be included in the alignment that were obtained from Genbank from other strains. The tufA network was divided into 4 groups by ABGD analysis (Figure 3 and Supplementary Figures 10–12). One group was populated by sequences from only two E. huxleyi strains from the Southeast (SEPA25 and SEPA584) and another group by a single sequence of G. oceanica (RCC1316). Within the strains from the Southeast Pacific, there was a very strong association between mitochondrial and chloroplast haplogroups (c² = 184.8, df = 4, p < 0.0001) (Table 2). Only a single cox haplogroup Ib isolate was assigned to tufA haplogroup II rather than I (of 48 sequenced for both mitochondrial and chloroplast markers), only a single cox haplogroup II isolate was assigned to tufA haplogroup I (out of 29). Thus, for the rest of the analysis we focus only on cox haplogroups.

Morphological classification was obtained or previously available for a total of 181 strains isolated from the Southeast Pacific between 2011 and 2015. Examples of each morphotype from 2015 to 2016 are shown in Supplementary Figure 13 and morphological analysis identifying SEPA565 as a new isolate of G. muellerae by comparison to Bendif et al. (2015, 2016) is provided in Supplementary Table 1. There was a significant difference among haplogroups in the distribution of morphotypes (c² = 135.5, df = 2, p < 0.0001) (Table 3). Substantial numbers of all three haplogroups were in the broad B morphotype class, which here may include B, B/C, and O morphotypes that were not distinguished here. Almost all of the haplogroup Ib strains were B morphotype. The A-OC morphotype included substantial numbers of both haplogroups I than II. The two morphotypes with over or hypercalcification (closed central area or fusion of distal shield elements) showed contrasting distributions among haplogroups: The A-CC morphotype, with closed central area but separated distal shield elements, was almost exclusively haplogroup I, while the HC/R morphotype (with fused distal shield elements) was almost exclusively haplogroup II. Nevertheless, there was no case where a morphotype was exclusively populated by one mitochondrial haplogroup.

At a global scale, haplogroup I (without separating off Ib, detected only in Chilean waters) was the most widespread (Figure 4A). Haplogroup II became more prevalent at higher latitudes in both hemispheres. All 10 strains originating from latitude >55° were haplogroup II (5 North Sea, 1 Baltic Sea, 2 northern central Atlantic, 1 northern central Pacific, 1 northwestern Pacific), whereas haplogroup II was mostly excluded from latitudes less than 30° N or S. Haplogroup III (G. oceanica) strains came from sub-tropical and tropical latitudes.

This latitudinal pattern held within the regional scale in the Southeast Pacific (Figure 4B) but appeared to be modified in relationship to patterns of upwelling. Haplogroup II dominated in the cooler nearshore waters from the strong coastal upwelling center north of Punta Lavapie (37°–36° S) to the strong upwelling center near Punta Lengua de Vaca (30° S), but was represented only 9 out of 33 strains isolated from the oceanic waters near the Juan Fernandez Archipelago (33° S, 78° W), and only 4 of 10 strains from the northern Chiloe outer coast at 42° S, south of the region where coastal upwelling is dominant.

There was a significant negative association of haplogroup II (vs. haplogroup I) with increasing SST (Figure 4C), SSS and minimum distance from mainland coast (MDMC), while MLD was not significant (Tables 4, 5). In contrast, haplogroup Ib was not significantly associated with SST but was significantly associated negatively with SSS and MDMC and positively with MLD, but this group was only found in a small subset of locations in coastal Chile and nowhere else. Haplogroup III (vs. haplogroup I) was positively correlated with SST and SSS and negatively with MLD and MDMC, but the p-value for SST was just above the a-threshold for significance (p = 0.055). Also, SST was significantly correlated with SSS and MDMC.

When only strains from the SE Pacific isolated after 2011 were considered, haplogroup II (vs. haplogroup I) was significantly negatively associated with both SST and SSS from direct measurements of the surface water of origin (Figure 4D and Tables 6, 7).

The mitochondrial phylogenies generated in this study, both of cox1 and cox3 sequences separately and of concatenated cox1–3, were largely consistent with expectations from previous published mitochondrial phylogenies of the Emiliania-Gephyrocapsa complex (Hagino et al., 2011; Bendif et al., 2014, 2016). However, with the inclusion of a larger number of strains we detected a new mitochondrial sub-clade not observed in previously studies. We assigned this new clade as haplogroup Ib. Its phylogenetic association with haplogroup I was supported by the ABGD analysis of concatenated cox1–cox3 sequences in network analysis, which divided mitochondrial sequences into either three or four clades. When cox1–cox3 sequences were divided into only three clades, haplogroup Ib was within the clade of haplogroup I. Trees constructed by both Maximum Likelihood and Bayesian Inference also placed Ib as basal to haplogroup I, although with only moderate support, and the posterior probability for this placement in Bayesian Inference was low (0.75). Haplogroup Ib has not yet been observed in any strains isolated outside of Chilean coastal waters, and all but one of the Ib strains were isolated during or soon after the strong 2015–2016 El Niño event, when Ib was very common among strains isolated from 30° S and 27° S.

The culture-based approach here also allowed comparing mitochondrial and chloroplast phylogenies in a substantial set of strains. We were expecting to find a more complex relationship of mitochondrial (cox) and chloroplast (tufA) haplogroups based on the incongruence reported between organellar phylogenies (Bendif et al., 2016). Distinct mitochondrial haplogroups of E. huxleyi are expected to be able to interbreed or at least to have done so in the recent past, based on a phylogenomic study suggesting extensive mixing among nuclear genomes (Filatov, 2019), as well as the occurrence of strains where mitochondrial and chloroplast haplogroups contrast. However, such inter-organelle phylogenetic incongruence turned out to be rare in the E. huxleyi strains from the SE Pacific analyzed here. Patterns of organellar inheritance are not known in these organisms, so the pattern might imply inheritance of both organelles is preferentially from one parent following syngamy, or perhaps also a degree of post-zygotic cytonuclear incompatibility (Willett, 2011) could be important.

The fact that mitochondrial and chloroplast haplogroups mostly coincided permitted important simplifications in analysis with respect to environmental data and in experimental design. For example, as no strains exhibited cox haplogroup I with tufA haplogroup II, and only 1 strain (less than 1%) exhibited cox haplogroup II with tufA haplogroup I (Table 2), it was not possible or necessary to include such combinations in the experimental determination of thermal reaction norms.

There was only a partial correspondence between morphotype and haplogroup. All haplogroups exhibited morphotypes classed in the broader B morphotype class, although we did not attempt here to distinguish among B, B/C, C, and O in this class. Distinct A morphotypes are distinguished more easily, based on whether the distal shield central area is uncovered (the original A morphotype), or covered/closed by over-calcification (here termed A-CC) or whether the distal shield elements are fused (von Dassow et al., 2018; Díaz-Rosas et al., 2021). This last is considered to be a variant of the R morphotype (Beaufort et al., 2011), but which we term HC/R to indicate that the strains from this region have a tendency for variable overgrowth of the central area (von Dassow et al., 2018; Díaz-Rosas et al., 2021), a character distinct from the original description of the R morphotype (Young et al., 2003). All of these morphotypes contained both haplogroups I and II. The original study defining mitochondrial haplotypes in E. huxleyi found the R morphotypes only in haplogroup I (Hagino et al., 2011), but the HC/R morphotypes from Chilean waters were preferentially group II. Morphotype is likely nuclear-determined, and is tightly associated with the nuclear GPA gene that produces a protein related to coccolith-associated polysaccharides (Schroeder et al., 2005), however, morphotype and associated genetic polymorphisms in GPA do not correspond to microsatellite population genetics (Krueger-Hadfield et al., 2014), suggesting different morphotypes may interbreed. The imperfect correspondence between morphotypes and mitochondrial haplogroups observed here thus might also support the possibility of nuclear genetic exchange among E. huxleyi organelle haplogroups.

The global latitudinal pattern previously reported in E. huxleyi mitochondrial haplogroups was also reflected at the regional scale. The association of E. huxleyi II with decreasing SST, both on a global scale (using reconstructed SST) and on the regional and local scale (using direct SST measures), were consistent with temperature explaining the observed latitudinal patterns. However, the multinomial regression also detected significant associations with SSS and with distance from the coast in the global analysis. A caution is that salinity also shows a strong latitudinal gradient globally. Also, a large number of the E. huxleyi I cox1 or cox3 sequences are from the high salinity Mediterranean, but only one Mediterranean E. huxleyi strain was assigned as haplogroup II, and the strain origins that are farthest from continental shelves are in tropical and subtropical waters (with higher temperatures and salinities). Nevertheless, a statistically significant relationship with salinity was also suggested at the local scale, from direct measurements of salinity among samples collected near 30° S. Indeed, strain-specific variability in reaction norms with respect to salinity has been demonstrated (Gebühr et al., 2021). Temperature and salinity are the conservative tracers of water masses, and so other physiologically and ecologically important variables, such as nutrients or total chlorophyll, also tend to co-vary with temperature, salinity, and distance from the coast. This emphasizes the value in experimentally determining whether physiological reaction norms in fact are consistent with patterns of specific environmental parameters at ecological scales.

The environmental associations of the new haplogroup Ib are not yet clear. So far, it has only been identified from coastal waters of central Chile. It was almost exclusively observed during and immediately after an intense El Niño, when the northern Chilean coast typically experiences a tendency for deepening of the thermocline and warmer, more saline waters (Escribano et al., 2004; Aguilera et al., 2020). Measured SST at 30° S in Nov. 2015, when haplogroup Ib was most important, was 15.8°C, 1.7°C higher than the temperature in Nov. 2012, and 3.2°C higher than the average measured in Oct. 2011 (Supplementary File 2), but it was also prevalent at the end of Aug. 2016 when cooler temperatures were present, and was not found to the north in warmer waters. Continued sampling effort, both regionally and globally and using both cultures and culture-independent methods, would be necessary to determine the geographic distribution and oceanographic preferences of this clade.

Thermal reaction norms (or at least thermal performance curves for growth, when norms were not calculated) have previously been reported for two species studied here, including a total of 17 strains of E. huxleyi (Conte et al., 1998; Buitenhuis et al., 2008; Zhang et al., 2014; Listmann et al., 2016; Wang et al., 2019) and two strains of G. oceanica (Conte et al., 1998; Buitenhuis et al., 2008). The previous studies that included G. oceanica showed it to be less cold-tolerant than E. huxleyi, but did not permit estimating T_opt as determined here. In the present study, E. huxleyi strains exhibited similar thermal reaction norms to those previously reported for this species, with the exception of one strain from the SW Indian ocean previously reported to show relative cold-intolerance and two strains from subpolar waters reported to show a much lower T_opt, both by Conte et al. (1998). To these species we add data from G. muellerae and G. parvula/ericsonii, two species that were only recently isolated into culture.

Hypothesis 1, the possible role for ocean temperature in determining distributions among the different species in the Emiliania-Gephyrocapsa genus, was supported by the thermal reaction norms measured in the laboratory. The species restricted to cooler waters, G. muellerae, very clearly had lower temperature thermal reaction norms, both in terms of optimal temperatures and thermal limits for growth (although one of the three strains barely grew at 8°C). Consistent with their geographic distributions, the six strains of G. oceanica tended to show the highest optimal temperatures and thermal limits, closely followed by the five G. parvula/ericsonii strains. Finally, the norms for E. huxleyi were well above the norms for G. muellerae and overlapped with the lower end of the ranges exhibited by G. parvula/ericsonii and G. oceanica. Although the differences in thermal optima (both Norberg Z, the competitive optima, and T_opt) between G. parvula/ericsonii and E. huxleyi were not significant, all E. huxleyi tested in both experiments could grow at 8°C whereas none of the G. parvula/ericsonii could grow at that temperature, supporting that G. parvula/ericsonii is less cold-tolerant.

The competitive optimum temperature Z related linearly with the SST of origin among the Emiliania/Gephyrocapsa species. However, the slope was less than 1. This may reflect the larger trend among phytoplankton to exhibit higher optimum temperatures for growth than their temperature of origin except at the highest temperatures (Thomas et al., 2012). We caution that the pattern observed might be modified if there is adaptation in the lab to temperatures at which cultures are maintained prior to experiments. Nevertheless, we did not observe trends of Z or difference of Z and SST of origin vs. time in culture either within each E. huxleyi haplogroup or among the G. oceanica strains. The recent study of Anderson and Rynearson (2020) did not find a relationship between thermal optima for growth and isolation temperature among isolates of the diatom genus Skeletonema, despite a wider range of isolation temperatures (0 to 20°C). That difference might partly be explained by the fact that Skeletonema is a neritic species, and that earlier study focused on strains isolated from the Narraganset Bay. Such neritic organisms may be retained within in the bay and exposed through the year to much greater changes in temperature than coccolithophores in the Emiliania/Gephyrocapsa genus with a more oceanic habitat and no known benthic resting stages.

To avoid local adaptation complicating the ability to detect differences among haplogroups, all but two strains tested in the Experiment 2 were from the same water sample (and the observed differences held when those two strains were excluded). Even in the same water sample, cells may have different histories, particularly in zones of intense upwelling where there is mixing of different water masses. We minimized this effect by focusing on a water sample that was well away from strong upwelling centers. The three strongest centers of upwelling along the Chilean coast are at 37° S, 30° S, and 23° S, with several weaker upwelling centers between 37° S and 30° S and at 29° S (Thiel et al., 2007), and the site at 27°S was therefore well away from these. We also attempted to minimize the possible effects of in-culture evolution on thermal reaction norms as TPCs were fully completed within 18 months of isolation for all the E huxleyi strains used in the second test.

The results presented weak support for Hypothesis 2: Among E. huxleyi haplogroups, the thermal reaction norms of haplogroups I and II exhibited significant differences of the expected directions only in the strains isolated in the same 2 week period and tested within 18 months of isolation. In this set, haplogroup II tended to show cold-shifted thermal reaction norms compared to haplogroup I, whether in terms of temperature of maximum growth T_opt, the competitive optimal temperature (Norberg’s Z), or in growth at 27°C (but not in growth at 8°C).

Despite the phylogenetic association of E. huxleyi haplogroup Ib with haplogroup I, its thermal reaction norms were cold-shifted to be similar to haplogroup II in the second experiment. Although the environmental distribution of Ib remains to be defined, the narrower (more specialist) thermal reaction norms of this haplogroup suggests a more specialist type, which might explain its more sporadic occurrence and apparently more restricted distribution compared to E. huxleyi I and II.

The differences among E. huxleyi I, Ib and II in Experiment 2 were relatively small, on the order of 1°C. The differences were similar to the intra-specific differences in thermal optima reported in planktonic diatoms such as Skeletonema marinoi, S. pseudocostatum, and Thalassiosira rotula (Boyd et al., 2013; Anderson and Rynearson, 2020), but are much smaller than those reported within T. pseudonana or the dinoflagellate Akashiwo sanguinea (Boyd et al., 2013). In the study of E. huxleyi strains isolated from Bergen and the Azores, a 22 degree latitudinal difference, corresponding to a 6.3°C difference in maximum monthly mean temperature or a 9.6°C difference in minimum monthly temperatures, the mean difference among populations was only 1.1°C (Zhang et al., 2014). A higher difference was reported for strains from the Indian Ocean compared to subpolar North Atlantic and North Pacific (Conte et al., 1998). Clonal selection under asexual reproduction in the lab was reported to result in a 0.7°C shift in thermal optima in E. huxleyi within 1200 asexual generations (Listmann et al., 2016). In that study the selection temperature used, 26°C, was near the thermal limits and well outside the range of the environment from which that original strain was isolated (Bergen). Similarly, in the diatom T. pseudonana, selection at 31°C (vs. 16°C) was reported to drive a nearly 2°C change in T_opt within 450 asexual generations (O’Donnell et al., 2018). Both of those evolution experiments imposed a very strong and constant selective pressure, both in terms of an extreme selection temperature and for maximizing growth rates (high r-selection) compared to that expected in nature, where temperatures are lower and selection for maximum growth rate (over other traits) is unlikely to be constant. Nevertheless, this emphasizes how small were the differences determined in this study among E. huxleyi haplogroups.

The lack of difference between E. huxleyi haplogroups I and II in the first experiment is unlikely to only reflect in-culture evolution, even though the strains used had been in culture for 3–4 years by the time the preliminary experiments had started, and up to nearly 6 years when completed. Any in-culture evolution was not enough to have abolished the larger differences observed among species, where it was necessary to compare strains isolated in many different years (and including with time in other culture collections), and the single haplogroup Ib strain from 2011 that was included in the first test showed lower temperature thermal reaction norms, similar to Ib strains from 2015. Also, Conte et al. (1998) observed larger differences among five E. huxleyi despite many years in culture. At regional scales, temperature and salinity correlate as conservative properties of water masses (for the case of the Southeast Pacific, see, e.g., Silva et al., 2009), which also can determine or be correlated with variables such as nutrients or mixed layer depth, which may in turn influence productivity and biomass. It is possible that perhaps adaptation to these other variables, or the interaction of one or more these variables with temperature, might also influence E. huxleyi mitochondrial haplogroup distributions. For example, variations among E. huxleyi strains have been seen in reaction norms to salinity (Paasche, 2001; Gebühr et al., 2021) and in utilization of different N-sources (Strom and Bright, 2009).

Another intriguing possibility is that thermal adaptation is important, but there is nuclear gene flow between mitochondrial haplogroups I and II within E. huxleyi, as discussed above. Although the mitochondrial genome might contribute to thermal adaptation in other eukaryotes (Baker et al., 2019), the results suggest that the dominant control is nuclear-encoded in Emiliania-Gephyrocapsa. In that case, the present latitudinal distributions of the haplogroups might represent a historical association between organelle haplogroups and thermal adaptations, perhaps dating from earlier introgression events between E. huxleyi and other Gephyrocapsa species as hypothesized by Bendif et al. (2016). In regions where the haplogroups overlap, exchange of nuclear genes among E. huxleyi as well as local adaptation might then be erasing the association of organellar history with thermal traits.

Generalists are species which can occupy broad niches (diverse habitats) while specialists occupy more narrow niches, which may reflect wider or narrower physiological tolerance ranges, respectively (Futuyma and Moreno, 1988; Sriswasdi et al., 2017). Based on the wide geographical distribution of E. huxleyi (Paasche, 2001; Winter et al., 2014) and the much more restricted distributions of its closest relatives (Gephyrocapsa) in the present ocean (McIntyre et al., 1970; Bollmann and Klaas, 2008), E. huxleyi would be considered a generalist while the Gephyrocapsa may be specialists. Although it is only one of the multiple abiotic and biotic dimensions over which niches are defined (Colwell and Rangel, 2009), temperature is a dominant variable controlling biochemistry and physiology so is expected to be fundamental in the distributions of organisms, ranging from microorganisms to ectothermic vertebrates (Eppley, 1972; Angilletta et al., 2003; Norberg, 2004; Logan and Cox, 2020). One of the principal classes of evolutionary trade-offs considered to underly species distributions is the generalist-vs.-specialist trade-off, and it has long been expected that this trade-off should be apparent as an inverse relationship between the height and width of thermal performance curves (Angilletta et al., 2003; Izem and Kingsolver, 2005; Buckley and Kingsolver, 2021). The width of the thermal niche has been speculated to play a key role in phytoplankton distributions (Boyd et al., 2013).

Some support was provided for Hypothesis 3 (wider thermal niche of E. huxleyi explaining its broader latitudinal distribution): Haplogroup I of E. huxleyi did exhibit a broader thermal (fundamental) niche in the laboratory in comparison to the Gephyrocapsa species and has the widest latitudinal range compared to the other E. huxleyi haplogroups and the other species. Nevertheless, the newly detected mitochondrial clade Ib of E. huxleyi exhibited a narrower thermal niche in both experiments (although only a single strain was available in the first experiment), and so far it has been identified only between 30 °S and 27 °S on the Chilean coast. Also, although all the E. huxleyi strains reported here (as well as in most other studies) grew at 8°C, Conte et al. (1998) reported that a strain from the SW Indian Ocean could not grow at 9°C, while strains from the subpolar Pacific and from Norway had an apparent T_opt below 20°C and lower thermal tolerance compared to other E. huxleyi. That shows that some E. huxleyi strains can exhibit thermal reaction norms outside the ranges observed here and in the other studies cited, so regional thermal adaptation also contributes to the broad distribution of E. huxleyi. Likewise, some E. huxleyi, such as haplogroup Ib, do exhibit narrower thermal niche breadths. Thus, the broad latitudinal distribution of E. huxleyi can partly be explained by such thermal reaction norm variation within the species (Hypothesis 4).

There was no significant tendency for strains or species with narrower thermal performance curves (or narrower tolerances) to have higher maximum growth rates, either as absolute maximum growth rates (in units of 1/time), or normalized to the Eppley Curve to adjust for the fact that the highest possible maximum growth rates are expected to be higher at higher temperatures in the range of interest here (Eppley, 1972). Previous studies comparing phytoplankton thermal reaction norms across broad phylogenetic groups have also not detected such a relationship (Thomas et al., 2012; Boyd et al., 2013; Anderson and Rynearson, 2020). In part, it may be difficult to estimate well the width parameters (or difference between critical minimum and maximum temperatures), as discussed by Boyd et al. (2013). Nevertheless, as several studies have not detected such a tradeoff, temperature and thermal reaction norms alone may not be sufficient to explain wider or narrower species distributions in phytoplankton. Physiologically, interactions among temperature and other parameters such as salinity, light, nutrients or CO₂ might be important (Feng et al., 2008). Analysis of distributions of E. huxleyi morphotypes and related Gephyrocapsa species in the Southeast Pacific showed that realized niches could be wide with respect to some parameters (e.g., pH and pCO₂) but narrow with respect to others (e.g., temperature, salinity, and calcite saturation state) (Díaz-Rosas et al., 2021). Variable temperatures can negatively affect growth rate (Wang et al., 2019) and might be more important as a driver of adaptation in nature. Thus, although we did not detect a generalist-specialist trade-off in lab tests of reaction norms to a single parameter, such a trade-off in a multidimensional niche space might still be necessary to explain why there were times and locations where specialists appeared dominant, e.g., E. huxleyi haplogroup Ib at 27° S in 2015 (this study) or G. parvula at 20.7° S in 2013 (Bendif et al., 2016).

As the global ocean warms, the habitat thermally permissive for G. oceanica and G. parvula/ericsonii is expected to expand poleward. In the case of G. oceanica, its thermal optima are near or even below the SST at the warmer end of its range, so the low latitude end of its range should also move poleward. Gephyrocapsa parvula/ericsonii may partially follow the border of the South Pacific Gyre (McIntyre et al., 1970; Bollmann and Klaas, 2008), and this gyre is becoming larger (Polovina et al., 2008). In the case of G. muellerae, Eastern boundaries of ocean basins may experience intensification of upwelling-favorable coastal winds, resulting in lower temperatures (Sydeman et al., 2014), a phenomenon documented to be occurring in the Southeast Pacific (e.g., Schneider et al., 2017). This latter process might permit G. muellerae to maintain or even extend its range equatorward, confined to the narrow coastal zone with cooler upwelling water, while this trend continues. Nevertheless, in these waters, G. muellerae appears to exist in lower abundance compared to E. huxleyi (Díaz-Rosas et al., 2021). Overall, E. huxleyi behaves as a generalist and is capable of rapid adaptation, and the lack of evidence for strong differentiation in thermal niches between haplogroups suggest that gene flow may be high among E. huxleyi when populations with different adaptations mix. Thus, E. huxleyi might be expected to increase its dominance over its close relatives in the changing ocean.