Adaptation to environmental temperature in divergent clades of the nematode Pristionchus pacificus

Abstract

Because of ongoing climate change, populations of organisms are being subjected to stressful temperatures more often. This is especially problematic for ectothermic organisms, which are likely to be more sensitive to changes in temperature. Therefore, we need to know if ectotherms have adapted to environmental temperature and, if so, what are the evolutionary mechanisms behind such adaptation. Here, we use the nematode Pristionchus pacificus as a case study to investigate thermal adaptation on the Indian Ocean island of La Réunion, which experiences a range of temperatures from coast to summit. We study the evolution of high-temperature tolerance by constructing a phylogenetic tree of strains collected from many different thermal niches. We show that populations of P. pacificus at low altitudes have higher fertility at warmer temperatures. Most likely, this phenotype has arisen recently and at least twice independently, consistent with parallel evolution. We also studied low-temperature tolerance and showed that populations from high altitudes have increased their fertility at cooler temperatures. Together, these data indicate that P. pacificus strains on La Réunion are subject to divergent selection, adapting to hot and cold niches at the coast and summit of the volcano. Precisely defining these thermal niches provides essential information for models that predict the impact of future climate change on these populations.

Adaptation, divergent selection, parallel evolution, Pristionchus pacificus, temperature

Temperature is an important environmental variable, which affects the physiology of organisms (Huey and Stevenson 1979), and extremes of hot and cold adversely affect fertility (Walsh et al. 2019), ultimately resulting in sterility. The range of temperatures over which an organism remains fertile tends to match the environmental temperature of its niche (Addo-Bediako et al. 2000). Species are therefore rarely exposed to extreme temperatures in their optimal habitat. However, during colonization of a new habitat, an organism might be exposed to stressful temperatures more often, resulting in a decrease in its fitness and potential loss of fertility if it is unable to use behavioral responses to regulate its body temperature indirectly or to acclimate physiologically (Angilletta et al. 2002). If organisms cannot acclimate, then they might adapt to a warmer habitat by restoring their fertility at higher temperatures—a process referred to as thermal adaptation (Porcelli et al. 2016). Studying thermal adaptation is particularly important in the light of increasing evidence for climate change (Alexander et al. 2013), yet its underlying evolutionary mechanisms are poorly understood.

A useful framework to understand temperature adaptation is the thermal performance curve, which measures the response of fitness traits to changes in temperature. The optimum temperature of an organism, as well as the critical limits of fitness (i.e., the upper and lower limits of the range of temperatures where the organism remains fit), is typically shown (Angilletta 2009). In this context, temperature adaptation is the evolutionary process by which an organism adapts its performance curve to match the temperature profile of its habitat. When a species’ habitat experiences a range of temperatures, evolutionary pressure can result in isolated organisms adapting their performance curves differently. For example, the release of stabilizing selection can result in a broadening of the thermal performance curve (generalists) and directional or divergent selection can give rise to organisms with narrower curves and higher or lower optimum temperatures (specialists, Knies, et al. 2009). In the latter case, divergent selection might result in reproductive isolation if continued thermal adaptation gives rise to two subspecies that have nonoverlapping thermal preferences and inhabit niches with nonoverlapping temperature profiles (Keller and Seehausen 2012).

There is some evidence that thermal adaptation occurs in nematodes and flies (Trotta et al. 2006; Prasad et al. 2010). However, it is not known to what degree organisms adapt to the temperature of their niche by parallel or convergent evolution (Arendt and Reznick 2008). Examples of parallel and convergent evolution have been shown for cichlid fish, which have adopted similar body forms while adapting to the same ecological niche in the lakes of East Africa (Kocher et al. 1993; Elmer et al. 2010). Alternatively, changes in the temperature tolerance of organisms could be influenced by other evolutionary mechanisms, like the founder effect or drift, if temperature is not a strong selective factor (Lynch and Walsh 2007).

Thermal adaptation occurs in organisms like nematodes, many of which are amenable to laboratory analysis and therefore represent excellent model organisms for studying this phenomenon. As ectothermic organisms (those that adopt the temperature of their environment), their physiology is strongly affected by temperature, which also affects their fertility (Petrella 2014; Begasse et al. 2015; Leaver et al. 2016). For example, there is evidence that the nematode Caenorhabditis elegans and its relative Caenorhabditis briggsae have adapted to different temperature niches, the former being more prevalent in temperate regions and the latter being more prevalent in the tropics (Prasad et al. 2010 and references therein). Similarly, we have previously shown that there is temperature-dependent variation in the fertility of Pristionchus pacificus natural isolates (Leaver et al. 2016).

Pristionchus pacificus offers an opportunity to understand thermal adaptation, because much is known about its biology and ecology. It has a short life cycle of about 3.5 days that includes progressive development through four juvenile larval stages (Hong and Sommer 2006). Adults can produce ∼175 offspring at optimum temperatures, with this trait known to vary with temperature (Leaver et al. 2016). Pristionchus pacificus is a soil nematode that is reliably found in associations with scarab beetles, in a necromenic relationship—nematodes can be found on beetles where they remain in a dormant state until the death of the host, at which point they begin feeding on the bacteria growing on the decomposing beetle (Herrmann et al. 2010; Ragsdale et al. 2015). Knowledge of this relationship has enabled the collection of large numbers of natural isolates, that is, more than 40 Pristionchus species and more than 1000 P. pacificus strains (Rödelsperger et al. 2017, 2018). A rich source of P. pacificus isolates is the island of La Réunion in the Indian Ocean (Morgan et al. 2012; McGaughran et al. 2016). This island is an active shield volcano (Gillot and Nativel 2002) and its combination of climatic variables and complex geology has generated a diverse set of habitat types (McGaughran and Morgan 2015). These abiotic features make La Réunion an ideal case study of thermal adaptation because different nematode populations experience a range of temperatures, from hotter coastal regions to cooler locations at the volcano peak at approximately 3000 m a.s.l. (Jumaux et al. 2011). Thus, La Réunion offers a natural experimental template to test if different populations of P. pacificus have differentially adapted to local temperature regimes.

