Mitochondrial genotype influences the response to cold stress in the European green crab, Carcinus maenas

Hybrid zones are geographic regions where distinct lineages or species meet and mate, often over many generations (Harrison, 1993). The genetic novelty maintained in hybrid zones provides unique opportunities to examine the linkages between genotype and phenotype (Rieseberg et al., 1999). In particular, hybrid zones provide natural experiments to examine how interactions between novel combinations of genes affect physiological and ecological performance. These include interactions between nuclear genes that result from hybridization and backcrossing (Pereira et al., 2014), and interactions across genomes as a result of new cyto-nuclear combinations (Dowling et al., 2008).

Before the Anthropocene, most hybrid zones resulted from shifting range limits driven by Pleiostocene climate change that brought formerly non-overlapping lineages or recently evolved species into secondary contact (Hewitt, 2000). More recently, hybrid zones may also be established by repeated biological invasions into the same geographic region from different source populations. Such is the case of the European green crab, Carcinus maenas (Linnaeus 1758), which has been introduced into all the major continents of Earth, with the exception of Antarctica, over last two centuries (Darling et al., 2008). On the East Coast of North America, a double invasion has established two distinct genetic lineages across a north–south transect (Roman, 2006). The species was first sighted in coastal Massachusetts in 1817, after which its range slowly expanded north until ∼1960 when it appeared to reach its northern range limit near Halifax, Nova Scotia (Carlton and Cohen, 2003). In the early 1980s, this range limit began to move northward, when C. maenas was sighted further north and in much colder waters (Roman, 2006). Since 1980, C. maenas has expanded northward by roughly 725 km, now spanning the coastal waters of Nova Scotia, Prince Edward Island and a significant portion of Newfoundland (Blakeslee et al., 2010). Analysis of nuclear and mitochondrial DNA (mtDNA) has revealed that this rapid northern expansion was most likely facilitated by a second genetically distinct introduction (Roman, 2006). The first introduction has been traced back to a likely source population in southern or central Europe, while the second introduction consisted of individuals from a northern European source population. Over the last 10 years, the northern lineage has been spreading southward, while the southern lineage has been spreading northward, leading to a region of hybridization and admixture located along the eastern Gulf of Maine and eastern Nova Scotia (Pringle et al., 2011; Darling et al., 2014). The structure of this zone presents a number of interesting experimental opportunities to examine links between genotype and phenotype. For example, mitochondrial haplotypes are spreading more rapidly than nuclear markers, leading to strong mito-nuclear discordance outside the ‘nuclear’ hybrid zone, and novel mito-nuclear combinations within the region of nuclear admixture (Darling et al., 2014; Jeffery et al., 2017). Furthermore, the two lineages are mixing across an ecological seascape that varies in a number of physical and biological factors. Perhaps most importantly, seasonal temperature dynamics vary considerably from Cape Cod to the Gulf of St Lawrence – including a strong gradient in average winter temperature that reaches minima of 0–2°C in the Gulf of St Lawrence (Lehnert et al., 2018).

Previous analyses of hybridization dynamics in this complex of closely related populations have largely assumed that mitochondrial variation within and among populations is selectively neutral, and that a number of stochastic demographic factors can explain different rates of migration of genes located in the mitochondria versus the nucleus (Darling et al., 2014). However, to our knowledge there have been no direct studies of how genotypic variation (either mitochondrial or nuclear) might translate into phenotypic variation in this system. There is some indirect evidence that genetic variation may be linked to thermal stress phenotypes from population-level studies. Tepolt and Somero (2014) showed that C. maenas collected from different regions within Europe and the USA had distinct responses to heat and cold stress. These same populations showed fairly strong population differentiation in transcriptome-derived single nucleotide polymorphisms (SNPs) across thermal gradients, though sample sizes were too small to test individual-level genotype–phenotype associations conclusively (Tepolt, 2014; Tepolt and Palumbi, 2015). Other studies have found differences in foraging and competitive behavior between C. maenas populations from sites with different histories and presumably different genetic backgrounds, leading to the hypothesis that green crabs with ‘northern’ genetic backgrounds may outcompete their southern brethren (Haarr and Rochette, 2012; Rossong et al., 2012). Taken together, this work raises the possibility that thermal stress thresholds in C. maenas have genetic components.

