Evolutionary capacity of upper thermal limits: beyond single trait assessments

Thermal tolerance is a key factor contributing to distributional limits in many taxa (Cossins and Bowler, 1987). Temperature is increasingly likely to be a source of strong selective pressure for many organisms, with temperatures expected to increase across the globe over coming decades (IPCC, 2013). Importantly, the close association between environmental and body temperatures means that climate change is likely to impact ectotherms' distribution and abundance (Chown et al., 2010; Colwell et al., 2008; Parmesan and Yohe, 2003), metabolism (Dillon et al., 2010) and therefore risk of extinction (Deutsch et al., 2008; Huey et al., 2009; Sinervo et al., 2010).

Importantly, behavioural thermoregulation may have a limited ability to ameliorate the effects of climate warming in ectotherms (Huey and Pascual, 2009; Huey and Tewksbury, 2009; Kearney et al., 2009; Rego et al., 2010), and temperature is expected to impose significant selection pressures on both ectotherms and endotherms (e.g. Huey et al., 2012). However, whether organisms are able to modify upper thermal limits via evolutionary responses or phenotypic plasticity and mediate their extinction risk remains largely unknown. Inter-specific studies of Drosophila (Kellermann et al., 2012) suggest that some organisms will be unlikely to mediate the effects of global warming by evolving higher thermal limits. However, comprehensive intra-specific assessments of adaptive capacity for upper thermal limits, which provide important insight into contemporary microevolutionary responses to global change, are more limited, and provide mixed support for such conclusions (but see Kelly et al., 2013).

In addition, estimates of upper thermal limits may vary depending on the methodology used (Mitchell and Hoffmann, 2010; Rezende et al., 2011; Santos et al., 2011; Sgrò et al., 2010; Terblanche et al., 2007). Thermal tolerance can be assessed using either static (constant temperature) assays (Hoffmann et al., 2002; Hoffmann et al., 2003) or dynamic (variable temperature) assays that involve gradually heating or cooling an animal from a particular starting temperature until physiological failure, such as knockdown or loss of righting ability (Mitchell and Hoffmann, 2010; Overgaard et al., 2012; Sgrò et al., 2010; Terblanche et al., 2007; Terblanche et al., 2011). Ramping assays are argued to be more ecologically relevant because they are thought to better reflect changes in temperature in the field and because they indicate the activity range for a population under acute conditions experienced in nature. However, the rate of change in temperature used in these assays has been shown to affect predictions about upper thermal limits (Mitchell and Hoffmann, 2010; Sgrò et al., 2010), leading some to question the extent to which different measures adequately or accurately assess upper thermal limits (Rezende et al., 2011; Santos et al., 2011). It should be noted that inferences about upper thermal limits based on static measures have been linked to fitness in the field under hot conditions (Kristensen et al., 2007), suggesting that such measures also capture ecologically important aspects of Drosophila life history in response to extreme thermal stress. Repeated clines in static measures of heat tolerance (Sgrò et al., 2010) and other studies (Loeschcke et al., 2011) also suggest that the ability to mount responses to sudden increases in temperature is ecologically important. Such responses may become increasingly important under climate change, where the frequency and severity of extreme thermal stress is predicted to increase.

Importantly, estimates of adaptive capacity for upper thermal limits can also depend on the methodology used (Chown et al., 2009; Mitchell and Hoffmann, 2010; Rezende et al., 2011). Specifically, Mitchell and Hoffmann (Mitchell and Hoffmann, 2010) found that compared with static measures of heat tolerance, the narrow-sense heritability and evolvability of ramping heat tolerance was significantly reduced in two populations of D. melanogaster. These results imply a constrained evolutionary response to selection imposed by the gradual heating that may also be experienced in natural populations (Hoffmann, 2010; van Heerwaarden et al., 2012). This reduced adaptive capacity was driven by significantly lower levels of additive genetic variance, and not inflated environmental variance (see Chown et al., 2009; Rezende et al., 2011). Nonetheless, criticism of this and other studies (e.g. Sgrò et al., 2010; van Heerwaarden et al., 2012) that have used slower (0.06°C increase per minute) ramping rates to examine upper thermal limits has continued on the grounds that such rates are confounded by increased stochasticity, which will inflate the environmental variance and subsequently reduce narrow-sense heritability estimates, thereby artificially lowering estimates of adaptive capacity for upper thermal limits (Rezende et al., 2011; Santos et al., 2011; Santos et al., 2012).

