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.