In this study, we tested the hypothesis that mitochondrial, nuclear and mito-nuclear interactions can explain differential responses to cold temperature stress, using natural variation in the C. maenas hybrid zone in the Gulf of Maine. We chose to focus on mitochondrial variation because the mitochondrial genome codes for a number of key genes that control physiological processes. Variation in mtDNA has been shown to have effects on fitness, often linked to thermal tolerance, in a variety of different animal species. These include flies (Camus et al., 2017), tuna (Dalziel et al., 2006), hares (Ben Slimen et al., 2017) and humans (Balloux et al., 2009). In mammals, variance in the ATP6 mitochondrial gene has been connected to shifts in temperature tolerance. Mechanistically, ATP6 variants reduce the coupling efficiency of oxidative phosphorylation, which increases heat production while simultaneously increasing basal metabolic rate (Ballard and Whitlock, 2004). This mechanism results in a trade-off – an individual with an ATP6 variant has increased cold tolerance at the expense of a higher metabolic rate in more temperate environments – and illustrates how mtDNA can be under selection with thermal tolerance as the selective factor. Furthermore, the mitochondrial genome has a great deal of interaction with the nuclear genome. Mitochondrial proteins, which play an essential role in the oxidative phosphorylation (OXPHOS) system, also interact extensively with nuclear proteins (Burton and Barreto, 2012). The mitochondrial genome is also particularly prone to accumulating deleterious point mutations, and the nuclear genome can take over the function of mitochondrial genes that are knocked out by deleterious substitutions, which may explain the long-term shrinkage of the mitochondrial genome (Rand et al., 2004; Gray, 2011). This tight co-evolution between mitochondrial and nuclear genes can break down in hybrid zones when mtDNA and nuclear DNA are no longer performing complementary roles (Burton and Barreto, 2012). Thus, we expect mito-nuclear interactions to be important in a variety of traits linked to fitness in a broad variety of systems, and that these interactions are most likely to express themselves in hybrid zones, where new combinations of mitochondrial genomes and nuclear genes are generated.

Our nuclear markers target regions identified in ongoing work on temperature-associated differentiation in cardiac transcriptomes between northern and southern C. maenas populations in Europe and the East Coast (Tepolt and Palumbi, 2015). We designed primers for several ‘outlier loci’ identified in that study [loci with high F_ST (inbreeding coefficient of the subpopulation compared with the total population) values across these regions] and used PCR primers to sequence and genotype two SNP markers in nuclear regions that showed strong frequency differences between northern and southern Europe. Variation at one of these SNPs is characteristic of southern Europe, while variation at the other SNP is characteristic of northern Europe. Thus, by capitalizing on recombination that is naturally occurring within the C. maenas nuclear genome, and between nuclear and mitochondrial genes, we tested whether mitochondrial and selected nuclear genotypes were associated with cold stress phenotypes in wild-caught crabs collected along the hybrid zone on the east coast of North America.

To make sure that our experimental assay included genetic variation that was representative of the two ends and middle of the genetic cline, we conducted preliminary sampling and genotyping in Harpswell, ME, USA (N 43°47′27″, W 69°57′34″) and Kent Island, NB, Canada (N 44°33′00″, W 66°45′23″) and found that Harpswell had a high frequency of southern mt-CO1 haplotypes while Kent Island was near the cline center with mixtures of southern and northern haplotypes. To make sure we sampled from the northern end of the cline, we chose Pomquet, NS, Canada (N 45°38′27″, W 61°48′41″) based on previous studies that revealed a high frequency of northern mt-CO1 haplotypes at this site (Roman, 2006). In addition to the 251 crabs collected for the experimental assay, we sequenced 15 crabs from Halifax, NS, Canada (N 44°38′28″, W 63°55′45″) and 22 crabs from Isles of Shoals, ME, USA (N 42°59′23″, W 70°36′58″). These additional samples were used solely to establish geographical clines in our molecular markers and were not used in the experimental assay (Fig. 1). Crabs were collected either by hand in the low intertidal zone or via baited traps placed in the shallow subtidal zone. In Canada, crabs were collected under Fisheries and Oceans Canada Fishing License no. 342021 and 345689. In Maine, crabs were collected under Department of Marine Resources Special License no. ME 2018-24-00.