However, this argument has been framed entirely within a univariate context. This is important, because we know that the ability of a population to respond to selection for increasing heat tolerance will be determined by the patterns of genetic variation and covariation in the traits under selection (Blows and Hoffmann, 2005). If both static and dynamic measures of thermal tolerance significantly covary with one another, then reliance on univariate measures of adaptive capacity may provide inaccurate information about the evolution of upper thermal limits under climate change.

Mitchell and Hoffmann (Mitchell and Hoffmann, 2010) found that the static and ramping measures of heat tolerance were positively correlated across Drosophila species, suggesting perhaps some potential for correlated responses to selection. However, only one study to date has directly assessed the extent to which static versus ramping measures of heat tolerance share a genetic basis. Specifically, using a multivariate approach, we have recently shown (van Heerwaarden and Sgrò, 2013) that D. simulans harbours significant levels of additive genetic variation for both ramping and static measures of heat tolerance, and that both types of traits will contribute to selection for increased heat tolerance in this species via direct and correlated responses. These results contrast with those for D. melanogaster (Mitchell and Hoffmann, 2010), which were based on a univariate assessment of adaptive capacity.

The aim of the present study was therefore to examine the extent to which static and dynamic measures of heat tolerance in D. melanogaster share a genetic basis, and to determine whether the evolution of heat tolerance in natural populations of this widespread species might indeed be constrained by low additive genetic variances within or covariances between different measures of heat tolerance. We did this by estimating additive genetic variances and covariances for three commonly used measures of heat tolerance, static basal and hardened heat knockdown time and ramping (dynamic) heat knockdown time, in a single population of D. melanogaster. Previous work (Sgrò et al., 2010) has shown that all three measures of heat tolerance show parallel clines in D. melanogaster populations from along the east coast of Australia that reflect the action of selection and not neutral processes. Furthermore, we know that temperature is a significant selective agent for Drosophila (Huey and Pascual, 2009; Kellermann et al., 2012; Rego et al., 2010), and the observed clinal patterns in heat tolerance in D. melanogaster (Sgrò et al., 2010) reflect this. However, whether these parallel clinal patterns are the result of independent selection acting on all three measures of heat tolerance or reflect a shared genetic basis remains unknown. We therefore performed a half-sib–full-sib breeding design to empirically assess the additive genetic variance for, and additive genetic covariances between, all three measures of heat tolerance. This allowed us to test whether the predictions arising from previous work (Mitchell and Hoffmann, 2010) in D. melanogaster, where univariate estimates of adaptive capacity suggest a limit for adapting to the gradual increases in temperature that are commonly experienced in nature, hold true when an explicitly multivariate perspective is taken.

We detected significant levels of additive genetic variance only for the two static measures of heat tolerance (Table 1). The additive genetic covariance between basal and hardened knockdown time was positive and significantly different from zero. The covariances between ramping heat knockdown time and both static basal and hardened heat knockdown time were not significantly different from zero (Table 1).

One additive genetic correlation was significantly different from zero: a significant positive genetic correlation was found between basal and hardened heat knockdown time (Table 1). This genetic correlation was significantly different from zero but not one, implying that both traits will show correlated evolutionary responses to selection pressures. None of the genetic correlations between ramping and either static measure of heat tolerance were significantly different from zero (Table 1).

The genetic relationship between all three traits was investigated by examining the additive genetic variance–covariance (G) matrix. An eigen analysis of G was performed to determine how many genetically independent traits (eigenvectors) were represented by the three original traits measured, and how much genetic variance (eigenvalues) was associated with each eigenvector. The eigen analysis revealed that the genetic variance in G was distributed in a single dimension (Table 2). The first eigenvector g_max explained 88.71% of the total additive genetic variance in G. That is, most of the additive genetic variance in all three traits measured can be represented by the first eigenvector, g_max. Basal, hardened and ramping heat knockdown time all loaded positively to g_max. Basal knockdown time made the largest contribution to this vector (Table 2). Ramping knockdown time made a large positive contribution to the second eigenvector, g₂, with a much smaller positive contribution from static hardened heat knockdown time. Static basal heat knockdown made a small negative contribution to this eigenvector. However, this eigenvector only accounted for a small amount of the total additive genetic variance in G (Table 2). The third eigenvector, g₃, only accounted for 4% of the total additive genetic variance in G, with negative contributions from ramping and static basal knockdown time, and a large positive contribution from static hardened heat knockdown time.

