Macroevolution of thermal tolerance in intertidal crabs from Neotropical provinces: A phylogenetic comparative evaluation of critical limits

Abstract Thermal tolerance underpins most biogeographical patterns in ectothermic animals. Macroevolutionary patterns of thermal limits have been historically evaluated, but a role for the phylogenetic component in physiological variation has been neglected. Three marine zoogeographical provinces are recognized throughout the Neotropical region based on mean seawater temperature (T m): the Brazilian (T m = 26 °C), Argentinian (T m = 15 °C), and Magellanic (T m = 9 °C) provinces. Microhabitat temperature (MHT) was measured, and the upper (UL 50) and lower (LL 50) critical thermal limits were established for 12 eubrachyuran crab species from intertidal zones within these three provinces. A molecular phylogenetic analysis was performed by maximum likelihood using the 16S mitochondrial gene, also considering other representative species to enable comparative evaluations. We tested for: (1) phylogenetic pattern of MHT, UL 50, and LL 50; (2) effect of zoogeographical province on the evolution of both limits; and (3) evolutionary correlation between MHT and thermal limits. MHT and UL 50 showed strong phylogenetic signal at the species level while LL 50 was unrelated to phylogeny, suggesting a more plastic evolution. Province seems to have affected the evolution of thermal tolerance, and only UL 50 was dependent on MHT. UL 50 was similar between the two northern provinces compared to the southernmost while LL 50 differed markedly among provinces. Apparently, critical limits are subject to different environmental pressures and thus manifest unique evolutionary histories. An asymmetrical macroevolutionary scenario for eubrachyuran thermal tolerance seems likely, as the critical thermal limits are differentially inherited and environmentally driven.


| INTRODUCTION
Thermal tolerance in animals depends on systemic through subcellular functions of resistance, entailing complex physiological interactions that underpin most biogeographical distribution patterns of life on Earth. Such tolerance results from physiological mechanisms usually manifested under extreme environmental conditions, which can be considered as descriptors of spatial discontinuities among ectotherms. Critical temperatures are those at which the loss of essential motor capability occurs, resulting in disorganized locomotor activity, affecting an organism's ecological function (Cowles & Bogert, 1944;Lutterschmidt & Hutchison, 1997;Sunday, Bates, & Dulvy, 2010). Survival under such extreme conditions is constrained by the availability of energy generated anaerobically (Zielinski and Pörtner, 1996;Sommer, Klein, & Pörtner, 1997), in addition to the mobilization of heat-shock proteins and antioxidant defenses (Heise, Puntarulo, Nikinmaa, Abele, & Pörtner, 2006;Tomanek, 2005).
Upper critical thermal limits (UL 50 ) can thus predict latitude, vertical position in an intertidal environment, or the biotope occupied by different species, as they correlate positively with microhabitat temperature and are associated with mechanisms of thermal tolerance. In contrast, lower critical thermal limits (LL 50 ) are poorly explored in a latitudinal context, although there is a link between lower thermal resistance and higher latitudes (Demeusey, 1957;Tashian, 1956;Vernberg & Tashian, 1959;Vernberg & Vernberg, 1967). Lower limits seem to correlate better with ecological factors such as predation, competition, and substratum preference than with physiological limits (Jensen & Armstrong, 1991;Paine, 1974).
Most studies have characterized thermal tolerance in a speciesspecific manner, or have compared tolerance over a broad range of taxonomic levels, neglecting the influence of the phylogenetic component on physiological variation. Interspecific comparisons must consider the phylogenetic relationships among species due to their lack of statistical independence owing to shared ancestry (Felsenstein, 1985;Garland, Bennett, & Rezende, 2005;Garland & Ives, 2000).
Different lineages tend to evolve independently of one another, and physiological variation among species thus increases as a function of phylogenetic distance, rendering closely related species more similar for reasons of ancestry and not necessarily as a consequence of environmental pressures (Harvey & Pagel, 1991;Rezende & Diniz-Filho, 2012). While it is essential to retrieve the evolution of thermal tolerance, the use of phylogenetic information when evaluating critical thermal limits across a large geographical distribution is rarely encountered in crustacean studies (see Stillman, 2002;Stillman & Somero, 2000).
The American continent is classified into 16 marine zoogeographical provinces, which are defined as part of the neritic zone characterized by a narrow range of water temperatures, and containing a fairly constant decapod crustacean fauna (see Boschi, 2000a,b).
Specifically, the eastern coast of South America is divided into three provinces: (1) the Brazilian zoogeographical province, delimited by the mouth of Orinoco River, Venezuela (9° N), and extending to Cabo Frio/RJ, Brazil (23° S) (Briggs, 1974); (2) (Carcelles & Williamson, 1951). The Brazilian (22-30 °C) and Magellanic (4-15 °C) provinces are more stenothermic than the Argentinian province (8-23 °C), because the first is dominated by the South Equatorial Current (Thurman, Faria, & McNamara, 2013), and the second by the homogeneous mass of subantarctic waters (Boschi, 2000a). The Argentinian province is the most eurythermic region, characterized by a mixture of cold water from the Malvinas Current and warmer waters from the Brazil Current, and tends to include species more tolerant of temperature variation (Boschi, 2000a).
Here, we propose an evolutionary history of thermal tolerance in 12 intertidal, eubrachyuran crab species, selected based on their ample distribution across the three zoogeographical provinces along the eastern coast of South America. A phylogenetic analysis was performed using partial sequences of the 16S mitochondrial genes from the selected species, including other Brachyura and Anomura for comparative evaluation. We tested for: (1) phylogenetic patterns of microhabitat temperature (MHT), UL 50 , and LL 50 ; (2) an effect of zoogeographical province on the evolution of both critical thermal limits; and (3) an evolutionary correlation between MHT and these limits.
The macroevolutionary pattern of thermal tolerance is discussed in a biogeographical context, particularly regarding the phylogenetic and environmental components.