The evolutionary history of P. pacificus on La Réunion and its host interactions are very well understood. La Réunion has been colonized many times independently by representatives of the major clades of P. pacificus and their beetle hosts (Herrmann et al. 2010; Morgan et al. 2012; McGaughran et al. 2013). Originally, four clades of P. pacificus were distinguished using mitochondrial markers, designated clades A–D (Morgan et al. 2012). Subsequently, whole genome SNP markers further subdivided clade A into three distinct clades—called A1, A2, and A3—and showed that clades C and D are closely related—redefining them as C1 and C2, respectively (Rödelsperger et al. 2014, 2017). Four clades are predominantly found on La Réunion, namely, the SNP clades A2, B, C1, and C2. Strains from these clades are more prevalent at particular locations and on particular beetle hosts. For example, clade C1 strains are found spread across most of the western part of the island (Fig. 1a; Morgan et al., 2012) with the rhinoceros beetle Oryctes borbonicus representing a major host (Meyer et al. 2016). Maladera affinis beetles are host to clade A2 strains, whereas clade B strains form an isolated population at high altitudes and are associated with the beetle Amneidus godefroyi (Herrmann et al. 2010; Meyer et al. 2016; Moreno et al. 2016). Competition between recently arrived and well-established beetle hosts has presumably limited nematode dispersal, as suggested by the finding that the four clades have largely distinct distributions that only overlap occasionally. Since colonizing the island, as these populations adapt to local conditions they continue to accumulate genetic and phenotypic variations, for example, by recombination via outcrossing that happens within clades and, more rarely, between clades (Morgan et al. 2012; Rödelsperger et al. 2014) and a minor degree of host switching (McGaughran and Sommer 2013, McGaughran and Morgan 2015). This and other lines of evidence (Morgan et al. 2012; Rödelsperger et al. 2014) suggest that P. pacificus is undergoing a process of incipient speciation (Mayer et al. 2015; McGaughran et al. 2016). Although members of the different clades are cross-fertile, they are on the way to becoming different species, making it an ideal system to study if different lineages have adapted to different conditions on the island.

Figure 1

The occurrence of the Htt strains increases with decreasing altitude of their collection site. (a) The location of sampling sites on La Réunion. Also shown are the locations of the weather stations (Jumaux et al. 2011) and the peaks of the volcanoes (above approximately 2000 m altitude) and the east/west divide (Morgan et al. 2012). (b) The correlation between altitude and the temperature recorded by weather stations (red line: linear fit with 95% confidence intervals).

In this work, P. pacificus natural isolates from La Réunion were used to study thermal adaptation in a model organism amenable to laboratory study. To achieve this, fertility was used as a proxy for fitness. The tolerance of P. pacificus to temperature extremes was assayed by measuring the number of offspring produced at low temperatures (10°C); scoring the ability of animals to give rise to fertile offspring at high temperatures (30°C); and determining the highest temperature at which strains could remain fertile. We used these data to test whether temperature tolerance is subject to natural selection by assaying these temperature tolerance phenotypes for hundreds of strains collected from different sites on the island that experience different temperatures. We also tested whether different populations have adapted to hot and cold temperature niches, accounting in our analysis for the influence of phylogenetic history and host beetle of each strain. The evolution of temperature tolerance was studied in a detailed phylogeny made from whole-genome data and by using thermal performance curves as a framework to understand the evolutionary mechanisms behind thermal adaptation.

Results

POPULATIONS OF NEMATODES ON LA RÉUNION ARE SUBJECT TO NATURAL SELECTION

We define the high-temperature tolerant (Htt) phenotype as the ability to give rise to fertile offspring at 30°C (Leaver et al. 2016). If P. pacificus is adapting to temperature on La Réunion, then the Htt phenotype should be subject to natural section. To test this hypothesis, we looked for a relationship between the frequency of the trait and the selective factor (Endler 1986). Here, we used altitude as the selective factor, which is a good proxy for the maximum annual temperature on La Réunion. To do this, we used the existing collection of natural isolates of P. pacificus from La Réunion and collected additional strains from sites below 500 m (Fig. 1a and see Methods), thus sampling from sea level to the volcano plateau. These locations cover a range of maximum daily temperatures from 22 to 36°C (Jumaux et al. 2011), with strains encompassing all of the genomic clades that are found on the island (McGaughran et al. 2016).

Natural isolates were assayed for fertility at high temperature in the following way. Strains were maintained at 20°C. Next, individual larva of the third juvenile stage of development (J3) were placed onto fresh Escherichia coli-inoculated plates, which were then shifted to 30°C (methods). After 7 days of incubation at the test temperature, worms were inspected to determine if they had given rise to offspring. Strains were scored as Htt if they had given rise to offspring that were themselves fertile. Strains were scored as tolerating a mesothermal temperature (Mtt) if they gave no offspring or gave rise to only one generation (Methods). We had previously defined this class as low-temperature tolerant (Ltt; Leaver et al. 2016), but we now re-designate this term in the light of new data (see below). In total, 19.7% of the 310 tested strains were Htt (Table S1).