Factor analytic modelling tested the statistical significance of the genetic dimensions described by G. Of the three dimensions, only the first was significant (Table 3), and therefore only g_max can be considered to display significant genetic variance.

The first two dimensions were displayed in a biplot to visualise the genetic relationships between the three traits (Fig. 1). In the biplot, the squared length of a vector is the variance explained by the two dimensions, while the cosine angle between vectors is the genetic correlation between them in this two-dimensional space (Smith et al., 2001). Therefore, vectors orientated in the same direction have a high correlation. As the angle between vectors increases, the genetic correlation decreases. Therefore, basal and hardened heat knockdown time have vectors of similar direction and magnitude, reflecting the high pair-wise genetic correlation between them (Table 1). In contrast, ramping heat knockdown time has a vector in a different direction, and the genetic correlation between this trait and both basal and hardened heat knockdown time is reduced (Table 1).

Narrow-sense heritability estimates were only significant for static basal heat knockdown time (Table 4). The non-significant heritability for static hardened heat knockdown time seems to be driven by higher environmental variance, as significant levels of additive genetic variance were detected (Tables 1, 4). In contrast, very low levels of additive genetic variance, rather than inflated environmental variance, seem to be driving the non-significant heritability estimate for ramping heat knockdown time (Table 4). This interpretation is supported by the mean standardised estimates of variation (Table 4), where the coefficients of additive genetic variation and environmental variation were larger for static basal heat knockdown time, and very much smaller for ramping heat knockdown time. The lower evolvability estimate for ramping heat knockdown time suggests that the potential rate of univariate evolutionary change in this trait is less than for either static basal and hardened heat knockdown (Table 4).

Climate change is expected to impose increasing selection on upper thermal limits. There has subsequently been mounting interest in understanding the extent to which the evolution of upper thermal limits might be constrained by low additive genetic variance (Kellermann et al., 2012; Mitchell and Hoffmann, 2010). It has recently been suggested that inferences about an organisms' ability to evolve higher upper thermal limits will depend on how heat tolerance is measured (Chown et al., 2009; Mitchell and Hoffmann, 2010; Rezende et al., 2011; Santos et al., 2012). Specifically, estimates of heritability using ramping assays may be reduced largely because this measure incorporates more stochasticity, thereby increasing the environmental variance (Chown et al., 2009; Rezende et al., 2011; Santos et al., 2012). In the first direct test of these ideas, Mitchell and Hoffmann (Mitchell and Hoffmann, 2010) did indeed show that the narrow-sense heritability for ramping heat knockdown time was not significantly different from zero. However, they clearly showed that this was due to reduced additive genetic variance, and not inflated environmental variance.

However, this debate has occurred within a univariate framework, and ignores the fact that the ability of a population to respond to selection for increasing heat tolerance will be determined by the patterns of genetic variation and covariation in the traits under selection. The results of Mitchell and Hoffmann (Mitchell and Hoffmann, 2010) seem to be at odds with the parallel linear clines in static and ramping female heat tolerance reported for D. melanogaster populations from eastern Australia (Sgrò et al., 2010). Whether the cline in ramping heat knockdown time in D. melanogaster is the result of correlated responses to selection on an unmeasured trait could not be determined, as Mitchell and Hoffmann (Mitchell and Hoffmann, 2010) did not assess the additive genetic covariance between the different measures of heat tolerance.

Thus the motivation of our study was to take a multivariate approach to estimating adaptive capacity for upper thermal limits, and in doing so determine the extent to which adaptive capacity differs across three different, but commonly used, measures of heat tolerance in D. melanogaster. We first showed significant additive genetic variance only for the static measures of heat tolerance. We then demonstrated that the additive genetic variance covariance matrix, G, is of reduced rank, with only the first eigenvector, g_max, explaining most (88%) of the total additive genetic variance in G. Basal and hardened static heat knockdown time made the largest contributions to g_max, with a much smaller contribution by ramping heat knockdown time. This suggests that selection for increased heat tolerance (in the direction of g_max) should result in evolutionary increases in heat tolerance in the two static traits, with a very small response in the ramping measure.