| Crab species and laboratory maintenance
The crab species chosen were collected from three distinct thermal provinces along the eastern coast of South America: (1)  Approximately 80 adult, intermolt crabs of either sex from each of the 12 species were collected manually from mangroves, salt marshes, and sandy or rocky beaches along the eastern coast of South America ( Figure 1). Substrate temperature was measured using an infrared, digital thermometer (Icel TD-965, 0.1 °C precision) during crab collections (7 ≤ N ≤ 15 measurements per species). Specimens were obtained during the morning low tides at the end (February-March) of the summers of 2013 and 2014. In the laboratory, the crabs were held in plastic boxes containing a 4-to 8-mm-deep film of seawater in incubators or a water bath (Fanem 347 DCG; Solab SL 224; Polystat 12002, F I G U R E 1 Collecting sites within the three zoogeographical provinces on the eastern coast of South American (sensu Boschi, 2000a) (left) and a phylogenetic hypothesis for selected brachyuran crab species (right). The Brazilian zoogeographical province is delimited by the mouth of Orinoco River, Venezuela (9° N), and Cabo Frio/RJ, Brazil (23° S); the Argentinian province lies between Cabo Frio and Rawson, Argentina (43° S); and the Magellanic province extends from Rawson to Ushuaia, Argentina (55° S). Aratus pisonii, Cardisoma guanhumi, Goniopsis cruentata, Ocypode quadrata, Pachygrapsus transversus, Uca maracoani, and Ucides cordatus were collected from the Brazilian province (red, ≈7.8° S/34.8° W, Ilha de Itamaracá or Itapissuma, PE, Brazil); Armases rubripes, Neohelice granulata, and Uca uruguayensis from the Argentinian province (cyan, ≈32.1° S/52.1° W, Rio Grande, RS, Brazil); and Acanthocyclus albatrossis and Halicarcinus planatus from the Magellanic province (blue, ≈53.2° S/ 67.2° W, Ushuaia, Tierra del Fuego, Argentina). The molecular phylogeny was generated using a maximum likelihood search method employing a partial sequence of the mitochondrial 16S gene (592 base pairs, GTR + Γ + I model) for 36 brachyuran and anomuran (outgroup) species. The final alignment of the fragments consisted of 15 novel sequences and 21 sequences obtained from NCBI GenBank. The nodal values represent clade support (bootstrap) Palmer Instrument Company) at the annual mean seawater temperature of the respective provinces (sensu Boschi, 2000a,b, Boschi & Gavio, 2005 for 3 days prior to experiments, i. e. Brazilian province, 26 °C, 12-hr light/12-hr dark photoperiod; Argentinian province, 15 °C, 14-hr light/10-hr dark photoperiod; and Magellanic province, 9 °C, 14-hr light/10-hr dark photoperiod. were chosen, such that mortality was greater than 50% in at least one temperature and less than 50% in another.This 6-hr exposure period represents the mean emergence time of a crab at low tide and thus exposure to the greatest thermal amplitude occurring during the tidal cycle. Crabs were examined every 30 min and were considered "dead" if they could no longer right themselves after 30 sec. The critical thermal limits define the temperatures at which 50% of the crabs lose their key motor functions, being unable to right themselves when placed in a supine position, revealing disorganized locomotor activity. In this case, organismic functions are affected as "the animal loses its ability to escape from conditions that will promptly lead to its death" (Cowles & Bogert, 1944). Critical thermal limits were determined by probit analysis (Finney, 1971), adjusting the percentage survival in a linear regression model.

| Phylogenetic analysis
A partial fragment of the mitochondrial 16S ribosomal gene (≈590 base pairs) was used as a molecular marker. This gene has been widely used in studies of phylogenetic inferences in decapod crustaceans (e.g., Mantelatto, Pardo, Pileggi, & Felder, 2009;Mantelatto, Robles, Schubart, & Felder, 2009;Pileggi & Mantelatto, 2010;Schubart, Cuesta, Diesel, & Felder, 2000). Twelve new sequences were obtained, one for each species. All sequences obtained were submitted to NCBI GenBank (Table 1). Other species of Brachyura were included in the alignment to provide a more consistent analysis (Table 1) following previous phylogenetic methodologies (Tsang et al., 2014).
Some species of Anomura were used as outgroups (Table 1).
Phylogenetic reconstructions were performed using the maximum likelihood search method (Felsenstein, 1973(Felsenstein, , 1981 in RAxML 7.2.7 (Randomized A(x)ccelerated Maximum Likelihood) (Stamatakis, 2006) implemented on the CIPRES (Cyberinfrastructure for Phylogenetic Research) system (http://www.phylo.org). The substitution model used was GTR + ∫ + I [General Time Reversible (Tavaré, 1986) + Gama + invariables sites] as specified by RAxML. Topology consistency was measured using a rapid bootstrap method (1,000 replicates) (Stamatakis, Hoover, & Rougemont, 2008) and only confidence values greater than 50% were reported. The topologies were visualized and edited using MEGA 5.1 software. The phylogeny obtained is provided in Figure 1 and was pruned to match the number of species (12) for which thermal tolerances were characterized. Topology and branch lengths among species were maintained, as is standard practice for comparative evaluations employing nonmissing trait data (Kembel et al., 2010).

| Comparative analyses
Comparisons among closely related species, when performed using an appropriate phylogeny, real measure of branch lengths and good fitting model of evolution, reduce type I error, decreasing the probability of detecting spurious correlations owing to uncontrolled factors.
This also reduces type II error, adding stronger statistical power allowing the detection of correlated evolution between ecological/physi- Hypotheses regarding the effect of zoogeographical province and microhabitat temperature on the evolution of upper and lower critical thermal limits were tested using a phylogenetic generalized least squares (PGLS) model. Using categorical and quantitative predictors (PGLS ANCOVA), phylogenetically correlated residual variation among species (Garland & Ives, 2000;Grafen, 1989;Lavin, Karasov, Ives, Middleton, & Garland, 2008) is assumed, which is ideal for comparative studies concerning mean values for species. To avoid interference from uncontrolled environmental conditions (Wildt & Olli 1978), MHT was treated as a covariable, which removes its effects, corrects mean values, and improves the power of hypothesis testing regarding the effect of the discrete predictor (province). Simultaneously, we estimated selection strength (α) of the Ornstein-Uhlenbeck (O-U) model of evolution, in which the minimum and maximum values of a trait are limited, stabilizing them toward a central point of variation (Butler & King, 2004;Diniz-Filho, 2001;Nunn, 2011;Revell, 2010). When α is null, covariance among species is a consequence of stochastic evolutionary changes, or results from randomly shifting optima driven by natural selection, as modeled by Brownian motion (Butler & King, 2004;Hansen, Pienaar, & Orzack, 2008). The Holm-Bonferroni post hoc procedure was used to compare multiple means, P-values being corrected by phylogenetic simulations (Revell, 2012). Analyses were performed using the R environment (R Core Team 2015), and its nlme (Pinheiro, Bates, DebRoy, Sarkar & R Core Team, 2016) and ape (Paradis, Claude, & Strimmer, 2004) packages, with the minimum significance level set at p = .05.

| RESULTS AND DISCUSSION
The microhabitat temperatures (  T A B L E 2 Habitat characteristics, microhabitat temperatures (7 ≤ N ≤15), and lower (LL 50 ) and upper (UL 50 ) critical thermal limits for 12 eubrachyuran crab species collected from three zoogeographical provinces along the eastern coast of South America the evolution of the thermal niche of the eubrachyuran crabs examined here.
Both upper and lower critical thermal limits appear to display different phylogenetic patterns (Figure 2). LL 50 seemed to show neither similarities nor dissimilarities between pairs of closely or distantly related species, although there was a tendency toward negative autocorrelation. UL 50 appeared to be strongly structured at the more apical and deeper hierarchical levels: Closely related species tended to share similar upper limits, while more distant species exhibited significant differences ( Figure 2). This pattern of decreasing phylogenetic cor-  Figure 3).
The evolution of LL 50 is more plastic than that of UL 50 , as sug- Fiddler crabs (Uca) from tropical regions exhibit higher upper thermal limits than representatives from subtropical and temperate zones (Vernberg & Tashian, 1959;Vernberg & Vernberg, 1967). Differences F I G U R E 2 Phylogenetic correlogram over three distance classes, using Moran's I coefficients (± SEM), for microhabitat temperatures (MHT) and lower (LL 50 ) and upper (UL 50 ) critical thermal limits of 12 eubrachyuran crab species from three different zoogeographical provinces (see Figure 1). Horizontal line at I = 0.0 indicates no correlation with phylogeny. LL 50 does not correlate with phylogeny at any hierarchical level, revealing labile evolution; MHT and UL 50 exhibit strong phylogenetic signal at the species level, which decreases with phylogenetic distance, suggesting that closely related species share similar thermal niches and mechanisms of hightemperature tolerance F I G U R E 3 Lower (LL 50 ) and upper (UL 50 ) critical thermal limits for 12 eubrachyuran crab species from three zoogeographical provinces on the eastern coast of South America (see Figure 1).  (Cuculescu et al., 1995). The most complete macroevolutionary scenario for upper thermal limits proposed to date refers to anomurans (Petrolisthes): In species occupying the upper intertidal, and in those from tropical zones, the evolution of UL 50 is linked to surface water temperatures and maximal microhabitat temperatures (Stillman, 2002;Stillman & Somero, 2000). With regard to LL 50 , populations of fiddler crabs U. pugilator, U. pugnax, and U. uruguayensis from cooler microhabitats survive longer at lower temperatures than do populations from warmer regions, and than other congeneric species of exclusively tropical distribution (Demeusey, 1957;Tashian, 1956;Vernberg & Tashian, 1959;Vernberg & Vernberg, 1967).
The evolution of thermal tolerance in invertebrates depends on the proportion and types of fatty acids in their cell membranes (Hazel & Williams, 1990;Kates, Moldoveanu, & Stewart, 1993), and on a positive balance between the capability for oxygen supply and oxygen demand (Pörtner, 2001(Pörtner, , 2002. A higher percentage of polyunsaturated fatty acids in winter than in summer appears to underpin the physical properties of the membranes affecting regulation of intracellular osmolality, as transport enzyme kinetics are dependent on membrane fluidity and on intracellular ionic composition (Hazel & Williams, 1990;Murata & Wada, 1995). This phenomenon is known as the "homeoviscous adaptation" (Hazel & Williams, 1990;Kates et al., 1993). Further, maintenance of aerobic scope is linked to an efficient cardiorespiratory system in supplying sufficient oxygen for mitochondrial demands, in addition to the heat-shock response and antioxidant capability in providing protection against protein denaturation and reactive oxygen species (Peck, Webb, & Bailey, 2004;Pörtner, 2001;Zielinski & Pörtner, 1996). Thus, strong selection pressure at the cellular and subcellular levels may explain the lower LL 50 observed in the crab species from Magellanic province, as well as the more elevated UL 50 seen in the eubrachyurans from the two northernmost provinces. These findings reveal the homoplastic nature of the physiological mechanisms of thermal tolerance arrayed against acute thermal challenge.
Thus, the limits of thermal tolerance in eubrachyuran crabs appear to manifest distinct physiological constraints, are subject to different environmental pressures, and show unique evolutionary histories. The upper and lower critical thermal limits are informative descriptors of the discontinuous distribution of crabs throughout the intertidal zone of eastern South America, as they reflect different functional properties of an integrated system. The macroevolutionary landscape explored here suggests an asymmetrical scenario for eubrachyuran thermal tolerance, because the critical thermal limits are differentially inherited and environmentally driven.