As expected for a trait subject to natural selection, low-altitude locations tended to have more Htt strains. For example, 75.9% (n = 29) of strains from Saint Benoit (SB; 22 m altitude) were Htt, whereas no strains from locations above 1000 m a.s.l. were Htt (n = 74) (Table S1). To quantify this relationship, we looked at the maximum annual temperature recorded by weather stations at different locations on the island (published by the Bureau d'étude climatologique; Jumaux et al. 2011), as it reflects the highest temperature nematodes might experience over the year. A linear model between maximum annual temperature and the associated altitude of the weather stations provided a good fit to the data (linear regression model: Alt ∼ Temp, R² = 0.9336, F_1,23 = 324, P = 4.81 × 10^–15; Fig. 1b). Thus, maximum annual temperature was used to estimate the temperature of all nematode collection sites (Fig. 1b). Next, the occurrence of the Htt phenotype was plotted by altitude with maximum annual temperature as a secondary axis (Fig. 2a and see Methods for details). To do this, the data were partitioned into bins (see Methods) and the number of Htt strains per bin alongside the mean altitude for all strains in that bin was calculated. The percentage of strains giving rise to fertile offspring at 30°C was as high as 70% at the lowest altitude bin and only 6.7% for the bin centered at 897 m (Fig. 2a). For all bins above 1000 m a.s.l., the percentage of Htt strains was zero but, for strains below this altitude, the data fell onto a straight line (linear regression model: Alt ∼ Htt, R² = 0.6038, F_1,6 = 9.14, P = 0.0233). Thus, altitude is related to high-temperature tolerance. Analysis of variance (ANOVA, Methods) showed that variation in Htt was significantly apportioned across the altitudes at which strains were collected (ANOVA model for binary data: Htt ∼ alt, D² = 74.13, df = 1, residual deviance = 233.3, null deviance = 307.5, P = 2.2 × 10^–16). Because of the linear relationship between altitude and maximum annual temperature, these findings show that the occurrence of the Htt phenotype is related to environmental temperature (Fig. 1b).

Figure 2

The percentage of Htt strains differs depending on their collection altitude (as a proxy for temperature), beetle host, and clade. (a) The correlation between altitude and the occurrence of Htt strains. Strains were ranked by the altitude of their collection site then grouped into bins for all strains. The percentage of Htt strains in each bin was calculated and plotted against the average altitude of the sample site of all strains in that bin (Methods). The secondary axis is the maximum annual temperature at the corresponding altitude shown for convenience (as extrapolated from the linear model presented in Fig. 1b). The blue line is a fit for the bins centered between 0 and 1000 m. The red line is the temperature where the percentage of Htt strains reduced to zero (30.8°C). (b) Chart showing the range of altitudes where beetles were found for each beetle species (point and line represent the mean and the 5th to the 95th percentile of the data). (c) The correlation between altitude and the occurrence of Htt for strains in bins corresponding to their beetle host. Again, the maximum annual temperature at the corresponding altitude is shown on the secondary axis. For nematodes isolated from Adoretus sp., the data were split up into bins of 34 strains and fit with a straight line (yellow) and compared to the fit for the whole dataset (blue). (d) Bar chart showing the percentage of nematodes that are Htt for strains isolated from each beetle species. The dashed red line is frequency of Htt strains for the whole dataset. (e) The correlation between altitude and the occurrence of Htt for strains in bins corresponding to their clade (clade c and d were split into bins of 30 or 33 strains) with the maximum annual temperature at the corresponding altitude shown on the secondary axis. (f) Bar chart showing the percentage of nematodes that are Htt for strains belonging to each clade, plotted as in panel d.

Using the linear model, we predicted that the altitude where the percentage of Htt strains reduced to zero should be around 877 m (Fig. 2a). Interestingly, the maximum annual temperature at this altitude is 30.8°C, as estimated from the weather station data (Fig. 1b,a), which is very close to what we would expect from our phenotype assay at 30°C.

STRAINS ISOLATED FROM BEETLES FOUND AT LOWER ALTITUDES ARE MORE LIKELY TO BE Htt

Next, we investigated other factors that might influence the Htt phenotype. The main extrinsic factor for which we have systematic data is the beetle species from which the nematodes were isolated. These beetles have different habitat ranges corresponding to their own preferred altitudes (Fig. 2b). For example, A. godefroyi is found only at higher altitudes (approximately >2100 m a.s.l.; Moreno et al. 2016), whereas other species, including M. affinis, prefer lower altitudes (approximately <300 m a.s.l.; Ahrens 2003). We note that the temperature of the beetle is unlikely to influence the nematode as the beetles themselves are also ectothermic and adopt the surrounding temperature of the soil or air. Any heat generated by the beetle as it decays is likely to quickly dissipate. Considering these factors together, we argue that altitude (as a proxy for environmental temperature) is the variable that should influence both the distribution of beetle species and the percentage of Htt strains. Thus, the occurrence of Htt strains isolated from different beetle species might also differ. Indeed, Htt strains were found more frequently on M. affinis and Aphodius sublividus (found at low and intermediate altitudes 0–1425 m), but none were found on A. godefroyi and O. borbonicus (found at higher altitudes approximately 1350–2350 m a.s.l.; Fig. 2b,d). ANOVA also showed that the beetle species has a significant effect on the percentage of Htt strains (ANOVA model for binary data: Htt ∼ beetle, D² = 62.53, df = 8, residual deviance = 244.5, null deviance = 307.0, P = 1.5 × 10^–10).

To eliminate any influence that different beetle hosts might have on the occurrence of Htt strains, we plotted the data for strains isolated from different beetles separately (Fig. 2c). The largest number of nematode strains was isolated from Adoretus sp. (n = 200/310), which is found across a wide range of altitudes. Indeed, Adoretus-derived strains show a negative effect of altitude on the percentage of Htt strains (linear regression model: Alt ∼ Htt(Ado), F_1,4 = 5.26, P = 0.0834, R² = 0.5682; Fig. 2c, dashed line). The magnitude of the estimated effect among the Adoretus-derived strains is stronger than that estimated among the whole dataset, but we found no statistically significant difference between group (ANCOVA: Htt∼Alt+group vs. Htt∼Alt*group, F_1,10 = 0.181, P = 0.6163), although the correlation is weaker (see R² values above) meaning the relationship between altitude and the percentage of Htt strains is maintained in this subset of strains. Thus, different beetle species carry different frequencies of Htt strains and, more importantly, for strains from Adoretus, the correlation between altitude (as a proxy for environmental temperature) and the occurrence of Htt is robust.

STRAINS THAT BELONG TO CLADE A ARE MORE LIKELY TO BE Htt

Next, we considered the evolutionary relationship between strains of P. pacificus as another potential influencing factor on the frequency of Htt strains. For example, strains belonging to clade B might respond differently to temperature because they form an isolated population (Moreno et al. 2016) that might be susceptible to founder effects or nonadaptive forces, such as drift. Therefore, the percentage of Htt strains belonging to each of the four mitochondrial clades was calculated. Figure 1f shows that no clade B strains were Htt. In contrast, clade A strains, which are more prevalent at lower altitudes, have a higher percentage of Htt strains. ANOVA indicates that the percentage of Htt strains found in different clades deviated from that which would be expected purely by chance (ANOVA model for binary data: Htt ∼ clade, D² = 67.82, df = 4, residual deviance = 239.6, null deviance = 307.5, P = 6.6 × 10^–14). Thus, we conclude that a strain will be more or less likely to be Htt depending on which clade it belongs to.

TEMPERATURE IS THE MAIN FACTOR INFLUENCING THE FREQUENCY OF Htt

Finally, we attempted to distinguish which of the variables—beetle, clade, or altitude—best explained the occurrence of Htt. The influence of beetle can be eliminated by analyzing strains isolated from Adoretus sp. (Fig. 2c), but the influence of clade and altitude remain to be disentangled. To address this, we performed model selection analysis using Bayesian Information Criterion (BIC). BIC is a statistical tool that selects from a set of models the one that best explains the observed data using a likelihood function (Burnham 2004). The log-likelihood of each of the general linear models (Htt ∼ Alt, Htt ∼ beetle, Htt ∼ clade) was calculated and converted to weighted differences. The BIC analysis showed that altitude explained the variation in frequency of the Htt trait better than beetle or clade (0.99 for altitude, 7.8 × 10^–6 for clade, and 7.3 × 10^–12 for beetle: the largest value indicates the variable that best fits the data, see Methods for details). This result supports our hypothesis that altitude (as a proxy for temperature) is likely the main factor influencing the frequency of Htt phenotype.

PHYLOGENY OF WORLDWIDE STRAINS OF P. pacificus SHOWS MULTIPLE INSTANCES OF THE ACQUISITION OF THE Htt PHENOTYPE

To investigate the evolution of temperature tolerance in a larger collection of P. pacificus strains, we combined the current La Réunion dataset with strains from worldwide locations that we had analyzed previously (Leaver et al. 2016). Single-nucleotide polymorphisms (SNPs) were identified by comparing whole-genome data of these strains, and used to construct a maximum likelihood phylogenetic tree with ultrafast bootstrapping (Minh et al. 2013; Nguyen et al. 2015). The worldwide diversity of P. pacificus consists of six major clades based on genome-wide markers: clades A1, A2, A3, B, C1, and C2 consistent with previous work (Rödelsperger et al. 2014; Baskaran and Rödelsperger 2015; Leaver et al. 2016).

When the phenotypic data were overlaid on the phylogenetic tree, it showed that Htt strains are sporadically scattered throughout clade C1 (Fig. 3). This finding suggests recent and independent acquisitions of high-temperature tolerance in strains from La Réunion and Mauritius, as had been previously suggested for a smaller dataset (Leaver et al. 2016). Also, there were Htt strains in clades A1, A2, and A3, including subclades with many Htt strains or single strains from diverse locations (Fig. 3). To test if temperature tolerance has arisen multiple times independently, we inferred the phenotype of ancestral strains using a Bayesian algorithm. This analysis predicted that the last ancestor of all P. pacificus strains was likely to be Ltt (indicated at the nodes of the major branches of the tree; Fig. 4). Furthermore, the ancestor of deeply-rooted clade B is predicted to have been Mtt and the same is true for clades C1, C2, and A1 (Fig. 3). In contrast, the ancestor of clade A2 is predicted to have been Htt (Fig. 3). Thus, the Htt phenotype is likely a derived state and could have arisen independently as clade A2 diverged, and more recently within clades A1, A3, and C1.

Figure 3

Phylogenetic tree of P. pacificus strains shows one clade highly enriched with Htt strains. A maximum likelihood tree based on 873,932 genome-wide SNPs from 212 strains. The phenotype of each strain is indicated by a circle at the terminal node of the tree (red = Htt, green = Mtt). The sample site is indicated by a colored square. For the nodes of the major clades, the inferred ancestral phenotype is indicated by pie charts and the degree of bootstrap support shown at the internal nodes.

Figure 4

Increased fitness at high temperature comes at the cost of decreased fitness at low temperature. (a) Number of offspring for representative strains from clade A2 from low altitude, clade C1 from intermediate altitudes, and clade B from high altitude. Inset shows diagram of island with cold and hot niches; dashed line delineates the upper limit of the hot niche at 877 m. (b) Maximum thermal fertility limit (TF_max) in the same strains.

SWITCHING OF THE Htt PHENOTYPE BACK TO THE ANCESTRAL STATE

A closer examination of strains from clade A2 shows that there are many closely related strains that are Htt and form a “hot” subclade (Figs. 3, S1). Most of these strains were collected from Saint Benoit on La Réunion, which has an altitude of 22 m and experiences a maximum annual temperature of 34°C. Interestingly, there is one example of a reversion back to the Mtt ancestral state in this subclade: strain RSC028, which is Mtt despite being very closely related to two Htt strains, RS5417 and RS5424. There are further examples of reversions in clade A2, including an event that gave rise to a “cold” subclade (Fig. S1). Also, two strains in this subclade, JU138 and RS5264 (from Hawaii and Bolivia, respectively), have regained the Htt phenotype.

Using the available mutation rates from P. pacificus mutation accumulation lines, clade C1 strains seem to have adapted to 2°C higher temperatures within 5.6 × 10⁵ generations (Molnar et al. 2011; Weller et al. 2014). In comparison, the reversion and switching of phenotypes observed in clade A2 is likely to have happened in far fewer generations. We prefer to estimate the observed divergence in number of generations rather than absolute time frames because the exact generation time in the wild is hard to estimate: experimental studies showed that P. pacificus generations are likely to range between 3 days (25°C on agar plates in the laboratory) and 1 year (the maximum survival of dormant larvae under laboratory conditions) (Mayer and Sommer 2011).

A LOW TEMPERATURE NICHE AT HIGH ALTITUDES

Given that strains at low altitude had adapted to higher temperatures, we asked whether the opposite might be also true for strains at higher altitude. Therefore, we measured the number of offspring at 10°C (Methods) of a subset of clade B strains, which are found at altitudes higher than approximately 2100 m a.s.l. (Fig. 4a, inset). These strains had a fairly robust number of offspring at 10°C: six strains having between 46 and 133 offspring on average at 10°C, but two had very low fertility. Specifically, RSB001 gave on average only four offspring and RSB013 only two. This observation is counter to expectation, indicating that these two strains are maladapted to low temperature, perhaps because of inbreeding depression in the isolated clade B population (Moreno et al. 2016). However, for the other six strains, fertility was high compared to the P. pacificus reference strain (RS2333), which has on average 47 offspring. Clade C1 strains from intermediate altitudes (Trois Bassin at 1378 m) also had robust fertility at 10°C, giving between 57 and 142 offspring. Both clade B and clade C1 have much higher numbers of offspring than clade A2 strains from low altitudes (Saint Benoit at 22 m), which gave only one to six offspring.

We next measured the fitness at high temperature in the same set of strains by precisely determining the maximum thermal fertility limit (TF_max, the highest temperature where a strain is able to give rise to fertile offspring; Walsh et al. 2019). Interestingly, the data showed that clade B strains had a lower TF_max (either 26 or 27°C). Furthermore, clade C1 and A2 strains had TF_max values of 27–29°C and 29–30°C, respectively (Fig. 3b). The TF_max values of the three clades are statistically different (ANOVA model: TFmax ∼ clade, F_2,18 = 51.2, P = 3.73 × 10^–8). As only a subset of clade A, B, and C1 strains were tested, we cannot be sure that all clade B strains are Ltt or that some Ltt strains are present at intermediate altitudes. Nonetheless, the group of strains tested from the three clades covers a range of TF_max values and the measured TF_max values for strains from clades B and A2 do not overlap.

These two fitness assays show a correlation between fitness at extremes of temperature and the geographical distribution of clades. Clade A2 strains, which are found at low altitudes, have higher fertility at higher temperatures, having a higher TF_max (Fig. 4b). Conversely, clade A2 strains have a lower fertility at lower temperatures, as scored by the number of offspring at 10°C (Fig. 4a). Clade C1, which is found over a broader range of altitudes, also seems to maintain its fertility across a broader range of temperatures, with a robust number of offspring at 10°C (Fig. 4a) and a TF_max in between that of clade A2 and clade B strains. In contrast, clade B strains found at high altitudes have a higher number of offspring at 10°C (Fig. 4a), but a lower TF_max (Fig. 4b). In light of this, we define strains with a TF_max of 26 or 27°C as Ltt, meaning that some strains designated Mtt in the first screen are now redesignated as Ltt (Figs. S2, 3).

Discussion

In this work, we used natural isolates of the nematode P. pacificus to test if they were thermally adapted to the environmental temperature of natural habitat of the Indian Ocean island of La Réunion. We showed that there is a correlation between the altitude at which nematodes were collected and the occurrence of the Htt phenotype. This observation, together with the correlation between altitude and maximum annual temperature recorded at weather stations, led us to conclude that P. pacificus is subject to natural selection by temperature on La Réunion. This is consistent with previous work that showed P. pacificus strains isolated from different altitudes showed temperature-dependent differences in the rate of spontaneous male production, which can also be interpreted as evidence of thermal adaptation (Morgan et al. 2017). This analysis also predicts that there is a niche for P. pacificus Htt strains, corresponding to a zone from sea level up to 877 m a.s.l., where there appears to be a selective advantage to having the ability to give fertile offspring at 30°C.

Phylogenetic analysis of P. pacificus strains suggested that thermal adaptation has happened multiple times independently in the different clades, which is consistent with the phenotype arising many times in parallel. Alternatively, Htt could have evolved once and either remained polymorphic during the diversification of clades or transferred between clades by recombination. We note that the rare occurrence of outcrossing between clades (Rödelsperger et al. 2014) matches the rare occurrence of strains that gain the Htt phenotype, hinting at the possibility that Htt has arisen by adaptive introgression as has been observed in Drosophila (Turissini and Matute 2017). We cannot rule out these possibilities; however, two previous results suggest that they are less probable. First, analysis of mutation accumulation lines indicates that different clades are separated by millions of generations and it seems unlikely that neutral polymorphisms would be maintained for such long periods (Weller et al. 2014). In accordance with this, comparison of nucleotide diversity between clades shows very little (<7%) shared variation (Rödelsperger et al. 2014). Second, analysis of linkage disequilibrium and population structure shows that admixture between different clades seems to be rare (Morgan et al. 2012; Rödelsperger et al. 2014). Therefore, most likely, Htt has evolved multiple times independently from an ancestral strain that was Mtt, in an example of parallel evolution. This most likely happened once in clades A1–3 and once in clade C1. Because we do not have data about the genetic differences underlying this phenotype, we cannot make any conclusions regarding whether the causative mutations in clades A1–3 and clade C1 are in the same gene or genetic pathway, or whether these changes in temperature tolerance may have resulted from completely different mechanisms among clades (i.e., Arendt and Reznick 2008). Elucidating the mechanisms further will be important for future work.

In the phylogeny, we also observed cases of ancestors gaining the Htt phenotype then later losing it, for example, in clade A2. These cases of relatively rapid switching between phenotypes on shorter time scales represent examples of reversion. For example, strain RSC028 seems to have reverted from Htt to Mtt, while two closely related strains, RS5417 and RS5424, appear to have switched from Mtt to Htt. RSC028 was collected at Grand Etang (GE; 527 m) and RS4371 and RS5424 are from Saint Benoit (22 m). This suggests that the ancestor of strain RSC028 was part of an Htt population at Saint Benoit and was likely transported to GE, where it was subjected to lower selective pressure or balancing selection, and then reverted to Mtt. In this case, reversion is likely to have occurred by recombination transferring alleles between individuals of the same clades in the local population (in contrast to the examples of parallel evolution between clades), because recombination rates are high within clades (Morgan et al. 2012; Rödelsperger et al. 2014).

Measuring the altitudinal range at which Htt strains are found, we were able to delineate a high-temperature niche from approximately 0–870 m, to which P. pacificus seems to be adapting. In contrast, we present evidence of a low-temperature niche: clade B strains, known to be geographically restricted to high altitudes (Moreno et al. 2016), are adapted to low temperatures having high fertility at cold temperatures. This comes at the cost of a decreased fitness at high temperature for this lineage. Therefore, the geographic range of clade B strains (approximately 2100 m a.s.l.) also delineates a low-temperature niche.

If we consider the TF_max as a measure of the upper limit of P. pacificus isolates’ thermal fitness curves and the number of offspring at 10°C as a measure of fitness at the lower end of the curve, we can make the following conclusions. For the Htt clade A strains that we tested, the fitness curve is shifted higher and is perhaps narrower than the Mtt clade C strains, typical of “hotter is better” (Frazier et al, 2006; Knies et al. 2009). Clade A Htt strains probably have a lower critical limit around 10°C (because their fertility drops close to zero, indicating that the lower critical limit is close to 10°C) and an upper critical limit of 29–30°C. This compares to Mtt clade C strains, which have high fertility at 10°C—their lower critical temperature is much lower than this temperature and their upper critical limit of 28°C gives a rather broader curve. Such a broad curve suggests that a generalist-specialist trade-off might be at play; however, data on the optimum temperature are needed to confirm this (Huey and Hertz 1984). The Ltt clade B strains that we tested seem to have shifted their fitness curves lower, having both their upper critical limits lower than clade C Mtt strains and a similarly low lower critical limit. Therefore, we propose that thermal adaptation to high- and low-altitude-temperature niches, especially for clades A2 and B, leads to disruptive selection with increased fitness at extreme temperatures. Furthermore, this likely contributes to ongoing allopatric speciation to two geographically isolated temperature niches. However, the range of temperature for which individuals from each clade remain fertile still overlaps, indicating that they are not yet reproductively isolated. In summary, our data suggest that the landscape of La Réunion represents a continuum of temperature niches to which different populations of P. pacificus have adapted. This lays the ground work for future work to fully characterize the thermal fitness curves of strains from these populations/clades.

It is currently unclear if P. pacificus populations on La Réunion can adapt to ongoing climate change quickly enough to persist in their current geographic distribution. Our work provides useful information for models that could predict if this is the case (Jezkova and Wiens 2016). However, further work is needed for such assessments, including measurements of the generation time of P. pacificus in its natural environment to estimate how long it takes to adapt to higher temperatures in years; migration rates to better understand how quickly extant Htt strains might spread to higher altitudes; intraclade recombination in clade A2 to quantify how quickly the rapid switching of the temperature tolerance phenotype occurred; and finally the capacity of P. pacificus to further increase its TF_max under higher selective pressure in a future, hotter climate.

Summary

In a 2015 review, Cutter posed a number of open questions in the field of C. elegans ecology and evolution (Cutter 2015). Specifically, he posed the following outstanding questions for which this work provides answers, at least for P. pacificus clades. Are the primary dimensions of niche space differentiation biotic or abiotic? We identify temperature as an abiotic factor that seems to be in the process of driving P. pacificus to adapt to different geographically isolated niches. Which traits define species range limits? We have shown that the Htt phenotype and the maximum thermal fertility limit are two traits that limit the range of clade B and clade A2. How important in the accumulation of reproductive isolation are divergent adaptation, parallel adaptation to common environments in allopatry, sexual selection, and genomic conflict? We show that disruptive selection (leading to divergent evolution) might be underway in clades A2 and B adapting to high- or very-low-temperature niches and could in the future lead to reproductive isolation between these two populations.

Most excitingly, we can study thermal adaptation in P. pacificus strains that are most likely subject to natural selection in the wild and learn what molecular and cell biological processes are behind this in the laboratory studies.

Methods

COLLECTION OF NATURAL ISOLATES OF P. pacificus

Strains of P. pacificus had previously been collected by the Sommer lab (prefix RS), the Sternberg lab (prefix PS), and Félix lab (JU) from worldwide locations, including many from La Réunion (Tables S2 and S3). To increase the number of strains for analysis, more strains were collected (see Herrmann et al. [2006] for details of collection methods). Beetles were dissected along the midline and placed onto NGM plates seeded with Escherichia coli OP50. After 3 days of incubation, plates were inspected for gravid hermaphrodites using a dissecting microscope. The small subunit ribosomal RNA (SSU) gene was sequenced to confirm strains as P. pacificus.

SELECTION OF SAMPLING AT LOWER ALTITUDES

Preliminary analysis of a set of 149 strains collected from La Réunion by the Sommer lab from 2008 to 2010 showed that one location had a high percentage of Htt strains. However, this dataset was lacking sufficient locations from 28 to 700 m (inhabited by Adoretus sp. beetles) to properly test for a correlation between altitude and the occurrence of Htt strains. We identified locations around Plane d’ Affourches (PD5–PD7) and Saint Phillipe (SP) where Adoretus sp. could be found. Sampling over the subsequent years (2011–2015) collected enough beetles to give nematodes from locations covering the full range of altitudes.

STRAIN MAINTENANCE

Worms were propagated on OP50/NGM plates according to standard methods for C. elegans (Brenner 1974). To minimize epigenetic effects, all worms were kept at 20°C for at least two generations on abundant food before temperature experiments were performed. The incubators used for routine maintenance were Heraeus BK 6160 incubators (Thermo; accuracy: ±0.2°C).

TEMPERATURE ASSAYS

Carefully staged worms of the third juvenile stage of development (J3) from mixed stage plates incubated at 20°C were transferred singly to NGM plates, which were then shifted to 30°C. This stage of development was chosen for the assay because it is a time before the expansion of the gonad (Rudel et al. 2005) so that the temperature shift affects the production of sperm and oocytes. After 7 days of incubation at the test temperature, the plates were inspected using a dissecting microscope and scored for the ability of the founder worm to give rise to fertile offspring. A strain was designated as Htt if the founder worm gave rise to the first filial generation (F1), which subsequently gave rise to the second filial generation (F2)—a situation that can be easily scored because gravid adults (the F1) and a mixed population of eggs and juvenile worms (the F2) can be seen on the plate. A strain was designated Mtt if the founder was killed by heat; arrested during development; developed to adulthood but was sterile; or gave rise only to F₁ that were infertile (i.e., there was no sign of eggs or juvenile worms on the plates). For large-scale application of the Htt assay, a custom-made 30°C room (stability: ±0.5°C) was used to incubate the strains at the test temperature. Any strain that did not give reproducible results was eliminated from the final dataset, giving a final dataset of 310 strains.

The strains RS2333 and RS5194 gave a very robust response to the Htt assay. The result of the assay was reproducible for worms of different ages and for small differences in the assay temperature (±0.5°C). However, in such a large dataset, it is likely that some strains would have an upper temperature limit very close to that of the test temperature of 30°C and could be prone to small changes in the assay conditions. To minimize this error, any strain that did not give reproducible results was eliminated from the final dataset, which was the case for 16 strains. This resulted in a final dataset of 289 strains.

To determine the maximum thermal fertility limit (TF_max) precisely, this assay was repeated for a subset of strains for every temperature between 25 and 31°C. The TF_max was defined as the highest temperature at which a strain still gave rise to fertile offspring.

To determine the number of offspring at 10°C, we used the method to measure the life-time fecundity described by Leaver et al. (2016). Briefly, 10 J3 larva from plates grown at 20°C were picked onto fresh plates, then shifted to 10°C. Founder hermaphrodites were singled to fresh plates each day. After three further days, the number of offspring was counted on each plate and summed to give the total number of offspring for each worm over its life.

We note that we could not measure the number of offspring at 10°C or determine TF_max for all strains in the dataset because these two fitness assays are not as high throughput as the Htt assay. As a result, there are a number of strains for which we assigned the phenotype of Mtt but do not have data on their precise upper temperature limit or fertility at 10°C. However, we are confident in our classification of these strains because low-temperature tolerance seems to be a trait unique to clade B strains that has been shown to be an isolated population that has adopted other phenotypes associated with adaptation to high altitudes, namely, low oxygen (Moreno et al. 2016).

DATA ANALYSIS FOR PLOTTING THE RELATIONSHIP BETWEEN ALTITUDE AND THE PERCENTAGE OF Htt WORMS

Temperature data from weather stations are published in the Atlas Climatique de La Réunion (Jumaux et al., 2011). The annual maximum temperature for each location was plotted against its altitude and fit with a linear model. Then, strains were ranked by the altitude of the collection site and put into bins because, for some sampling locations, low numbers of nematodes were collected. Large numbers of strains are needed to accurately estimate the percentage of Htt strains, especially at intermediate altitudes where the frequency of Htt is low. A bin size of around 30 was chosen to give a good estimate of the frequency of Htt strain per bin. To most evenly distribute the strains across all bins, bin size was set for each condition to n = 30 for the whole dataset, n = 34 for strains from Adoretus sp., and n = 30 for clade C and n = 33 for clade D. For each bin, the mean altitude of the collection site of all the nematodes in that bin was calculated and plotted against the percentage of Htt strains in that bin. The temperature at the corresponding altitude, as extrapolated from the linear model of the weather station data, was plotted as a secondary axis. The map shown in Figure 1a is available under the Open Database License (for further information, visit https://www.openstreetmap.org/copyright).

DATA ANALYSIS FOR DETERMINING THE ALTITUDINAL RANGE FOR BEETLE SPECIES

For each sampling site visited on La Réunion, the number of beetles of each species collected and the GPS coordinates and altitude of each site were routinely noted for years 2010–2015. Strains from worldwide locations were not included in this analysis because we had not previously recorded precise GPS or altitude data as a matter of routine. We define the altitudinal range for a beetle as being the range of altitudes where 95% of beetles were found. Mean altitude of the collection site for all beetles of one species was also calculated.

DATA ANALYSIS FOR DETERMINING THE PERCENTAGE OF Htt STRAINS IN EACH CLADE

To calculate the percentage of Htt strains belonging to each clade, only the mitochondrial clades were used (i.e., clades A–D). We chose to do this because whole genome sequences are not available for all strains. If a strain had a clade assigned from both mitochondrial and SNP data, then the mitochondrial clade was chosen because it is possible to assign SNP clade A1, A2, and A3 strains to mitochondrial clade A. However, the converse is not possible because mitochondrial sequence data are not divergent enough to resolve these clades. SNP clade C1 corresponds to mitochondrial clade D. See column “consensus” in Table S2.

STATISTICAL ANALYSIS

Statistical analysis was performed with RStudio software (RStudio Team 2020) running R version 4.1.1 (R Core Team 2021). The temperature tolerance data were fit with a generalised linear model (binomial distribution) using the glm command. The three variables considered for each strain were the clade it belonged to, the altitude of its collection site, and the beetle species the strain was isolated from. See Table S2 for the full data. Analysis of variance (ANOVA) was performed using the anova command on the generalised linear model for a binomial distribution and a chi-square significance test. ANOVA was assessed for frequency of Htt strains for each variable individually: Htt ∼ Alt, Htt ∼ beetle, Htt ∼ clade. Htt is modeled as binary data, the beetle and clade data are categorical, and the altitude data are continuous. Unless otherwise stated, the analysis was performed on the full dataset with un-binned data. The P-value of the chi-squared test is presented, as well as degrees of freedom (df) and deviance (D²)—the degree of the deviance explained by the model (residual deviance) compared to the null model (null deviance). For ANOVA of FT_max ∼ clade, the anova command was used on the generalised linear model for a Gaussian distribution and the F-test. The P-value and the F-statistic for the degrees of freedom and the residual degrees of freedom (F_x,z) are presented. Analysis of covariance (ANCOVA) was used to test if there is a significant difference between Alt versus %Htt of the whole dataset compared to the subset of data from Adoretus. This was done using the following commands: anova(aov(Htt∼Alt*group); aov(Htt∼Alt +group)), where the dummy variable “group” is either the whole dataset or the Adoretus subset. BIC analysis was performed as follows. The generalised linear model was fit to each variable independently, then for each glm the log-likelihood was calculated using the bic command. The weighted difference between the log-likelihood values was calculated, such that the model that “best fits” the observed data has a value closest to 1. To compare the linear models for the data of percentage of Htt strains versus altitude for the whole dataset and the strains isolated from Adoretus sp. (the only beetle for which we had sufficient sample sizes), the following analysis was performed. A linear regression model was first fitted to the data using the lm command in R using the formula Allalt ∼ Allfreq * groups, where Allalt is the mean altitudes for the binned data for the whole dataset and Allfreq is the percentage of Htt strains for each bin and groups refer to the dummy variable: either the whole dataset or just the data of strains isolated from Adoretus sp. The significance of the interaction term Allfreq:groups was then examined using the summary function to give a P-value (Pr(>|t|)) which showed if the subset and whole dataset were significantly different from each other or not.

CONSTRUCTION OF A PHYLOGENETIC TREE

Whole-genome data for a subset of La Réunion strains and the worldwide (Tables S2 and S3) stains used in the study were already available (McGaughran et al. 2016). Variant bases were identified using a custom pipeline implemented in Perl and Unix as previously described (Leaver et al. 2016), and the genotype of all strains at all variant positions was extracted from the genome sequence covered by all the contigs in P. pacificus draft genome hybrid assembly 1 (Borchert et al. 2010). A total of 873,932 SNPs were selected and concatenated and used to construct a maximum likelihood tree using IQ-Tree 1.5.0 (Nguyen et al. 2015). The substitution rate model GTR+ASC was applied (where ASC corresponds to use of a SNP ascertainment bias model) and the tree was bootstrapped 1000 times using the ultrafast bootstrapping method (Minh et al. 2013).

INFERRED ANCESTRAL STATES

Bayesian analysis was performed using the ace (ancestral character estimation) command (Yang et al. 1995) implemented in the R phytools package (Revell 2011). A single rate model was applied and probabilities for either phenotype for each ancestor were calculated. The data were plotted with the R ggtree package (Yu, 2020) with the probabilities at the internal nodes depicted as pie charts and the figure was assembled in Adobe illustrator.

AUTHOR CONTRIBUTIONS

ML, AM, RJS, and AAH conceptualized the idea of the study. ML and CR designed methodology and performed visualization. ML and AM performed validation. ML, EM, MK, AM, and RJS performed investigation. ML wrote the original draft. ML, EM, AM, CR, RJS, and AAH reviewed and edited the manuscript. ML, RJS, and AAH performed supervision. RJS and AAH acquired funding.

ACKNOWLEDGMENTS

We would like to acknowledge: Katy Morgan, Jan M. Falcke (ne Meyer), Josephine Piffaretti, Jacques Rochat and especially Matthias Herrmann for help conducting fieldwork; Praveen Baskaran for help with statistical analysis; Susanne Ernst, Gabi Eberhardt, Heike Haussmann and Waltraud Röseler for technical support. We would also like the thank Ralf Schnabel and Geraldine Seydoux for their help and support.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

DATA ARCHIVING

All strain and phenotype data are included in the Supporting Information. Additional data and scripts are available at https://doi.org/10.5061/dryad.0vt4b8h1x and https://doi.org/10.5281/zenodo.6447995.

Associate Editor: C. Burch

Handling Editor: A. G. McAdam