We also showed that static basal and hardened knockdown time were positively genetically correlated to each other. In contrast, the genetic correlations between ramping and both static measures of heat knockdown time were not significantly different from zero. Although Mitchell and Hoffmann (Mitchell and Hoffmann, 2010) suggested that static basal and ramping heat knockdown time might be correlated based on inter-specific correlations between the two traits, our data indicate that adult responses to static and ramping thermal stress are largely genetically independent, which is consistent with recent work in D. simulans (van Heerwaarden and Sgrò, 2013). We also show that ramping heat knockdown time and static hardened heat knockdown time are genetically independent of each other at least in the population of D. melanogaster examined here. These results are also consistent with recent work in D. simulans (van Heerwaarden and Sgrò, 2013). As ramping involves a gradual increase in temperature, it has been suggested that responses to ramping heat stress likely involves hardening responses (Chown et al., 2009; Rezende et al., 2011; Sgrò et al., 2010). However, our results indicate that the two measures are not genetically related.

Our results therefore suggest that the parallel clines observed in all three traits in populations of D. melanogaster from eastern Australia (Sgrò et al., 2010) are likely the result of both direct and correlated responses to selection. This is possible in spite of the non-significant univariate estimates of additive genetic variance for ramping heat knockdown time because this trait still made a positive, although small, contribution to g_max (Lynch and Walsh, 1998), and g_max was the only vector (combination of independent traits) that displayed significant additive genetic variance in multivariate space. It is also possible that the cline in ramping heat knockdown time is the result of historical directional selection that resulted in the depletion of additive genetic variance for this trait over time (Blows and Hoffmann, 2005).

Our univariate analysis was consistent with the results of Mitchell and Hoffmann (Mitchell and Hoffmann, 2010). Narrow-sense heritability estimates for ramping heat knockdown time were not significantly different from zero because of low levels of additive genetic variance for this trait, not because of increased environmental variance (see Santos et al., 2012). Indeed, the environmental variance for ramping heat knockdown time was smaller than for either of the static measures of heat tolerance, as was the mean standardised coefficient of environmental variation. These results are also consistent with those reported for the same measures of heat tolerance in D. simulans (van Heerwaarden and Sgrò, 2013). Taken together, these results suggest that estimates of adaptive capacity for ramping measures of heat tolerance using the slower rates of temperature increase (0.06°C min⁻¹) are not downwardly biased by increased levels of stochasticity as argued by Santos et al. (Santos et al., 2012). Nor are they confounded by starvation or desiccation stress (Overgaard et al., 2012). Why D. melanogaster displays non-significant levels of additive genetic variance for ramping heat knockdown time, in contrast to the similarly widespread D. simulans, is puzzling. Finally, the non-significant narrow-sense heritability for hardened heat knockdown time reported here reflects an increased environmental variance for this trait, because it did display significant additive genetic variance. These results again contrast with those reported for hardened heat knockdown time in D. simulans, where both narrow-sense heritability and the additive genetic variance were significantly different from zero. Whether these results extend to other Drosophila species or insects is not yet known.

It is important to note that a true empirical test of the extent to which the evolution of heat tolerance in D. melanogaster may or may not be constrained by genetic variances or covariances requires an estimate of both the additive genetic variance–covariance matrix (G) and the vector of directional selection gradient, β, for all traits (Lynch and Walsh, 1998). While we have estimated the former, we do not have direct estimates of β for any of the traits examined. We can only infer the role of natural selection from clinal studies of D. melanogaster from eastern Australia that demonstrate parallel clines in all three heat tolerance traits (Sgrò et al., 2010) as well as clines in starvation resistance and body (wing) size (James et al., 1995) that have been shown to result from selection rather than genetic drift. Whether on-going selection for increased heat tolerance in D. melanogaster will result in unconstrained long-term evolutionary responses, or whether the multivariate genetic variance for upper thermal limits might be exhausted, resulting in a selection response plateau (Gilchrist and Huey, 1999), remains to be investigated. Further studies that take an explicitly multivariate perspective to this question are required.

In conclusion, we have shown that D. melanogaster does indeed harbour reduced additive genetic variation for a dynamic measure of heat tolerance. However, our multivariate analysis shows that selection for increased heat tolerance will be possible, largely as the result of direct and correlated responses in static measures of heat tolerances, with only a small contribution from the dynamic measure. Further empirical studies that take an explicitly multivariate perspective to the evolution of thermal tolerance across more taxa are needed to obtain a clearer understanding of adaptive capacity for upper thermal limits.