The righting response – the time taken for a crab to right itself after being placed on its carapace, ventral-side up – has been shown to be an excellent proxy for other physiological stress variables in a variety of decapods, including C. maenas (Cuculescu et al., 1998; Ern et al., 2015). It is also quick to assay, making it easy to test large numbers of crabs. For these reasons, righting response was chosen as the measure of cold tolerance. The cold tolerance assay was conducted at the Schiller Coastal Studies Center in Harpswell, ME, USA. To control for acclimatization, which can have a strong impact on cold tolerance (Tepolt and Somero, 2014), crabs were acclimated in a lab-based flowing seawater system for 3–5 weeks before experimentation at 15.5±0.1°C with at a flowrate of 3.79 l min⁻¹. Each tank held crabs from one site at a density of 28–60 crabs per 379 l tank. Tanks were exposed to a 12 h:12 h light:dark cycle with overhead fluorescent light fixtures, and crabs were fed with herring or shrimp three times a week. To ensure that differential food consumption did not impact behavior, crabs were starved 4–5 days before the start of the testing.

To assay the righting response, crabs were placed in a large plastic tote filled with aerated seawater cooled to 4.5°C in a Percival incubator (Percival Scientific, Perry, IA, USA). For each trial, the individual crabs were subjected to a 5 min cooling period, after which a metal ruler was used to flip the crab onto its back, and the time to right was measured with a stopwatch. All cold tolerance assays were conducted between 07:00 h and 13:00 h on a single day. The time to right was defined as the period between when the animal was flipped on its dorsal side and when the rearmost dactyls maintained contact with the surface of the tote. As the majority of crabs either flipped rapidly within the 2 min trial or failed altogether to right themselves, we used a binary variable for righting response over the continuous variable righting time. After righting response was scored, a variety of morphological parameters – sex, carapace width, number of legs and morph type (red or green) – were recorded for each crab and used as co-variables in statistical models. A single leg of each crab was preserved in 95% ethanol for DNA analysis.

