Targeted approach to identify genetic loci associated with evolved dioxin tolerance in Atlantic Killifish (Fundulus heteroclitus)
BMC Evolutionary Biology volume 14, Article number: 7 (2014)
The most toxic aromatic hydrocarbon pollutants are categorized as dioxin-like compounds (DLCs) to which extreme tolerance has evolved independently and contemporaneously in (at least) four populations of Atlantic killifish (Fundulus heteroclitus). Surprisingly, the magnitude and phenotype of DLC tolerance is similar among these killifish populations that have adapted to varied, but highly aromatic hydrocarbon-contaminated urban/industrialized estuaries of the US Atlantic coast. Multiple tolerant and neighboring sensitive killifish populations were compared with the expectation that genetic loci associated with DLC tolerance would be revealed.
Since the aryl hydrocarbon receptor (AHR) pathway partly or fully mediates DLC toxicity in vertebrates, single nucleotide polymorphisms (SNPs) from 42 genes associated with the AHR pathway were identified to serve as targeted markers. Wild fish (N = 36/37) from four highly tolerant killifish populations and four nearby sensitive populations were genotyped using 59 SNP markers. Similar to other killifish population genetic analyses, strong genetic differentiation among populations was detected, consistent with isolation by distance models. When DLC-sensitive populations were pooled and compared to pooled DLC-tolerant populations, multi-locus analyses did not distinguish the two groups. However, pairwise comparisons of nearby tolerant and sensitive populations revealed high differentiation among sensitive and tolerant populations at these specific loci: AHR 1 and 2, cathepsin Z, the cytochrome P450s (CYP1A and 3A30), and the NADH dehydrogenase subunits. In addition, significant shifts in minor allele frequency were observed at AHR2 and CYP1A loci across most sensitive/tolerant pairs, but only AHR2 exhibited shifts in the same direction across all pairs.
The observed differences in allelic composition at the AHR2 and CYP1A SNP loci were identified as significant among paired sensitive/tolerant populations of Atlantic killifish with multiple statistical tests. The genetic patterns reported here lend support to the argument that AHR2 and CYP1A play a role in the adaptive response to extreme DLC contamination. Additional functional assays are required to isolate the exact mechanism of DLC tolerance.
Dioxin-like contaminants (DLCs), such as some polychlorinated biphenyls (PCBs), are highly toxic aromatic hydrocarbon pollutants whose ubiquitous occurrence presents global ecological and human health risks. The early life stages of fish are particularly sensitive to these toxic DLCs, and the Atlantic killifish, Fundulus heteroclitus, is one of the more sensitive fish species . Despite this species’ relative sensitivity to DLC exposure, several wild killifish populations residing in heavily contaminated North American Atlantic coast estuaries have recently and independently evolved dramatic, heritable, and adaptive tolerance to DLCs , for which the mechanistic basis has yet to be fully explained. To address this issue, a targeted, candidate gene scan was performed to reveal genetic variation associated with tolerance in four wild DLC-adapted killifish populations.
Considerable effort has been spent attempting to identify the genetic and biochemical mechanisms underlying inter- and intra-specific variation in DLC sensitivity in vertebrates. Multiple lines of evidence support the crucial role of the aryl hydrocarbon receptor (AHR) pathway in DLC toxicity in mammals. Polymorphisms in the ligand binding domain of the AHR among mouse strains result in differences in ligand binding affinity, and low binding affinity appears to protect against all toxic responses to DLCs. In rats, DLC tolerance is associated with variation in the transactivation domain of the AHR, yet functional consequences of the variation are less predictable . A candidate gene approach effectively identified two amino acid substitutions in the AHR among avian wildlife and a consistent relationship between the amino acid residues present and DLC sensitivity was observed at the species level [4, 5]. In fish, the striking difference in DLC sensitivity among wild Atlantic tomcod populations has been attributed to a six base pair deletion in the AHR2 gene that results in a five-fold decrease in ligand binding affinity and reduced ability to promote expression of detoxification enzymes targeted by the AHR pathway .
Within the killifish AHR pathway, several non-synonymous single nucleotide polymorphisms (SNPs) have been identified in two AH receptor genes (AHR1 and AHR2), but patterns of genetic variation at these loci do not unequivocally reflect differences in DLC sensitivity among populations and no functional consequences (i.e., ligand binding and the ability to interact with xenobiotic metabolizing enzymes) were associated with AHR1 variants [7–9]. As an alternative to the candidate gene approach, Williams and Oleksiak [10, 11] performed whole genome scans, whereby patterns of variation across hundreds of genetic loci were contrasted between killifish populations from polluted and reference sites in order to identify genes under selection with respect to DLC contamination. A handful of selectively important genetic markers were identified in each of three separate comparisons between populations residing in polluted habitats and their respective reference populations, and a single marker (in the CYP1A promoter) was found to be selectively important in all comparisons.
Both the single candidate gene approach and genome-wide scans (typically with sets of anonymous genetic markers) have led to great success in elucidating the genetic basis for many adaptive phenotypes [12, 13], but neither has offered a comprehensive link to the observed variation in DLC sensitivity among wild killifish populations. A ‘candidate gene scan’ approach, which targets a relatively large set of expressed genes with known physiological function, should increase the probability of isolating genes that are under selection and relevant with respect to traits of interest . Thus, this approach was adopted to maximize efficiency in the identification of genes associated with DLC tolerance, and complement previous and ongoing work investigating the mechanism(s) involved in the repeated adaptation to DLCs in wild killifish populations.
In this study, patterns of genetic variation at SNP markers distributed across genes that are components of, or whose expression is affected by the AHR and interacting pathways were examined among eight killifish populations that vary genetically in their responsiveness to DLCs. Four relatively uncontaminated populations whose sensitivities to a prototypical DLC (3,3’,4,4’,5- pentachlorobiphenyl, PCB126) range from 20–38 ng/L (reviewed in ) were chosen as sensitive populations. Each of these populations is located near one of four EPA-designated urban/industrial estuarine Superfund sites. Killifish resident to these sites are dramatically tolerant to the effects of PCB126, displaying LC20 values ranging from ~400 to ~ 8000 times higher than those for the sensitive killifish (Table 1). A companion study  provides a fine-scaled examination of genetic variation in three AHR-related loci (AHR1, AHR2, and AHRR) among killifish populations residing in uncontaminated and polluted habitats of the North Atlantic coast of the US, including some of the same populations examined in this study.
The question of whether the genetic variants observed at the targeted SNP markers in this study could explain the stark phenotypic differences between DLC-adapted (−tolerant) and -sensitive killifish populations was addressed in several ways. Measures of genetic diversity were compared between DLC-sensitive and DLC- tolerant groups and patterns of genetic differentiation among populations were tested against expectations of isolation by distance. In addition, by contrasting the behavior of individual loci among DLC-sensitive and DLC-tolerant populations, specific loci under selection were identified.
SNP marker development and preliminary screening
An extensive literature search was conducted to identify genes and biochemical pathways with demonstrated and potential involvement in the toxic responses to DLCs. A list of over 150 genes was compiled, which included components of the AHR pathway, nuclear receptors known to ‘cross talk’ with the AHR pathway (i.e., estrogen receptors, retinoic acid receptors, hypoxia inducible factors), cytochrome p450s, genes involved in cardiac development, cathepsins, and genes having oxidoreductase activity (e.g. , and references therein). The gene list was filtered at several stages of the marker development process. Sequence information for many, but not all, of the genes listed in Additional file 1 was retrieved from the F. heteroclitus unigene database in GenBank (http://www.ncib.nlm.nih.go/GenBank/). Additional unpublished sequences were kindly provided by Sibel Karchner (personal communication). Putative SNPs were detected with the QualitySNP pipeline . In QualitySNP, sequences with > 95% similarity were assembled into contiguous sequences (contigs) using the sequence assembly program CAP3. Available contigs, those containing ≥ 4 overlapping sequences, were then evaluated for polymorphisms. Available contigs with SNPs were subjected to further scrutiny to assess the reliability of the putative SNPs identified. Specifically, SNPs were considered ‘true’ if they were represented in ≥ 2 ESTs, paralogous sequences within the cluster could be distinguished and haplotypes identified, and they were located in high quality sequence regions . PCR primers and melting temperature (Tm)-shift genotyping assays  were then designed for the ‘true’ SNPs with suitable flanking sequence (e.g. no SNPs, low sequence complimentarity in priming sites). A test panel of eight killifish from two populations (BI and NBH) was genotyped at 128 of the loci deemed highly reliable by QualitySNP. A SNP maker was considered valid if amplification product of the appropriate size was generated and polymorphism was observed among the eight individuals in the test panel. Outcomes for each stage of the SNP marker development process are detailed in Additional file 1.
During the summer and fall months between 2008 and 2011, adult Fundulus heteroclitus were collected using baited minnow traps from eight estuarine sites spanning approximately 600 km of the Atlantic Coast of the USA (Figure 1). These specific killifish populations had already been characterized as either DLC-tolerant or DLC-sensitive, based on early life stage sensitivity to PCB126 [2, 18, 19]. To better assess genetic differences between tolerant and sensitive populations absent geographic influences, each DLC-tolerant population was paired with a nearby DLC-sensitive population. These pairs are located within the same or adjacent states but are separated by discontinuous shoreline and deep channels that limit fish migration and contaminant spread. Latitudinal distances between killifish residence sites (Table 1) were used to provide a consistent proxy to assess regional influences but were not intended to convey ‘as-the-fish-swims’ distances. Sixty to 100 live or sacrificed/frozen fish from each population were transported to the US Environmental Protection Agency (EPA) laboratory in Narragansett, RI, USA. Those that were transported live were sacrificed immediately upon return to the laboratory. Whole fish were stored at either −20 or −80°C prior to DNA extraction.
Sample preparation and population genotyping
Genomic DNA was extracted from approximately 20 mg of caudal fin tissue from 36 or 37 of the 60–100 archived individuals per population according to the QIAGEN DNeasy protocol for animal tissue (optional RNase treatment included), quantified with the PicoGreen dsDNA assay (Invitrogen), and diluted to a standard concentration of 20 ng/μl. Diluted DNA extracts were submitted to the University of Minnesota’s BioMedical Genomics Center in a 96-well format for sample quality assessment and SNP genotyping using the Sequenom MassARRAY® technology. Three multiplex assays (containing 32, 27, and 12 SNPs respectively) were designed using MassARRAY® Designer software. Sampled F. heteroclitus for which DNA was extracted were genotyped with the first two plexes (because all SNP containing genes were represented) following the iPLEX assay protocol. Genotypes for each individual at each locus were called using the Sequenom System Typer Analysis package.
Routine population genetics analyses were conducted using the freely available software package Arlequin ver. 3.5 . Standard diversity indices, including the percentage of polymorphic loci (PO), average observed and expected heterozygosities (HO and HE), and the within population fixation index (FIS), were determined for each sampled population. Indices were based on data from loci with < 10% missing data. Loci were also tested for departures from Hardy-Weinberg equilibrium (HWE) with 100,000 permutations and the percentage of loci in HWE was calculated for each population. The ability of this suite of SNP markers to detect genetic differentiation among populations was assessed by computing pairwise multi-locus FST estimates. To test whether the assumption of independence among loci would be violated by including multiple SNPs per gene in the analysis, data from a subset of SNPs (representing one SNP per gene) were also analyzed. Results did not differ significantly between the complete and limited datasets; therefore, only results from the complete dataset are reported here.
Genetic structure, where populations were assigned to groups defined by DLC sensitivity, was also included in the analysis and the pairwise multi-locus FST between the two groups was estimated with the ‘Compute pairwise FST’ option in Arlequin ver. 3.5. An analysis of molecular variance (AMOVA) was performed to better understand how the observed genetic variance was partitioned within populations, among populations within each group, and between groups. In addition, as in , Student’s T-tests were used to detect significant differences in genetic diversity measures (PO and HO) among the two groups.
For species like F. heteroclitus that are non-migratory, genetic differences at neutral loci should accumulate over time and generate a pattern of isolation by distance (IBD), whereby differentiation among populations increases with geographic distance . Deviations from IBD patterns can be attributed to responses to local selection pressures . Mantel tests and reduced major axis (RMA) regression analyses were performed in Isolation By Distance Web Service (IBDWS)  with 10,000 permutations to test for significant correlations between genetic distance (pairwise multi-locus FST values) and latitudinal distance. Mantel tests were also used to test for significant correlations between pairwise genetic differences and relative differences in sensitivity to PCB126 (LC20 as reported in ); moreover, partial Mantel tests were conducted to determine if there was a significant relationship between sensitivity to PCB126 and genetic divergence after taking latitudinal distance into account.
Although most of the genes included in this analysis were chosen based on prior work suggesting they might be responsive to DLCs, whether any or all contribute to the adaptive phenotype remains largely unresolved. To identify selectively important SNPs among the genes surveyed, FST values were calculated at each locus separately for each sensitive/tolerant comparison with the AMOVA function in Arlequin ver. 3.5. The statistical significance of each FST was determined through permutation testing with 10,000 iterations. An FST modeling approach similar to that described in  as implemented in Arlequin ver. 3.5 was also used to detect outliers. FST distributions conforming to a neutral model were simulated with 20,000 iterations, heterozygosities computed from empirical data, and assuming 10 demes. SNPs with FST values exceeding the 95th percentile of the null distribution were considered to be strong candidates for natural selection. In addition, as in , minor allele frequencies were calculated at each locus for each population with the expectation that loci associated with the adaptive phenotype would display shifts in allele frequency for each sensitive/tolerant pair. Consistent shifts in the same direction across all four comparisons were considered further evidence for selection acting at the loci.
Killifish sequences representing 125 of the 150 candidate genes were mined from the GenBank nucleotide and unigene databases. CAP3 assembled those sequences into 183 contigs, of which 105 were available for further analysis (met the criterion of containing four or more overlapping sequences). Approximately 500 highly reliable putative SNP loci distributed among 50 genes were identified with the QualitySNP pipeline. Tm-shift genotyping assays were designed for 128 of the highly reliable putative SNPs and the validity of each locus was tested by genotyping eight killifish collected from two populations (BI and NBH) (Additional file 1). Twenty-four of the 128 putative SNPs assayed (19%) did not amplify or produced uninterpretable melting curves. An additional 29 of the SNPs tested (23%) appeared to be monomorphic. Ultimately, 75 (59%) of the putative SNPs screened were polymorphic and 59 (those represented in the first two iplex assays) were used in the population genetic analysis (Table 2).
Standard diversity indices
Substantial genetic variation was observed in each killifish population examined. The percentage of polymorphic loci (PO) ranged from 51% to 73% across populations, with a greater percentage of monomorphic loci in NBH and the Virginia populations. A similar pattern was detected when observed heterozygosity (HO) was used to measure diversity; ER, KC and NBH populations were the least heterozygous. According to the average within population fixation index (FIS), the loss of heterozygosity in NBH, ER, KC, and SH was statistically significant (Table 3).
A majority of the 59 SNP loci assayed (64%) were considered to be common SNPs, with minor allele frequencies greater than 0.1, when data from all populations were pooled. For each population, the percentage of common SNPs ranged from 54% to 76% with BP and ER at the two extremes. Within each population the minor allele frequency, averaged across all polymorphic loci, varied between 0.16 and 0.22. The mean minor allele frequency was lowest in BP and BI fish and highest in ER fish (Table 3).
Most of the SNP loci in this study were found to be in Hardy-Weinberg equilibrium within each population. Again, the population with the fewest loci conforming to Hardy-Weinberg expectations was ER (Table 3). In total, 14 loci deviated from HWE, but no single deviation was consistent across all eight populations or across the four tolerant populations.
Mean pairwise FST values for all population comparisons were found to be statistically significant and suggest moderate to very high levels of genetic differentiation. Overall, markedly higher genetic differentiation was detected between the two Virginia populations and all others (Figure 2). A standard analysis of molecular variance (AMOVA) confirmed that a substantial proportion (23%) of the observed molecular variation can be attributed to differences among populations.
No significant differentiation emerged when populations were grouped by DLC sensitivity. The expectation that, on average, DLC-adapted populations are less diverse (as quantified by PO and HO) than DLC-sensitive populations, was not supported (Table 4). The genetic differentiation between the two groups as measured by FST was 0.0212; however, the hierarchical AMOVA results suggest that the groups were not significantly different (Table 5).
Isolation by distance
A regression of the genetic and latitudinal distance matrices when all populations were included in the analysis was statistically significant (r =0.8607, p < 0.001) reflecting a pattern of IBD (Figure 3A). Given that populations were sampled from three distinct geographic regions differentially impacted by the Pleistocene glacial retreat, the observed relationship is most likely driven by long-term history and demography rather than contemporary forces . Not surprisingly, the pattern persisted when only sensitive populations were considered (r = 0.9223, p = 0.0372). A clear positive relationship was also apparent when only tolerant populations were included in the analysis (r = 0.9443), but the trend was not significant (p = 0.0856) (Figure 3C). The lack of a significant IBD pattern among DLC-adapted populations could be indicative of local selection (in response to chemical contamination) counteracting the effects of history, demography, migration, and drift, but is more likely a consequence of small sample size. It was hypothesized that if selection is a major force shaping patterns of genetic variation among the sampled killifish populations, a strong correlation between pairwise genetic differences and relative differences in sensitivity to PCB126 should be evident. Additional Mantel tests found no such relationship (data not shown).
Loci under selection
The number of loci with significant FST values greater than 0.10, suggesting they are greatly differentiated among populations , varied among the four sensitive/tolerant population pairings. Of the four comparisons, the largest number of highly differentiated loci was detected between BI and NBH (Table 6). Included among these loci were the aryl hydrocarbon receptors (AHR) 1 and 2, cathepsin Z, the cytochrome P450s (CYP 1A and 3A30), and the NADH ubiquinone oxidoreductase MLRQ subunit. Only the CYP1A locus exhibited a consistently high FST value across all population pairs while both AHR genes were highly differentiated in all comparisons except for that between SH and NWK. While great differences were detected between paired populations at the CYP1A locus, a hierarchical AMOVA did not confirm that those differences are associated with DLC sensitivity (Table 5). In contrast, variation at the AHR2 locus could distinguish DLC-sensitive and tolerant groups (Table 5).
Further exploration of associations between specific loci and DLC sensitivity using an FST modeling approach produced results similar to those derived from the locus-by-locus AMOVA. Again, the CYP1A locus was identified as a significant outlier with respect to the simulated FST null distribution in three out of the four pairwise comparisons between each sensitive population and its DLC-tolerant counterpart. Moreover, when data for all sensitive populations were pooled and compared to the pooled data for all tolerant populations, the AHR2 locus emerged as a significant outlier (Figure 4). These findings lend support to the hypothesis that the CYP1A and AHR2 loci may be involved in the evolution of DLC tolerance.
Evidence of selection can also be gleaned from subtle shifts in allele frequencies in response to environmental variables after controlling for population structure . By reviewing differences in minor allele frequency between sensitive and DLC-adapted population pairs, a strong signal was apparent for two SNPs: AHR2_1929 and CYP1A_2140 (Table 7). Upon further examination, the shift in minor allele frequency for the AHR2 SNP varied in magnitude but was consistently in the same direction for all four population comparisons, providing yet another line of evidence linking AHR2 with the DLC-tolerant phenotype. This was not the case at the CYP1A locus. Although substantial differences in minor allele frequency were observed between each sensitive and tolerant pair, the frequency shift between KC and ER was in the direction opposite of the other three pairings (Figure 5).
Four independent DLC-adapted killifish populations were contrasted with neighboring sensitive counterparts in an attempt to reveal genetic loci associated with intra-specific DLC tolerance in the wild. The AHR pathway is known to mediate toxic responses to DLCs in all vertebrates and several studies involving killifish have implicated altered AHR pathway function in eliciting the tolerant phenotype [29, 30, 37]. Therefore, a ‘candidate gene scan’ approach, focused on components and targets of the AHR and interacting pathways represented among available killifish genetic resources, was applied to gain a better understanding of the genetic basis for the observed phenotypic variation in DLC sensitivity among killifish populations.
Genetic diversity and population structure of tolerant killifish populations were examined with the expectation that patterns of historical stress would be revealed. Levels of genetic diversity and population structure in killifish have been evaluated in several studies; some with the specific aim of addressing the impacts of pollutant-driven selection [9, 22, 38–41]. Results of the multi-locus analyses conducted here are in agreement with previous reports. Comparably high levels of genetic variation were estimated for all populations, rendering the hypothesis that a genetic bottleneck facilitated the emergence of the DLC-adaptive phenotype highly unlikely. With respect to population structure, each population included in this study was moderately to highly different from all others and observed FST values were comparably higher than earlier measures for killifish populations spanning a similar geographic range (e.g. FST ranged from 0.03 - 0.2 in  and from 0.01 - 0.24 in ). The difference in the extent of population structure detected among independent studies could be a consequence of the type of markers used (targeted SNPs vs. putatively neutral microsatellites); however, the overall patterns are consistent: genetic differences appear to be driven by geographic distance (latitudinal or shoreline) rather than DLC sensitivity. Moreover, a hierarchical AMOVA did not attribute any of the existing molecular variation to differences between DLC-sensitive and DLC-tolerant groups.
The failure to detect a significant relationship between multi-locus measures of genetic diversity/differentiation and increased tolerance to DLCs may be because a majority of the markers used in the current and previous analyses are selectively neutral. It has been suggested that neutral markers best reflect the effects of anthropogenically-mediated environmental change when populations are in decline and genetic exchange among populations is restricted . There is no evidence that either sensitive or DLC-tolerant killifish populations have experienced a reduction in population size [22, 38, 39]. A more plausible explanation for the rapid adaptation to DLC contamination is that the trait in question is controlled by a small number of loci [39, 43]. This explanation is consistent with theoretical models that predict single allelic differences of large effect dominate adaptive shifts when environmental change is sudden and selection is intense . The prediction has been tested and verified in insect populations reacting to pesticides , benthic marine invertebrates subject to heavy metal toxicity , and fish exposed to DLCs .
Phenotypic similarities among tolerant populations exposed to a prototypical DLC (e.g. LC20)  supported an expectation of common loci under selection; however, genetic loci associated with tolerance might vary across populations due to differences in the selective agents present in their native habitats. The presumed selective agents include large classes of aromatic hydrocarbons whose toxic effects are mediated fully (DLCs) or partly (polycyclic aromatic hydrocarbons, PAHs) through the AHR pathway. While urban contamination includes moderate levels of both PAHs and PCBs, toxic levels of DLCs have been measured at NBH and NWK, PAHs at ER, and both at BP [2, 46]. Therefore, genetic analyses were conducted to identify common loci under selection in tolerant populations, and differences between sensitive/tolerant paired killifish populations.
Since functional variation in AHR-ligand binding initiates the AHR pathway cascade (and largely explains some intra- and inter-species differences in DLC sensitivity, e.g.  and references therein), evidence for selection acting on AHR loci was of obvious interest in comparing tolerant versus sensitive killifish populations. As in other fish species, killifish possess at least two distinct AH receptor genes, AHR1 and AHR2, and the expression of AHR2 predominates in most tissues . Strong signals of selection were detected for an AHR1 (AHR1_1530) and AHR2 SNP (AHR2_1929) included in this analysis in three of the four pairwise comparisons (SH/NWK excluded) when locus-by-locus FST were considered. Although these two SNPs are synonymous and do not result in amino acid changes, they are both located in the transactivation domain of their respective AHR genes and in close proximity to non-synonymous SNPs found to be under selection in Reitzel et al. . Minor allele frequencies (equivalent to base frequencies referred to in ) of the AHR2_1929 SNP were consistently higher in the DLC-tolerant populations for all four population pairs. Although the FST modeling approach identified the AHR1_1530 SNP as a significant outlier in only one of the comparisons (BI/NBH), when genotype data from all DLC-sensitive populations were pooled and compared to that of all DLC-tolerant populations AHR2_1929 was the only locus found to deviate from neutral expectations. Similar patterns of genetic variation in AHR1 and AHR2 loci is not surprising given that these two genes are arranged in tandem (and therefore linked) within the killifish genome [9, 48]. To determine whether variation in both loci, a single locus, or neither locus underlies the DLC-tolerant phenotype, population genetic data must be accompanied by functional assays. A focused examination of allelic variation in killifish AHR1 and AHR2 revealed SNPs under positive selection in both genes; however, AHR1 variants were not responsible for alterations in receptor function and AHR2 variations have yet to be tested [7, 9]. Studies investigating the genetic basis for DLC-tolerance in zebrafish and Atlantic tomcod have isolated an AHR2 gene as a key player in mediating DLC sensitivity [6, 48]; while AHR1 and AHR2 seem to play functional roles in dioxin toxicity in red seabream . Taken together, these results suggest that AHR2 variation likely plays a strong role in DLC sensitivity and tolerance in killifish, but more complex interactions may be revealed as new AHR paralogs are being identified and characterized .
Consistent with the refractory induction patterns for CYP1A (an early and sensitive marker of AHR pathway activation) observed among DLC-tolerant killifish populations in laboratory studies [2, 29, 49], the CYP1A SNP emerged as a strong candidate for selection in all tolerant populations. High FST values and large differences in minor allele frequencies were observed at the CYP1A locus in all four population comparisons and the FST modeling approach identified CYP1A as a significant outlier in three out of the four comparisons (BI/NBH excluded). However, the shift in allele frequency, although significant, was not in the same direction across all comparisons. A full genome scan analysis of three of the same tolerant killifish populations (excluding BP) and their sensitive counterparts also identified a SNP marker in the CYP1A promoter region as the only locus (out of 354 screened) under selection in all three DLC-tolerant populations surveyed, but again, the direction of the allele frequency shift between DLC-sensitive and tolerant populations was not uniform . In a follow-up study, Williams and Oleksiak  found that CYP1A promoter variants derived from a DLC-tolerant population (NBH) resulted in elevated expression of CYP1A in vitro, contradicting the well-documented refractory response of CYP1A to DLC exposure among tolerant populations in vivo. These previously reported results, coupled with the knowledge that the CYP1A SNP included in this analysis is located in the 3’ untranslated region of the gene  make the exact functional role (if any) of CYP1A SNPs in the DLC-tolerant phenotype unclear. It may be that additional, yet to be discovered factors associated with CYP1A regulation are also important.
The interpretation of CYP1A as a candidate for selection must also take into account that tolerance has evolved in response to different AHR ligands. Specifically, the role of CYP1A in AHR-mediated toxicity varies by ligand class. Classically described, AHR agonists induce the production of enzymes (predominantly CYP1A) that metabolize PAHs, but not DLCs. Unlike DLC toxicity, PAH toxicity is self-limiting, due to AHR-enhanced PAH elimination, and includes components that are not AHR-mediated. For example, CYP1A knockdown studies in zebrafish embryos have demonstrated the protective value of CYP1A during developmental PAH exposures . That tolerant killifish populations from varied selective environments (PAH- versus DLC-dominated pollution profiles) show similar, poorly-inducible CYP1A phenotypes warrants the consideration of the adaptive value of the CYP1A recalcitrant phenotype more broadly, e.g., as an energy conservation strategy associated with chronic pollution exposures. As in other species [53, 54], variation in the CYP1A sequence among killifish populations may be beneficial if associated with conditional fitness under transient, low level PAH exposures. Alternatively, variation in CYP1A may be related to its position ‘downstream’ in the AHR pathway, and secondary to changes in ‘upstream’ loci, causally associated with tolerance. Regardless of mechanism, variation at the CYP1A locus differentiates tolerant from sensitive killifish.
Among nearby sensitive/tolerant killifish population comparisons, the BI/NBH pairing was the most genetically differentiated based on multi-locus FST; many additional loci were identified as candidates for selection (i.e., FST > 0.1 and significant allele frequency shifts). In addition to the AHRs and CYP1A, these populations appear to be highly differentiated at cathepsin F, cathepsin Z, CYP3A30, and the NADH ubiquinone oxidoreductase MLRQ subunit loci. Cathepsins are a large group of proteolytic enzymes that have been implicated in cardiomyopathies and cardiovascular disease, ultimately resulting in impaired pump function [55, 56]. Given that the cardiovascular system is a main target of DLC toxicity in all vertebrates , it is reasonable to propose that alterations in the cathepsin coding sequence could contribute to existing differences in DLC sensitivity. CYP3A30 is an abundant xenobiotic metabolizing enzyme in killifish livers responsible for the breakdown and clearance of a wide array of anthropogenically derived pollutants  and has been identified as a target of the AHR pathway . In killifish, expression of this gene was found to be significantly higher in field caught ER females relative to those collected from KC ; however, it was not differentially expressed in killifish embryos derived from the same DLC-tolerant and sensitive populations included in this study when exposed to PCB126 under controlled laboratory conditions [29, 37]. Mutations in the NADH ubiquinone oxidoreductase MLRQ subunit are associated with metabolic diseases . Again, expression of this gene was found to be significantly higher in field caught ER females relative to those collected from KC . Moreover, expression of another component of the NADH ubiquinone oxidoreductase enzyme, NDUB2, was found to be higher in the brains of NBH, NWK, and ER adult fish . The repeated association of NADH ubiquinone oxidoreductase with DLC-tolerance suggests that biochemical pathways other than the AHR could be involved in this chemically-induced stress response.
It is not known whether the identification of additional candidate loci for strong selection in the BI/NBH pair only suggests unique tolerance-related or biologically-relevant differences (BI is the only oceanic rather than estuarine site included in this analysis) or technical artifacts (statistical power). While similar phenotypes among the four tolerant populations examined suggest a conserved biochemical basis for intra-specific DLC tolerance in killifish, whether that similarity is constrained to identical nucleotide changes remains to be seen. The genetic mechanisms of adaptation could vary among DLC-tolerant populations. Alternatively, the differences detected may be a reflection of additional unique stressors encountered by each population pair.
The purpose of this study was to identify genetic polymorphisms associated with DLC sensitivity in Atlantic killifish from a suite of candidate loci involved in the AHR and interacting pathways; whether the polymorphisms in and of themselves are responsible for the drastic differences in DLC sensitivity among populations was beyond the scope of this work. Two loci, AHR2 and CYP1A, displayed patterns of variation consistent with selection in each of four pairwise comparisons. Additional loci with specific alleles significantly overrepresented in some, but not all, DLC tolerant populations underscore the possibility that the genetic variation in each population may have been shaped by similar yet unique selection pressures. Although the intention was to include all components of the AHR and interacting pathways in this population screen, the underrepresentation or absence of key genes (e.g., aryl hydrocarbon receptor nuclear translocator (ARNT)) in the F. heteroclitus unigene database resulted in an enriched but incomplete set of candidate genetic markers for DLC toxicity and (presumably) tolerance. Genetic resources being made available through the Fundulus genome consortium as well as a parallel, unbiased, Quantitative Trait Locus (QTL) approach to the discovery of the genetic mechanism of DLC tolerance in killifish will greatly increase our understanding of this dramatic example of anthropogenically induced, rapid adaptation in the wild.
Availability of supporting data
The sequences for the SNPs analyzed in this article are available through the Dryad Digital Depository and can be accessed by searching for the data package title Data from: Targeted approach to identify genetic loci associated with evolved dioxin tolerance in Atlantic Killifish (Fundulus heteroclitus), doi:10.5061/dryad.2355t Data files: Proestou_etal_Killifish_SNP_sequences.
Aryl hydrocarbon receptor
Polycyclic aromatic hydrocarbon
New bedford harbor
Van Veld PA, Nacci D: Chemical tolerance: Aclimation and adaptations to chemical stress. The Toxicology of Fishes. Edited by: Diguilio RT, Hinton DE. 2008, Washington, DC: Taylor & Francis, 597-644.
Nacci DE, Champlin D, Jayaraman S: Adaptation of the Estuarine Fish Fundulus heteroclitus (Atlantic Killifish) to Polychlorinated Biphenyls (PCBs). Estuaries Coast. 2010, 33: 853-864. 10.1007/s12237-009-9257-6.
Pohjanvirta R, Korkalainen M, Moffat ID, Boutros PC, Okey AB: Role of the AHR and its structure in TCDD toxicity. The AH Receptor in Biology and Toxicology. Edited by: Pohjanvirta R. 2012, New Jersey: John Wiley & Sons, 181-196. 1
Head JA, Hahn ME, Kennedy SW: Key amino acids in the aryl hydrocarbon receptor predict dioxin sensitivity in avian species. Environ Sci Technol. 2008, 42: 7535-7541. 10.1021/es801082a.
Farmahin R, Manning GE, Crump D, Wu D, Lukas JM, Jones SP, Hahn ME, Karchner SI, Giesy JP, Bursian SJ, Zwiernik SJ, Frederick TB, Kennedy SW: Amino acid sequence of the ligand-binding domain of the aryl hydrocarbon receptor 1 predicts sensitivity of wild birds to effects of dioxin-like compounds. Toxicol Sci. 2013, 131: 139-152. 10.1093/toxsci/kfs259.
Wirgin I, Roy NK, Loftus M, Chambers RC, Franks DG, Hahn ME: Mechanistic basis of resistance to PCBs in Atlantic tomcod from the Hudson river. Science. 2011, 331: 1322-1325. 10.1126/science.1197296.
Hahn ME, Karchner SI, Franks DG, Merson RR: Aryl hydrocarbon receptor polymorphisms and dioxin resistance in Atlantic killifish (Fundulus heteroclitus). Pharmacogenetics. 2004, 14: 1-13. 10.1097/00008571-200401000-00001.
Hahn ME, Karchner SI, Franks DG, Evans BR, Nacci D, Champlin D, Cohen S: Mechanism of PCB-and Dioxin-Resistance in Fish in the Hudson River Estuary: Role of Receptor Polymorphisms. Final Report, Hudson River Foundation Grant, Final Report, Hudson River Foundation Grant 004/02Ahttp://hudsonriver.org/ls/ 2005
Reitzel AM, Karchner SI, Franks DG, Evans BR, Nacci D, Champlin D, Vieira VM, Hahn M: Genetic diversity at aryl hydrocarbon receptor loci in PCB sensitive and PCB resistant populations of Atlantic killifish (Fundulus heteroclitus). BMC Evol Biol. 2014, 14: 6-10.1186/1471-2148-14-6.
Williams LM, Oleksiak MF: Signatures of selection in natural populations adapted to chronic pollution. BMC Evol Biol. 2008, 8: 282-10.1186/1471-2148-8-282.
Williams LM, Oleksiak MF: Ecologically and evolutionarily important SNPs identified in natural populations. Mol Biol Evol. 2011, 28: 1817-1826. 10.1093/molbev/msr004.
Rosenblum EB, Hoekstra HE, Nachman MW: Adaptive reptile color variation and the evolution of the MClR gene. Evolution. 2004, 58: 1794-1808.
Jaquiery J, Stoeckal S, Nouhaud P, Mieuzet L, Mahéo F, Legeai F, Bernard N, Bonvoisin A, Vitalis R, Simon JC: Genome scans reveal candidate regions involved in the adaptation to host plant in the pea aphid complex. Mol Ecol. 2012, 21: 5251-5264. 10.1111/mec.12048.
Hemmer-Hansen J, Nielsen EE, Meldrup D, Mittelholzer C: Identification of single nucleotide polymorphisms in candidate genes for growth and reproduction in a nonmodel organism; the Atlantic cod Gadus morhua. Mol Ecol Resour. 2011, 11 (Suppl. 1): 71-80.
Denison MS, Soshilov AA, He G, DeGroot DE, Zhao B: Exactly the same but different: promiscuity and diversity in the molecular mechanisms of action of the Aryl Hydrocarbon (Dioxin) Receptor. Toxicol Sci. 2011, 124: 1-22. 10.1093/toxsci/kfr218.
Tang J, Vosman B, Voorrips RE, van der Linden CG, Leunissen JAM: QualitySNP: a pipeline for detecting single nucleotide polymorphisms and insertions/deletions in EST data from diploid and polyploid species. BMC Bioinforma. 2006, 7: 438-10.1186/1471-2105-7-438.
Wang J, Chuang K, Ahluwalia M, Patel S, Umblas N, Mirel D, Higuchi R, Germer S: High-throughput SNP genotyping by single-tube PCR with Tm-shift primers. Biotechniques. 2005, 39: 885-892. 10.2144/000112028.
Nacci D, Coiro L, Champlin D, Jayaraman S, McKinney R, Gleason TR, Munns WR, Specker JL, Cooper KR: Adaptations of wild populations of the estuarine fish Fundulus heteroclitus to persistent environmental contaminants. Mar Biol. 1999, 134: 9-17. 10.1007/s002270050520.
Nacci DE, Champlin D, Coiro L, McKinney R, Jayaraman S: Predicting the occurrence of genetic adaptation to dioxin-like compounds in populations of the estuarine fish Fundulus heteroclitus. Environ Toxicol Chem. 2002, 21: 1525-1532. 10.1897/1551-5028(2002)021<1525:PTOOGA>2.0.CO;2.
Excoffier L, Lischer HEL: Arlequin suite ver 3.5: a new series of programs to perform population genetics analysis under Linux and Windows. Mol Ecol Resour. 2010, 10: 564-567. 10.1111/j.1755-0998.2010.02847.x.
Wright S: Isolation by distance. Genetics. 1943, 28: 114-138.
Duvernell DD, Lindmeier JB, Faust KE, Whitehead A: Relative influences of historical and contemporary forces shaping the distribution of genetic variation in the Atlantic killifish, Fundulus heteroclitus. Mol Ecol. 2008, 17: 1344-1360. 10.1111/j.1365-294X.2007.03648.x.
Jensen JL, Bohonak AJ, Kelley ST: Isolation by distance, web service. BMC Genet. 2005, 6: 13-
Beaumont MA, Nichols RA: Evaluating loci for use in the genetic analysis of population structure. Proc R Soc Lond B Biol Sci. 1996, 263: 1619-1626. 10.1098/rspb.1996.0237.
Karchner SI, Franks DG, Powell WH, Hahn ME: Regulatory interactions among three members of the vertebrate aryl hydrocarbon receptor family: AHR repressor, AHR1 and AHR2. J Biol Chem. 2002, 277: 6949-6959. 10.1074/jbc.M110779200.
Carlson EA, Roy NK, Wirgin II: Microarray analysis of polychlorinated biphenol mixture-induced changes in gene epression among Atlantic tom cod populations displaying differential sensitivity to halogenated aromatic hydrocarbons. Environ Toxicol. 2009, 28: 759-771. 10.1897/08-195R.1.
Sartor MA, Schnekenburger M, Marlowe JL, Reichard JF, Wang Y, Fan Y, Ma C, Karyala S, Halbleib D, Liu X, Medvedovic M, Puga A: Genomewide analysis of aryl hydrocarbon receptor binding targets reveals an extensive array of gene clusters that control morphogenic and developmental programs. Environ Health Perspect. 2009, 117: 1139-1146. 10.1289/ehp.0800485.
Meyer JN, Volz DC, Freedman JF, DiGiulio RT: Differential display of hepatic mRNA from Fundulus heteroclitus inhabiting a superfund estuary. Aquat Toxicol. 2005, 73: 327-341. 10.1016/j.aquatox.2005.03.022.
Whitehead A, Triant DA, Champlin D, Nacci D: Comparative transcriptomics implicates mechanisms of evolved pollution tolerance in a killifish population. Mol Ecol. 2010, 19: 5186-5203. 10.1111/j.1365-294X.2010.04829.x.
Whitehead A, Pilcher W, Champlin D, Nacci D: Common mechanism underlies repeated evolution of extreme pollution tolerance. Proc R Soc Lond B Biol Sci. 2012, 279: 427-433. 10.1098/rspb.2011.0847.
Carney SA, Chen J, Burns CJ, Xiong KM, Peterson RE, Heideman W: Aryl hydrocarbon receptor activation produces heart-specific transcriptional and toxic responses in developing zebrafish. Mol Pharmacol. 2006, 70: 549-561. 10.1124/mol.106.025304.
Greytak SR, Calllard GV: Cloning of three estrogen receptors (ER) from killifish (Fundulus heteroclitus): differences in populations from polluted and reference environments. Gen Comp Endocrinol. 2007, 150: 174-188. 10.1016/j.ygcen.2006.07.017.
Fisher MA, Oleksiak MF: Convergence and divergence in gene expression among natural populations exposed to pollution. BMC Genomics. 2007, 8: 108-10.1186/1471-2164-8-108.
Paetzold SC, Ross NW, Richards RC, Jones M, Hellou J, Bard SM: Up-regulation of hepatic ABCC2, ABCG2, CYP1A1, and GST in multixenobiotic-resistant killifish (fundulus heteroclitus) from the Sydney Tar ponds, nova Scotia, Canada. Mar Environ Res. 2009, 68: 37-47. 10.1016/j.marenvres.2009.04.002.
Waits ER, Nebert DW: Genetic architecture of susceptibility to PCB-126 induced developmental cardiotoxicity in zebrafish. Toxicol Sci. 2011, 122: 466-475. 10.1093/toxsci/kfr136.
Wright S: Evolution and the genetics of populations. Variability within and among natural populations. Vol. 4. 1978, Chicago: University of Chicago Press
Hancock AM, Alkorta-Aranburu G, Witonsky DB, Di Rienzo A: Adaptations to new environments in humans: the role of subtle allele frequency shifts. Philos Trans R Soc Lond B Biol Sci. 2010, 365: 2459-2468. 10.1098/rstb.2010.0032.
Roark SA, Nacci D, Coiro L, Champlin D, Guttman SI: Population genetic structure of a nonmigratory estuarine fish (Fundulus heteroclitus) across a strong gradient of polychlorinated biphenyl contamination. Environ Toxicol and Chem. 2005, 24: 717-725. 10.1897/03-687.1.
Adams SM, Lindmeier JB, Duvernell DD: Microsatellite analysis of the phylogeography, Pleistocene history and secondary contact hypotheses for the killifish, fundulus heteroclitus. Mol Ecol. 2006, 15: 1109-1123. 10.1111/j.1365-294X.2006.02859.x.
McMillan AM, Bagley MJ, Jackson SA, Nacci DE: Genetic diversity and structure of an estuarine fish (Fundulus heteroclitus) indigenous to sites associated with a highly contaminated urban harbor. Ecotoxicology. 2006, 15: 539-548. 10.1007/s10646-006-0090-4.
Williams LM, Ma X, Boyko AR, Bustamante CD, Oleksiak MF: SNP identification, verification and utility for population genetics in a non-model genus. BMC Genet. 2010, 11: 32-
Whitehead A, Crawford DL: Neutral and adaptive variation in gene expression. Proc Natl Acad Sci U S A. 2006, 103: 5425-5430. 10.1073/pnas.0507648103.
Oleksiak MF, Karchner SI, Jenny MJ, Franks DG, Welch DBM, Hahn ME: Transcriptomic assessment of resistance to effects of an aryl hydrocarbon receptor (AHR) agonist in embryos of Atlantic killifish (Fundulus heteroclitus) from a marine superfund site. BMC Genomics. 2011, 12: 263-10.1186/1471-2164-12-263.
Hoffmann AA, Willi Y: Detecting genetic responses to environmental change. Nat Rev Genet. 2008, 9: 421-432.
Gahan L, Gould F, Heckel DG: Identification of a gene associated with Bt resistance in Heliothis virescens. Science. 2001, 293: 857-860. 10.1126/science.1060949.
Martinez DE, Levinton J: Adaptation to heavy metals in the aquatic oligochaete Limnodrilus hoffmeisteri: evidence for control by one gene. Evolution. 1996, 50: 1339-1343. 10.2307/2410674.
Meyer JN, Di Giulio RT: Heritable adaptation and associated fitness costs in killifish (Fundulus heteroclitus) inhabiting a contaminated estuary. Ecol Appl. 2003, 13: 490-503. 10.1890/1051-0761(2003)013[0490:HAAFCI]2.0.CO;2.
Hahn ME, Karchner SI, Evans BR, Franks DG, Merson RR, Lapseritis JM: Unexpected diversity of aryl hydrocarbon receptors in non-mammalian vertebrates: Insights from comparative genomics. J Exp Zool A Ecol Genet Physiol. 2006, 305A: 693-706.
Bak S-M, Iilda M, Hirano M, Iwato H, Kim E-Y: Potencies of Red seabream AHR1- and AHR2-mediated transactivation by dioxins: implicationsof both AHRs in dioxin toxicity. Environ Sci Technol. 2013, 47: 2877-2885. 10.1021/es304423w.
Mandel PK: Dioxin: a review of its environmental effects and its aryl hydrocarbon receptor biology. J Comp Physiol B Biochem Syst EnvironPhysiol. 2005, 175: 221-230. 10.1007/s00360-005-0483-3.
Williams LM, Oleksiak MF: Evolutionary and functional analyses of cytochrome P4501A promoter polymorphisms in natural populations. Mol Ecol. 2011, 20: 5236-5247. 10.1111/j.1365-294X.2011.05360.x.
Morrison HG, Weil EJ, Karchner SI, Sogin ML, Stegeman JJ: Molecular cloning of CYP1A from the estuarine fish Fundulus heteroclitus and phylogenetic analysis of CYP1 genes: update with new sequences. Comp Biochem Physiol Part C Toxicol Pharmacol. 1998, 121: 231-240.
Billiard SM, Timme-Laragy AR, Wessenberg DM, Cockman C, Di Giulio RT: The role of the aryl hydrocarbon receptor pathway in mediating synergistic developmental toxicity of polycyclic aromatic hydrocarbons to zebrafish. Toxicol Sci. 2006, 92: 526-536. 10.1093/toxsci/kfl011.
Goldstone H, Stegeman J: A revised evolutionary history of the CYP1A subfamily: gene duplication, gene conversion and positive selection. J Mol Evol. 2006, 62: 708-717. 10.1007/s00239-005-0134-z.
Zhou H, Qu Y, Wu H, Liao C, Zheng J, Diao X, Xue Q: Molecular phylogenies and evolutionary behavior of AhR (aryl hydrocarbon receptor) pathway genes in aquatic animals: implications for the toxicology mechanism of some persistent pollutants (POPs). Chemosphere. 2010, 87: 193-205.
Lutgens SPM, Cleutjens KBJM, Daemen MJAP, Heeneman S: Cathepsin cysteine proteases in cardiovascular disease. FASEB J. 2007, 21: 3029-3041. 10.1096/fj.06-7924com.
Reiser J, Adair B, Reinheckel T: Specialized roles for cysteine cathepsins in health and disease. J Clin Investig. 2010, 120: 3421-3431. 10.1172/JCI42918.
Hegelund T, Celandar MC: Hepatic versus extrahepatic expression of CYP3A30 and CYP3A56 in adult killifish (Fundulus heteroclitus). Aquat Toxicol. 2003, 64: 277-291. 10.1016/S0166-445X(03)00057-2.
Kim JW, Lee Y, Kang HB, Choe YK, Chung TW, Chang SY, Kwang SL, Choe IS: Cloning of the human cDNA sequence encoding the NADH: ubiquinone oxidoreductase MLRQ subunit. Biochem Mol Biol Int. 1997, 43: 669-675.
We owe many thanks to M. Hahn, S. Karchner, and D. Franks for sharing unpublished sequence and SNP data and for useful discussions. Helpful discussions with Eric Waits, technical support from Jason Grear, John Martinson, Mike Charpentier, Ian Kirby, Ashley Bertrand, and James Heltshe, and the advice from reviewers of earlier drafts, Mark Hahn, Sibel Karchner, Jeffrey Markert, John Martinson, Anne Kuhn, and 3 anonymous reviewers are also much appreciated. This is ORD tracking number ORD-005088 of the U.S. Environmental Protection Agency, Office of Research and Development, National Health and Environmental Effects Research Laboratory, Atlantic Ecology Division. This manuscript has been reviewed and approved for publication by the U.S. EPA. Approval does not signify that the contents necessarily reflect the views and policies of the U.S EPA. Mention of trade names, products, or services does not convey, and should not be interpreted as conveying official U.S. EPA approval, endorsement, or recommendation.
The authors declare they have no competing interests.
DAP, PF, DC, and DN contributed to study design, sample collection and processing. DP performed data analysis. DP and DN wrote the manuscript. All authors read and approved the final manuscript.
Patrick Flight, Denise Champlin contributed equally to this work.
Electronic supplementary material
Additional file 1: List of all genes considered for SNP analysis and their fate during the marker development process. Sequences provided by Sibel Karchner were derived from either 26 (AHR1, AHR2, and AHRR) or 10 individuals (not cDNA libraries) collected from Scorton Creek, MA. In some cases, Qualitiy SNP identified multiple available contigs from a single unigene entry. The number of SNPs genotyped does not equal the number of SNPs in the analysis because 8 of the markers failed to provide quality genotype information. N = No, Y = Yes, N/A = Not available, U = Unknown. (XLSX 23 KB)
About this article
Cite this article
Proestou, D.A., Flight, P., Champlin, D. et al. Targeted approach to identify genetic loci associated with evolved dioxin tolerance in Atlantic Killifish (Fundulus heteroclitus). BMC Evol Biol 14, 7 (2014). https://doi.org/10.1186/1471-2148-14-7