Drosophila melanogaster Meigen 1830 were collected from Coffs Harbour, a mid-latitude (latitude 30.30°S, longitude 153.12°E) site along the east coast of Australia, using banana baits in March 2012. Fifty field-inseminated females were collected and used to establish iso-female lines in the laboratory at 25°C under a 12 h:12 h light:dark cycle at 95–100% humidity. Species identification was confirmed in the F1 generation by examining males from each iso-female line, to ensure all lines were in fact D. melanogaster. These iso-female lines were then allowed an additional generation in the laboratory to ensure large population sizes in each line prior to setting up a mass-bred population. In the second generation after collection (F2), a mass-bred population was founded with 10 males and 10 females from each of the 50 iso-female lines. This mass-bred population was kept at 25°C under a 12 h:12 h light:dark cycle at 95–100% humidity in three 250 ml bottles containing 60 ml of potato, yeast and sucrose medium.

After five generations of mass rearing (to minimise laboratory adaptation and to ensure the break-up of linkage disequilibrium resulting from crossing iso-female lines), we used a paternal half-sibling breeding design to estimate the additive genetic variances and covariances underlying the different measures of heat tolerance. The parents of the focal flies were reared at controlled densities of 40 eggs per vial, and were collected within 6 h of emergence. Virgin females and males were separated using light (less than 2 min exposure) CO₂ anaesthesia and held in separate vials by sex, at a density of ~20–30 individuals per vial until 4 days old. One-hundred and seventy virgin males (sires) were randomly selected from all holding vials. Each sire was placed in a vial containing 6 ml of food medium and ad libitum live yeast, with five virgin females (dams) and left to mate for 3 days. After this time, each dam was placed individually in a separate vial and allowed to lay eggs for 6 to 8 h to control larval density to no more than 20 larvae per vial, then moved to a fresh vial and allowed to lay eggs for a further 6 h. Densities were ~300–350 flies per bottle to ensure a census population size of 900+ individuals.

Generation six individuals – the focal offspring – were collected within 1 day of emerging in vials and held together in a fresh food vial for a further 48 h to ensure females were mated. After 48 h, females were separated using light CO₂ anaesthesia and allowed to recover for a further 48 h. Two females from each vial were measured for each heat tolerance assay (for a total 100 sires, each mated to five dams, with two offspring per dam measured for each assay of heat tolerance). Flies used in the heat tolerance assays were 5 to 6 days old. The half-sibling breeding experiment was performed at 25°C under a 12 h:12 h light:dark cycle at 95–100% humidity.

Females were placed individually in 5 ml glass vials, and exposed acutely to 38.5°C by immersion in a preheated re-circulating water bath. The hardening treatment involved exposure of flies to 35°C for 30 min followed by recovery at 25°C for 3 h prior to the knockdown assay being performed (van Heerwaarden et al., 2012; van Heerwaarden and Sgrò, 2013). Basal and hardened flies were tested simultaneously. Heat knockdown time was scored as the time taken for all flies to be knocked down and immobilised. Heat knockdown time was assessed over 2 days, with six runs performed each day. The same two people (observers) performed all of the heat knockdown assays across all runs.

Individual females were placed in 5 ml glass vials, which were submerged into a water bath heated to 28°C. The temperature was increased gradually by 0.06°C min⁻¹, which is representative of maximal rates of temperature increase in southeastern Australia (van Heerwaarden et al., 2012) and other parts of the world (Nyamukondiwa and Terblanche, 2010; Terblanche et al., 2011). A data logger (Maxim Integrated i-button DS1923) was submerged into the heated water bath (along with the flies) to record the temperature of the heated water bath throughout the experiment. Resistance was scored as the time taken for all flies to be knocked down, ensuring that all three traits were measured on the same scale. Note that this measure is equivalent to CT_max, because the temperature was increased at the rate of 0.06°C min⁻¹ until all flies had succumbed to heat stress. The ramping assays were performed over 2 days with two runs per day. The same two people (observers) performed all of the heat knockdown assays across all runs.