We sequenced a portion of the mitochondrial CO1 gene (mt-CO1) as a marker for mitochondrial haplotype using the primers and PCR conditions published in Roman and Palumbi (2004). The two nuclear markers used targeted coding regions and flanking regulatory regions isolated and sequenced by Tepolt and Palumbi (2015). Three nuclear markers were initially developed from contigs that were identified as outliers in a gene scan between populations located across latitudinal clines in Europe and North America (Dryad repository: https://doi.org/10.5061/dryad.g8b96/1). These three initial outlier SNPs were selected as likely deriving from the second C. maenas introduction to the east coast. In samples collected in 2011–2012, these SNPs were variable in the Canadian province of Newfoundland and in Norway but fixed for the major allele in Portugal and in the American states of New Jersey and central Maine (Tepolt and Palumbi, 2015). We designed PCR primers to target SNPs within these contigs, Sanger sequenced a set of samples from collection sites that spanned the Gulf of Maine, and found that all three SNP markers had clinal structures with cline centers located around Penobscot Bay in the Gulf of Maine (Voss, 2016). We could reliably annotate one of these three loci as ubiquitin-conjugating enzyme E2 H or UBE2H, and chose it for additional sequencing and tests for association with the cold stress phenotype. The second nuclear marker we used for this study was also identified in transcriptome-based outlier analyses, and was selected as showing a clinal pattern across Europe and North America. In contrast to UBE2H, this SNP appears to have been present in the initial introduction to the east coast; in 2011–2012, it was variable in Portugal, New Jersey and Maine but fixed in Newfoundland and Norway (Tepolt and Palumbi, 2015). This SNP was located in a gene annotated as structural maintenance of chromosomes protein 3, or SMC. Primer sequences and PCR conditions are listed in Table S1.

DNA extraction was performed with a Qiagen DNEasy Blood & Tissue Kit (Valencia, CA, USA) following the manufacturer's instructions. Amplified products were prepared for sequencing using an exonuclease and shrimp alkaline phosphatase incubation (Exo-SAP) and Sanger sequenced using the GeneWiz sequencing laboratory (South Plainfield, NJ, USA). The resulting chromatograms were edited and aligned using the software package Geneious (Biomatters Ltd, Auckland, New Zealand). The nuclear sequences were scored for SNPs by inspecting the chromatograms for a single versus double peak, the latter indicative of heterozygous status. For the mt-CO1 gene, five unique haplotypes were identified from all samples. To identify the evolutionary relationships between our sample sequences and those from Europe, reference mt-CO1 sequences of C. maenas within their native range were downloaded from GenBank (accession nos: AY616437.1–AY616445.1, JQ305941.1–JQ305943.1, JQ306002.1–JQ306003.1, KF369118.1) and combined with our sequence data to construct a TCS haplotype network using the software PopART (http://popart.otago.ac.nz/).

For the two nuclear markers, we tested for Hardy–Weinberg expectations within loci, and linkage between loci, using Genepop (http://genepop.curtin.edu.au/), which implements the estimators of Raymond and Rousset (1995).

To analyze the effects of both categorical (sex, mitochondrial haplotype, nuclear genotype, morph color, number of missing limbs) and continuous (carapace width) independent variables on the righting response, we used a multiple logistic regression model implemented by the R-package PerformanceAnalytics (https://CRAN.R-project.org/package=PerformanceAnalytics). As explained earlier, the righting response was categorical, with crabs either righting or failing to right within the 2 min trial period. Two models were run. The first included all reproductive, ecological and morphological variables (sex, morph color type, number of legs, sampling location and carapace width) to determine whether any of these variables were associated with the righting response. The results of this model were used to select a smaller set of relevant variables for the final model, which included sex, sampling location and three genetic terms for mitochondrial and nuclear variation. Mitochondrial variation was binned into three haplogroups (A, B and C; see Fig. 2). Nuclear variation was binned into three categories corresponding to the three genotypic states for each biallelic SNP locus.

Five unique mt-CO1 haplotypes were sampled from the three Canadian and two USA sites and binned into three haplogroups based on evolutionary similarity (Fig. 2). The two sites in Nova Scotia had a high frequency of two haplotypes closely related or identical to haplotypes commonly found in northern Europe (B1, B2). These two haplotypes were binned within haplogroup B as ‘northern’ haplotypes. The Nova Scotia samples also had a high frequency of a third haplotype related to a haplotype previously observed in Portugal, and believed to have come to the east coast with the second C. maenas introduction (Roman, 2006). This is a pan-European haplotype found to comprise roughly 30–40% of Nova Scotia haplotypes in a previous study (Darling et al., 2008), which agrees with the frequency observed in this study. However, as it was not closely related to either of the northern haplotypes (B1 and B2) it was assigned to haplogroup C. The three haplotypes of the B and C haplogroups dominated the two Nova Scotia sites. The central (Kent Island) and two more southern sites (Harpswell and Isles of Shoals) were dominated by the A1 haplotype, identical to a GenBank sample from France (accession no. AY616437), and a defining haplotype of the initial C. maenas introduction to the east coast (Roman, 2006). The remaining sampled haplotype (A2) differed by one mutation from A1, and therefore was included in haplogroup A. Translation of these five haplotypes revealed that all mutations were synonymous, yielding identical protein sequences from these partial mt-CO1 DNA sequences.

Within each of the five sampling sites, we could not reject Hardy–Weinberg expectations within loci for either nuclear marker, and there was no statistical evidence for linkage between UBE2H and SMC (Table S2). We did observe slightly positive F_IS (inbreeding coefficient of an individual relative to the subpopulation) values for UBE2H among all sites, suggesting mild heterozygote deficits at this locus, but these values did not depart from null Hardy–Weinberg expectations.

For mt-CO1, we observed a strong cline in the frequency of haplotypes (Figs 1 and 3A). The A haplotypes dominated the two southern sites (Isle of Shoals and Harpswell) and declined in frequency beyond Kent Island at the two sites in Nova Scotia. In fact, this pattern was dominated by the A1 haplotype, which comprised 67%, 84%, 70%, 13% and 7% of the samples, moving from south to north, respectively. Conversely, the B1, B2 and C haplotypes increased in frequency at the two sites in Nova Scotia, so that ∼75% of the samples carried one of these three haplotypes. Interestingly, we sampled a few ‘northern’ B1 and C haplotypes at Isle of Shoals, a southern site that is offshore; however, the sample size was too small (n=6) for more precise estimates of haplotype frequency.

The two nuclear loci showed contrasting clinal patterns (Fig. 3B). The UBE2H A nucleotide, likely derived from the second introduction from a northern European source, increased in frequency moving from south to north, while the SMC T nucleotide, likely present in the original introduction from a south-central European source, reached the highest frequency on Kent Island, but decreased at both the southern and northern ends of the distribution.

In our model, which included only ecological, morphological and sexual variables, only site had a significant effect on the righting response (P=0.0046), while sex approached marginal significance (P=0.149; Table 1). We therefore included both of these terms in the final model, which also had terms for variation at mt-CO1 and nuclear markers. In this second model, sex was marginally significant (P=0.052), site was not significant (P=0.819) and mt-CO1 haplogroup was highly significant (P=0.003; Table 2). Somewhat surprisingly, variation at either nuclear marker had no effect on righting response in this model (P≥0.573). The effects of sex and mt-CO1 terms are graphically presented in Fig. 4. In males, the northern (B) and pan-European (C) haplotypes increased the probability of righting at 4.5°C by ∼20% compared with the southern haplotypes (A) (Fig. 4A). In fact, every single male with a northern or pan-European haplotype was able to successfully right itself at 4.5°C. While the righting response in females did not appear to be affected by mt-CO1 haplogroup, sample sizes were smaller compared with those of males, and were therefore unlikely to uncover small effect sizes, if they indeed exist. The effects of sex and mt-CO1 haplogroup were consistent across the three sites. While sample sizes of some haplogroup × sex combinations were necessarily small because of the clinal structure of mitochondrial variation, males that carried northern or pan-European haplotypes (B,C) benefited from a ∼0.20 increase in the probability of flipping at 4.5°C over males that carried southern haplotypes (A) at both Kent Island and Pomquet (Fig. 4B). Moreover, the effect of carrying an A haplotype was very similar across all three sites: males and females had a ∼0.75 probability of flipping at 4.5°C, regardless of site of collection.

Our results indicate a strong and sex-specific effect of mitochondrial haplotype on the ability of green crabs to tolerate cold stress. In our logistic regression model that included terms for site, sex and genetic factors, mitochondrial haplotype was highly significant and had a large effect on the righting response under cold stress (Table 2, Fig. 4). At each site, males carrying B or C haplotypes had a ∼20% greater probability of righting at 4.5°C than did their counterparts carrying A haplotypes. It appears that differences in mitochondrial genomes play an important role in differences between sites in terms of cold tolerance: in the model without genetic factors, site of origin is highly significant, but when both site and haplogroup are included, site becomes non-significant. We note that our results do not suggest that the mitochondrial CO1 gene specifically plays a role in cold tolerance. All five haplotypes yielded an identical protein sequence, and as mitochondrial genes do not undergo recombination and the entire mitochondrial genome is inherited as a single haploid locus, our mitochondrial marker is linked to mutations throughout the mitochondrial genome that may affect physiological performance.

Nearshore winter temperatures at the Pomquet site were much colder (<0°C) than those at either Kent Island or Harpswell (∼5°C) (Fig. 5). The Pomquet winter temperatures are well below our 4.5°C assay temperature, and suggest there could be strong selection for cold-tolerant phenotypes along the coast of Nova Scotia and northward. The fact that male cold tolerance was connected to mtDNA haplotype while female cold tolerance was not presents an intriguing evolutionary paradox. Because the mitochondrial genome is typically uniparentally inherited through females, mutations arising in the mitochondrial genome that affect males, but not females, are ‘invisible’ to natural selection (Beekman et al., 2014). Two caveats to this general result are situations where there is either high inbreeding or opportunities for kin selection (Wade and Brandvain, 2009). In C. maenas, where the potential for long-distance larval dispersal will likely preclude any substantial inbreeding or kin selection, we see two potential solutions to this paradox. The first posits strong mito-nuclear co-evolution in response to spatially varying temperature stress. Mitochondrial genes are either involved in mtDNA translation or produce proteins needed for the electron transport chain (Ballard and Whitlock, 2004). However, while there are fewer than 40 mitochondrial genes, there are roughly 1500 nuclear genes involved in OXPHOS (Gershoni et al., 2009). Thus, OXPHOS is likely to be strongly dependent on mito-nuclear interactions (Wolff et al., 2014). These protein interactions typically rely on a ‘lock-and-key’ principle, in which the nuclear and mitochondrial proteins must match each other. Support for the general idea that the evolution of stress tolerance requires substantial co-evolution between nuclear and mitochondrial loci comes from a handful of experimental insect studies that have crossed mtDNA variants into divergent nuclear backgrounds (Arnqvist et al., 2010; Koevoets et al., 2012; Hoekstra et al., 2013). These studies have shown that metabolic and life-history traits in F2 hybrids have increased temperature sensitivity when compared with non-recombinant parental lines. In contrast, F2 hybrids between divergent populations of the marine copepod Tigriopus californicus have shown patterns of transgressive segregation and novel thermal phenotypes that are inconsistent with mito-nuclear co-evolution (Willett, 2010; Pereira et al., 2014). A second mechanism that would facilitate selection on male traits caused by mtDNA variation are pleiotropic interactions between male and female traits that affect fitness. Camus et al. (2015) have shown strong, antagonistic pleiotropic effects of mtDNA mutations on male fertility and female longevity in Drosophila melanogaster, illustrating the importance of documenting the role of mtDNA variation in a broad array of male and female traits to understand how selective gradients are shaping mitochondrial and nuclear variation.

As discussed in the Introduction, there are a handful of other examples of species where mitochondrial variation plays an important role in thermal phenotype. While the majority of these examples are mammalian, there is also accumulating correlative evidence from ectotherms (including insects and fishes) that shows clinal variation in mitochondrial haplotypes maps onto the thermal environment, or that molecular signatures of natural selection are accelerated in comparisons between populations living in different thermal environments (reviewed in Camus et al., 2017). To our knowledge, our study is the first to establish a direct link between mitochondrial variation and temperature tolerance in a marine ectotherm. Copepods are a likely second candidate. Harada et al. (2019) have shown that variation in mitochondrial function can explain differences in tolerance to acute heat stress among populations distributed across a temperature gradient spanning coastal California in the intertidal copepod Tigriopus californicus. While mtDNA was not sequenced in this study, these same populations have previously been shown to have between 9.5% and 26.5% sequence divergence among populations. Given that temperature tolerance varies among populations and closely related species that partition distinct thermal environments in the intertidal zone (Somero, 2002), we predict a broader role for mitochondrial variation in thermal adaptation in the marine realm.

A genotypic effect on phenotype was not seen for the two nuclear genes chosen in this study, indicating that these particular loci are not linked to genes involved in the righting response in the cold or that our sample sizes were not large enough to detect smaller effect sizes. However, both nuclear markers were outliers in analyses of transcriptomic variation between northern and southern sites, and SMC was specifically associated with site temperature (Tepolt and Palumbi, 2015). An independent restriction site-associated DNA sequencing (RAD-seq)-based study has also uncovered evidence for potential selection in response to winter temperatures in east coast C. maenas populations (Jeffery et al., 2018) Given this, we suggest that future studies should target a wider range of nuclear markers and physiological tests to more fully explore the role of the nuclear genome in setting thermal limits in the system, particularly in the context of mitochondrial lineage-associated differences in cold tolerance.

The main mechanistic hypothesis for how the mitochondrial genome impacts cold tolerance in ectotherms is oxygen and capacity-limited thermal tolerance (OCLTT) (Pörtner and Knust, 2007; Ballard and Melvin, 2010; Ern et al., 2015). According to this hypothesis, cold-adapted ectotherms have a higher mitochondrial capacity but also an increased metabolic rate. Increased cold tolerance is a trade-off, with lower heat tolerance as its result. Direct evidence for a OCLTT trade-off comes from a mtDNA cline in Australian Drosophila melanogaster, where laboratory populations generated from isofemale lines sampled from the tropics (Queensland) carried mitochondrial haplotypes with greater heat tolerance and less cold tolerance than mitochondrial haplotypes originating from the temperate location (Brisbane), which had greater cold tolerance and less heat tolerance (Camus et al., 2017). There is some support for a similar OCLTT trade-off between heat and cold tolerance in C. maenas populations sampled from different geographical regions (Tepolt and Somero, 2014). Within their native range and the east coast of North America, northern populations dominated by B and C haplogroups are more cold tolerant than southern populations dominated by the A haplogroup (Roman and Palumbi, 2004; Tepolt and Somero, 2014). In turn, the southern populations are more heat tolerant than the northern populations. This apparent trade-off between heat and cold tolerance provides preliminary evidence that mitochondrial adaptation in C. maenas may follow the OCLTT model of thermal tolerance. However, the strongest test of the model and its predicted trade-off will be to experimentally subject crabs carrying different mitochondrial haplotypes to both cold and heat stress.

Both of our nuclear markers differ in their clinal patterns from CO1. Previous sampling revealed that UBE2H variation, which likely derives from the second C. maenas introduction, changed in allele frequency further south than CO1, near Penobscot Bay (Voss, 2016) and that the frequency of the T nucleotide does not increase to >0.50 at our two sites in Nova Scotia. In contrast, the SMC T nucleotide, deriving from the first C. maenas introduction, reaches its highest frequency of 0.33 at Kent Island, but declines in frequency at the more southern and northern sites. Multiple genome-wide studies have shown complex patterns of genomic admixture along the east coast. Initial work with microsatellites in 2007 suggested that mitochondrial DNA was introgressing faster both up and down the coast than nuclear DNA, with ‘southern’ haplotypes found with ‘northern’ nuclear backgrounds outside of the primary introgression zone, and vice versa (‘northern’ haplotypes in ‘southern’ backgrounds) (Darling et al., 2014). In 2011, transcriptome sequencing showed extensive admixture in the recent Newfoundland expansion, which is believed to derive from an admixed Nova Scotia population, with no evidence for any nuclear introgression in ‘southern’ populations in New Jersey and central Maine (Blakeslee et al., 2010; Tepolt and Palumbi, 2015). A geographically broad RAD-seq-based study showed extensive nuclear introgression in southeastern Nova Scotia and southeastern Newfoundland, with most other east coast populations composed primarily or exclusively of animals from a primarily northern or southern nuclear background (Jeffery et al., 2017). Given both this prior work and the nuclear patterns we see here, it is likely that our Harpswell crabs had primarily southern nuclear backgrounds; the Kent Island crabs had admixed nuclear backgrounds; and the Pomquet crabs had northern backgrounds. Thus, the association we observed between cold tolerance and mitochondrial lineage holds true across very different nuclear genetic backgrounds (Fig. 4B). While these mito-nuclear patterns might at first appear to reject the hypothesis of tightly coupled mito-nuclear interactions, to our knowledge none of the previous studies have explicitly tested for linkage between mitochondrial haplotypes and large panels of nuclear SNPs, and it seems likely that these associations would be overwhelmed by the dominant signal from non-linked markers. Future work in extensively admixed populations in southeastern Nova Scotia and Newfoundland that targets genes involved in OXPHOS will help to more fully explore the role of the nuclear genome in mediating the observed effect of mitochondrial lineage on cold tolerance.

The C. maenas hybrid zone presents an interesting system in which to study the role of sexual selection in shaping a variety of traits. Sexual dimorphism in morphology (males reach much larger sizes than females) and behavior is thought to result from male–male conflicts and sexual selection (Berrill and Arsenault, 1982). In addition, there is strong intra-sexual selection for male aggressive behavior to guard vulnerable females from other males before and after copulation, which may interact with natural selection to drive sexual dimorphism in cold tolerance. Any energetic advantage of substitutions in the mitochondrial genome, coupled with complementary nuclear variation, might amplify this behavioral difference between sexes. Mito-nuclear co-evolution will facilitate inter-sexual selection because it requires reliable matching of the mitochondrial genes with nuclear genes (Hill, 2015). Under this paradigm, female preference will evolve in concert with male signals that increase mitochondrial–nuclear matching. Areas where there is extensive natural mixing of mtDNA haplotypes and nuclear backgrounds provide an interesting field test of this hypothesis. While the genetic basis of sex determination in C. maenas is unknown, studies in a few other brachyuran decapods have documented both XY (in Charybdis feriata, also a portunid) and ZW (in Erochier sinensis, a distantly related varunid) systems (reviewed by Chandler et al., 2018). An intriguing possibility emerges in light of the findings that nuclear-encoded genes that function in the mitochondria (N-mt genes) are under-represented on the X chromosome and over-represented on autosomes in XY systems (Drown et al., 2012), but that N-mt genes whose expression has been shown to be sensitive to mtDNA variation are over-represented on the Y chromosome of D. melanogaster (Rogell et al., 2014). These latter genes on the Y chromosome are predicted to be responding to mito-nuclear antagonisms that result from maternal inheritance, by compensating for male-harming mtDNA mutations. Could we have experimentally uncovered an ecologically important male phenotype that has roots in sexually antagonistic evolution? Clearly, sex determination, the location of N-mt genes in C. maenas and their role in variable cold tolerance are all of considerable interest in this hybrid zone.

Previous studies of the dynamics of the green crab hybrid zone have focused on a variety of hypotheses that assume mitochondrial CO1 is a neutral proxy for dispersal and is not under selection (Roman, 2006; Pringle et al., 2011; Darling et al., 2014). More broadly, the field of phylogeography was founded on studies that relied on mitochondrial variation and the assumption of effective neutrality (Avise, 2004). We have established a clear link between mitochondrial haplotype and male cold tolerance in C. maenas, which suggests that differences in cline structure between mitochondrial markers and other, putatively neutral, nuclear markers (e.g. microsatellites and SNPs; Darling et al., 2014; Jeffery et al., 2017) are likely related to a shifting balance between selection and stochastic demographic processes that is marker dependent. With regards to gradients in temperature in the Gulf of Maine and the Canadian Maritimes, the influence of the cold Labrador current is particularly strong as it moves southward along the coast of Nova Scotia and into the upper reaches of the Gulf of Maine, essentially splitting the eastern and western Gulf into two different thermal regimes (Pettigrew et al., 2005). Thus, the thermal transition from the stable and colder waters of the eastern Gulf of Maine to the seasonally warmer temperatures in the western Gulf of Maine appears to play a key role in the structure and persistence of the mitochondrial cline, and we predict similar clines in nuclear loci that play key roles in mitochondrial functions.

Further, the functional differences among mtDNA haplotypes imply that adaptation has played a key role in the limitation of the original C. maenas invasion to the Gulf of Maine, with the species' expansion into the Canadian Maritime provinces coming only after the introduction of presumably cold-adapted genetic variants (Roman, 2006). Assuming mtDNA differences and their likely co-adapted nuclear-encoded mitochondrial genes represent standing genetic variation originally shaped by the physical differences between northern and southern Europe in its native range over the last 500,000 years (Roman and Palumbi, 2004; Tepolt and Somero, 2014; Jeffery et al., 2018), then the success of C. maenas in North America and elsewhere may be related to pre-adaptation to specific thermal environments. Observational studies have suggested a role for mitochondrial preadaptation in invasion success in other systems, notably the bryozoan Watersipora subtorquata (Mackie et al., 2012), but ours presents the first mechanistic support for this possibility. Predicting the success of new invasions may therefore depend on the invading genotypes as well as their interactions with the native community. Furthermore, such preadaptation would have important implications for invasion management, suggesting that there may be advantages to preventing further introductions of invasive species that have already become established in a new region.

Our finding of an association between mtDNA haplotypes and thermal tolerance in invasive green crabs may help to explain the discordance between mtDNA and presumably neutral nuclear markers along the northeastern coast of the USA. The increase in frequency of cold-adapted mitochondrial haplotypes along the coast of Nova Scotia and within the Gulf of St Lawrence coincides with a shift towards greater seasonal extremes in winter temperatures. This finding adds to a growing number of marine and terrestrial studies indicating an important role for mtDNA variation in shaping organismal responses to temperature stress, but also suggests complexity in the genomic architecture of stress-related phenotypes. The male-specific effect of mtDNA variation we observed can be viewed as an evolutionary paradox, whose resolution may reside in mito-nuclear co-adaptation to thermal stress, or in pleiotropic effects of mtDNA variation on male and female traits. With the advent of a rich and increasingly affordable array of genomic tools for natural systems (Sherman et al., 2016), C. maenas is an excellent model with which to explore the genomic drivers of a wide variety of physiological, population and ecosystem processes in a super-invader that is profoundly changing benthic seascapes throughout the world. This work highlights the importance of understanding the genetic makeup of invading populations and how these invading genotypes affect ecological phenotypes when predicting the success and ecological outcome of invasions.

We would like to thank T. Suskiewicz and L. Johnson for help collecting crabs from Halifax, NS, and Robin Seeley for collections from the Isle of Shoals. We thank Timothy Fuller for designing and testing the SMC primers. We thank Mark Murray for facilitating a productive stay on Kent Island, and Nick Keeney for assistance with animal care at the Schiller Coastal Studies Center. This is publication no. 5 from the Bowdoin Marine Laboratory.