It has recently been suggested that observer error will affect the estimation of variance components for traits such as heat knockdown time (Castañeda et al., 2012) and that multiple measurements for every individual should be taken, and repeatability statistics reported for thermotolerance. This is not feasible in experiments such as those described here. Instead, we checked for a significant effect of observer on the phenotypic variance of all three measures of heat knockdown time prior to estimating the variance components using Levene's test for equal variances. Observer had a significant effect on the variance of ramping (F_1,932=5.86, P<0.05) and static basal (F_1,964=5.83, P<0.05) heat knockdown time, but not static hardened knockdown time (not shown). To ensure that this did not bias estimates of the variance and covariance components, we variance standardised all of the heat knockdown time data by observer prior to the analyses described below. Run did not affect the variances (not shown), but did affect the means, so was treated as a fixed effect in the model described above.

The additive genetic variance for each trait was first estimated using a univariate model. Log likelihood ratio tests were performed, where the final model for each trait was compared with a model specifying σ²_S to be zero, to determine whether levels of additive genetic variance for each trait were significantly different from zero (Littell et al., 1996; McGuigan et al., 2011; Simonsen and Stinchcombe, 2010). We then estimated the unconstrained G matrix. In both cases, the variance at both the sire, d_s, and the dam, d_d, levels was modelled using an unstructured covariance matrix. The additive genetic variance and covariance components of G were individually tested for significance from zero by performing log likelihood ratio tests where the final models for each trait were compared with models specifying σ²_S and the sire-level covariances (COV_S) to be zero (Littell et al., 1996; McGuigan et al., 2011; Simonsen and Stinchcombe, 2010).

To examine the distribution of genetic variance in multivariate space, two complimentary approaches were utilised to estimate the number of dimensions of G.

To determine how many genetically independent traits (eigenvectors) were represented by the original traits (phenotypes) actually measured, and how much genetic variance (eigenvalues) was associated with each independent set of eigenvectors (traits), eigen analysis of the unconstrained additive genetic variance covariance matrix, G, was performed using the matrix analysis option implemented in the Microsoft Excel add-in PopTools (Hood, 2010). The eigenvector with the largest eigenvalue (g_max) (Schluter, 1996) is the vector explaining most of the additive genetic variance in the G matrix. The traits (eigenvectors) that load onto g_max represent independent trait combinations that together account for the most additive genetic variance in G.

Factor analytic modelling was then used to determine how many factors (dimensions) of G are required to describe all of the genetic variation across the three traits examined (Hine and Blows, 2006; McGuigan and Blows, 2010). Factor analytic modelling reports the fit of the model for the number of factors (dimensions) specified using a factor analytic covariance structure rather than the unstructured covariance at the sire level. In order to determine the number of significant genetic dimensions of G, the fit of the model was determined for one through to three genetic dimensions. The significance of the nth dimension was then tested by comparing the model with a model with n−1 dimensions in a hierarchical fashion (McGuigan and Blows, 2010). For example, to test for the significance of the second dimension of G, the model with two dimensions was compared with the model with one dimension. The χ² value is the difference in −2 log-likelihood scores between the two models, and the degrees of freedom are the difference in the number of covariance parameters associated with the two models. If the test produced a significant χ², then adding the second dimension produced a significantly improved model. If the test produced a non-significant χ², then adding the second dimension produced a significantly inferior model.

To complement the multivariate methods described above, we estimated the additive genetic correlation between all three traits using the MIXED procedure of SAS (SAS Institute). Log likelihood ratio tests were used to test whether any of the additive genetic correlations were significantly different from both zero and one (Littell et al., 1996; Simonsen and Stinchcombe, 2010).

To directly assess the extent to which univariate measures of evolvability for heat tolerance reflected the multivariate analyses described above, we also estimated the narrow-sense heritability for each trait. Narrow-sense heritability for each trait was estimated as the additive genetic variance (V_A) divided by the total phenotypic variance (V_P) (Falconer and Mackay, 1996; Lynch and Walsh, 1998). The evolvability, I_A, of each trait was estimated as the additive genetic variance divided by the square of the trait mean following Hansen et al. (Hansen et al., 2011), as its numerical estimate can be interpreted as the percent change in a trait under a unit strength of selection.

We would like to thank Fiona Cockerell, Richard Foo Heng Lee and Allannah Clemson for technical support, and Ary Hoffmann for comments on an earlier version.