Independently evolved pollution resistance in four killifish populations is largely explained by few variants of large effect

2.1. Breeding design

We created mapping families from the same source populations as those described in Reid et al. (2016), which have previously been characterized for sensitivities to PCB‐126 (Nacci et al., 2010). DLC‐resistant source populations originated from three sites from the northern portion of the species' range: New Bedford Harbor, MA, USA (NBH), Bridgeport Harbor, CT, USA (BRP), Newark, NJ, USA (NEW), and one site from the southern portion of their range: the Atlantic Wood site of the Elizabeth River, Virginia, USA (ELR) (Figure 1). The DLC‐sensitive source population originated from Block Island, RI, USA (BLI). Sediment DLC and/or PAH concentrations are orders of magnitude higher in the four polluted sites compared to the clean site, and magnitude of resistance correlates with the magnitude of pollution (Nacci et al., 2010) (Table S1).

We created QTL mapping families using wild killifish that were maintained in the laboratory (Narragansett EPA), where individuals from BLI (sensitive population) were mated with individuals from each resistant population (NBH, BRP, NEW, or ELR) to produce F1 hybrid families. F1 hybrid individuals were subsequently paired to produce the F2 mapping families. The specific identities of mapping families were NBH 1105, BRP 812, NEW 1413, and ELR 1433. Over the course of this study, we tried slightly differing breeding strategies to maximize the number of progeny with divergent phenotypes (Figure S2; also, see Figure S3 for the distributions of phenotypes per family). The NBH and ELR F2 hybrid mapping families descended from mating one wild‐caught female killifish from the resistant population (NBH or ELR) and one wild‐caught male killifish from BLI. The NEW mapping family descended from mating one laboratory‐bred female killifish from BLI and one laboratory‐bred male killifish from NEW, while the BRP‐mapping family descended from two lines, each produced by mating one wild‐caught BLI female and one wild‐caught BRP male. Thus, resistance was inherited through females in the NBH and ELR families, and through males in the NEW and BRP families (Figure S2). Previous studies have demonstrated that inheritance of resistance is not sex linked (Nacci et al., 1999). Furthermore, there was a single tolerant progenitor for NBH, ELR, and NEW families (sib: sib matings), while the BRP family was derived from two tolerant progenitors (nonsib mating) (Figure S2). To characterize tolerance, we challenged embryos from each of the four F2 mapping families with a discriminating dose of the model DLC (PCB‐126), scored individuals for developmental abnormalities associated with DLC toxicity, and then archived them for later genotyping (described below).

2.2. Phenotyping

We exposed embryos from all families to acetone‐solubilized PCB‐126 dissolved in sea water at a concentration (200 ng/L; ~0.6 nM) that discriminates sensitive from resistant killifish. PCB‐126 was selected because it is a highly toxic PCB congener (Doering et al., 2023) that is abundant at polluted sites (Nacci et al., 2010), and exposures to PCB‐126 reveal large differences in sensitivity between populations from clean and polluted sites (e.g., Nacci et al., 2010). The exposure lasted from 1 day post‐fertilization (dpf) to 7 dpf, followed by transfer to clean seawater. At 10 dpf, we observed the embryos microscopically and scored them for the presence or absence of a suite of developmental deformities that are characteristic of DLC‐mediated toxicity in fish embryos, where the appearance of multiple abnormalities is tightly correlated with reduced probability of larval survival and is, therefore, a good proxy for DLC sensitivity. Each abnormality is given a score of 1 (present) or 0 (absent), and each embryo was scored as the sum of the observed abnormalities, from 0 (no deformities) to 5 (most deformities), which represents the full range of sensitivity displayed by originating killifish populations (representative images of normal and severely affected individuals, Figure S4). Some aspects of phenotype scoring are subjective as to the presence or severity; phenotyping was, therefore, conducted with multiple experienced observers who had each trained together to minimize observer bias. The abnormalities commonly observed and the method of scoring each were as follows. Hemorrhaging: Most commonly hemorrhaging was found in the head and/or tail, occasionally in other areas such as the yolk sac or pericardial space; each occurrence of a hemorrhage in a different body area was given a score of 1 (i.e., a hemorrhage in both the tail and the head would result in 2 points). Reduced size: At 10 dpf, a normally developed embryo fills the chorion such that the tail has curled around in front of the fish and overlaps the eye, possibly even extending further; body size was scored as reduced and given a score of 1 if the caudal fin did not reach to overlap the fish's eye. Heart abnormality: At 10 dpf, the heart of normally developing embryos sits tightly in the space between the underside of the embryo's jaw and the yolk sac; they exhibit a fully looped heart with an imaginary line between the centers of the ventricle and atrium sitting roughly parallel to a line between the centers of the embryo's eyes. A hallmark of fish embryonic DLC toxicity is the loss of this structure and occurrence of what is commonly known as a “stringy heart”. Embryos were scored for the presence of a heart abnormality if they exhibited a loss of obvious differentiation of the two heart chambers, stretching of the heart such that the chambers appear vertically stacked (in a line with the inflow and outflow tracts as opposed to side by side), or in severe cases, the full loss of heart structure such that all that remained was a stringy heart tube. Pericardial edema: As described previously, the heart of a normally developing 10 dpf embryo sits tightly between the underside of the fish's jaw and the yolk sac and the pericardial space is small (sometimes it can be observed bulging out by the sides of the embryo's head), with the beating of the heart typically depressing the yolk sac slightly. Pericardial edema has historically been considered a sensitive hallmark of fish embryonic DLC toxicity, so it contributes to the deformity score as either 1 point for minor presence or 2 for severe edema. Embryos scored as a 1 for edema exhibited a larger space around the heart such that the embryo's head is lifted off of the yolk sac due to the visible enlargement of the pericardial space and the beating of the heart no longer depresses the yolk sac. In an observation of severe edema, scored as a 2, the pericardial space is greatly enlarged and the embryo's head is lifted far from the yolk sac. At this level of severity, the embryo can be turned on its side and the observer can typically look fully through the pericardial space and out the other side of the embryo with no other structures observed due to the amount of fluid in the pericardium. Finally, the phenotyping also included space for an “Other” observation which also resulted in a score of 1 and allowed for rare occurrences to be included in scoring (e.g., severe spinal curvature, severe developmental delay, edema of the yolk sac). Although the number of abnormalities described above could result in an individual with a higher score, the abnormality score was capped at 5 (4 for BRP); in our experience, any individual scored 4 or higher does not survive past 2 weeks post‐hatch. After phenotyping, we flash froze embryos and stored them at −80°C until DNA extraction. Phenotype scores for all individuals selected for genotyping are in Appendix S1.

2.3. Genotyping

We used balanced selective genotyping (or two‐tail genotyping) (Lander & Botstein, 1989; Lebowitz et al., 1987), where we selected ~48 individuals from each of the most sensitive and most resistant ends of the phenotypic distribution for each family (we did not genotype all individuals from each family because of limited resources). We used RAD‐seq (Miller et al., 2007) for genotyping these individuals from F2 intercross families and mapping family founders. DNA was extracted from frozen embryos with a proteinase‐K digestion and the Qiagen DNAeasy kit. Sequencing libraries were prepared for RAD‐seq by ligating individual barcodes and NEBnext Illumina oligonucleotides to genomic DNA at SbfI cut sites. Four lanes of paired‐end (PE‐100) sequence data (Illumina HiSeq 2500) were collected (one lane per plate of 96 samples) to enable genotyping of offspring and founders at RAD sites. After de‐multiplexing by barcode and evaluating the quality of each sample with FASTQC, sequence reads were aligned to the linkage‐mapped F. heteroclitus assembly (Miller et al., 2019; EBI BioStudies accession S‐BSST163; this map orders scaffolds from the F. heteroclitus reference genome assembly Fundulus_heteroclitus‐3.0.2, NCBI BioProject PRJNA177717). We assigned read group information and aligned reads with BWA‐MEM using default parameters, marked duplicates with SAMBLASTER, sorted reads with SAMTOOLS, and flagged improperly paired reads with BAMTOOLS (Barnett et al., 2011; Faust & Hall, 2014). We used FreeBayes to call genotypes on all populations simultaneously (Garrison & Marth, 2012). We clustered sample genotypes with metric MDS for QAQC and to cluster genotypes at chromosome 5 to assign sex to the immature embryos (Figure S5) so that sex could be included as a covariate in the QTL search and model fit. Families were each split into separate datasets and filtered for genotypes from repetitive sequence alignments and family‐specific invariant sites with PLINK 1.9 (Chang et al., 2015).

Genotypes for each mapping family were formatted to be loaded as independent crosses in R/qtl (Broman et al., 2003). In R/qtl, offspring genotypes were filtered based on the number of genotypes per locus (>12.5% missing data) and the level of segregation distortion (p < 0.001) from a 1:2:1 ratio following R/qtl recommendations (Broman & Sen, 2009). Outcrosses in R/qtl are assumed to be from inbred founders, so we limited our analysis to markers that were homozygous for alternate alleles in the founders (double heterozygotes in F1 parents). Our segregation distortion filter was conservative because selective genotyping is expected to lead to some distortion, particularly at QTL of large effect (Xu, 2008). After filtration, markers that were genotyped in the founders and F2 embryos were assigned to a consistent allele (A or B) within and among linkage groups. Any markers that we could not confirm with founder genotypes were filtered out so that those remaining were clearly linked with other markers in the linkage group (formLinkageGroups in R/qtl with a recombination frequency <0.1 and LOD score of 10). We initially anchored markers to their physical position along chromosomes, but then used recombination frequency to reorder (assign a different order than the initial meiotic map) markers and estimated mapping distance in this order with the Kosambi mapping algorithm (Kosambi, 1943). Differences from the initial meiotic map may be the result of genome structural variation within a mapping population or read mapping and genome assembly errors. We then checked whether the final filtered marker set sufficiently covered the reference genome physical map and provided a sensible mapping order by visualizing the pairwise recombination frequency and LOD linkage matrix for each chromosome, as well as the relationship between genetic and physical distances along chromosomes (Figure S6). After filtering, the number of linked markers were: 37,304 for the NBH family; 16,725 for the BRP family; 27,119 for the NEW family; and 57,489 for the ELR family (Table S2).

In addition to genome‐wide genotyping with RAD‐seq, we also sought to confirm genotypes at two discrete loci: at the loci that encode aryl hydrocarbon receptor 2a/1a (AHR2a/1a) on chromosome 1 and aryl hydrocarbon receptor‐interacting protein (AIP) on chromosome 2. These loci encode proteins that are core components of the signaling pathway that is adaptively de‐sensitized in resistant fish (Whitehead et al., 2012), were implicated as a QTL in the NBH population (Nacci et al., 2016), and show among the strongest signatures of natural selection in resistant populations (Reid et al., 2016). Since the adaptive haplotypes are not fixed in resistant populations (Reid et al., 2016), we sought to confirm whether adaptive haplotypes were segregating in our mapping families. We genotyped ELR offspring with sensitive (scores 4 and 5) and resistant (scores 0 and 1) phenotypes for a deletion that spans the last exon of AHR2a and the first six exons of neighboring AHR1a (99.7 kb deletion; Figure S7) which exists at 81% frequency in wild ELR fish but is absent in sensitive fish from a nearby reference population (methods described in Reid et al., 2016). The genotypes at the AHR2a/1a locus (presence/absence of a deletion; Table S3) were added to the QTL analysis with RAD‐Tag markers to test for linkage with other chromosomal markers and for an association with the resistant phenotype. Two AIP nonsynonymous SNPs at amino acid positions 224 and 252 were genotyped in a subset of mapping family individuals (n = 8 for each family, 4 each sensitive/resistant, except for n = 35 in ELR for SNP252; Tables S4 and S5) and in the founders of each cross. We amplified a 1.4 kb genomic fragment with PCR primers AIP3F (5′‐GGCGCTATACCCGCTCGTGTCC‐3′) and AIP5R2 (5′‐CTTCATATTTGAAGACGAGGGAGG‐3′) using 10 ng genomic DNA and Advantage DNA polymerase (Clontech) with the following cycling conditions: [94°C, 1 min]; [94°C, 5 s; 68°C, 2 min] 35×; [68°C, 5 min]. The amplified product was direct‐sequenced with the AIP5R2 primer for SNP analysis.

2.4. QTL analysis

Prior to interval mapping, we performed a marker regression on un‐filtered markers to set the priors for exploring QTL model‐space for interval mapping (Figure S8). We then performed interval mapping on the filtered markers (see above) in R/qtl to test the single QTL model on each chromosome via multiple imputation (n = 500) and Haley–Knott regression, which estimates the most likely interval genotypes and QTL position between RAD markers from the log posterior distribution. If the LOD score for the single QTL test exceeded our permuted threshold (top 15% of the null distribution) in a single QTL scan, the position of the highest LOD score was added to the full QTL model. We compared multiple single and full‐QTL models between methods (normal, binary, transformed data with parametric and nonparametric models). We examined models that included only the initial two major QTL discovered by marker regression, as well as more inclusive models that included the minor effect QTL that were discovered by interval mapping. QTL were dropped if they did not significantly change the fit of the full‐model in the drop‐one‐qtl analysis. The full model in R/qtl estimated the additive and dominance effect of each of the QTL in the model. We compared their relative effect sizes to illuminate the genetic architecture of resistance in each of the mapping families.

2.5. Comparison of methods to previous QTL experiment

Our methods were similar but not identical to those of (Nacci et al., 2016). Here, we contrast the approaches of these two studies to facilitate interpretation of overlap between our findings. The 2016 study created multiple families with different progenitor individuals from one resistant population (NBH), whereas the current study created one family from each of four resistant populations. The 2016 study genotyped individuals with microsatellites and SNPs resulting in ~240 genetic markers (Waits et al., 2016), whereas the current study genotyped using between 16,000 and 57,000 SNP markers (RAD‐seq) depending on the family. Both studies used the same breeding design (F2 intercross), the same exposure and phenotype scoring regime, both used two‐tailed genotyping, and similar statistical models. Note that the genetic map upon which the 2016 QTL study was based (Waits et al., 2016) had some assembly mistakes (Appendix S2); most notably it was missing chromosome 18, thereby precluding the 2016 study from revealing the large‐effect QTL that we observe on Chr18 (see Section 3).

2.6. Comparison to selection scan data

To further refine our understanding of the genetic variation that contributes to resistance, we tested for overlap between previously identified signatures of selection from population genomics studies (Reid et al., 2016) and the QTL identified here by aligning both datasets to the same chromosomal coordinates. This approach has been effective for connecting QTLs to strong selection from domestication and breeding (Rubin et al., 2012). Signatures of selection were identified from whole‐genome scans of diversity and divergence between each of the four resistant populations and nearby nonresistant populations (Reid et al., 2016). In the population genomics studies, genomic scaffolds were scanned for extreme values of p _i, F _ST, and Tajimas D in 5 kb windows to identify signatures of strong recent natural selection. The selection signature thresholds for each statistic were set using a demographic model to simulate a null distribution of each statistic. These thresholds were used to convert the statistics into an aggregate z‐score. The Z‐scores (which included the width and magnitude of the signature region) formed a composite rank for each signature of selection in each polluted population. F _ST, delta p _i, and composite ranks on the scaffold‐level genome assembly were lifted over to a physical position on the killifish chromosomal assembly (Miller et al., 2019, Fundulus_heteroclitus‐3.0.2, NCBI BioProject PRJNA177717) with liftover using a chain‐file generated by ALLMAPS (Hinrichs et al., 2006; Tang et al., 2015). QTL intervals were evaluated against the 50 and 200 top‐ranked composite signatures of selection from Reid et al. (2016).

3.1. Phenotypes

We observed the full spectrum of sensitive to resistant phenotypes in mapping families, indicating that the founders from each of the populations harbored genetic variation for resistance to PCB‐126. Frequency distributions of phenotypic variation for the families that were chosen for mapping (Figure 1b) were representative of phenotype distributions from other F2 intercross families created from these same populations (Figure S3). In our analyses, we grouped individuals with phenotypes of 4 and 5 as sensitive and individuals with phenotypes of 0 and 1 as resistant and modeled the resistance phenotype as binary.

3.2. QTL analysis

We identified a small number of loci that account for a large proportion of phenotypic variation in PCB sensitivity in QTL mapping families from all four resistant populations. Genome‐wide marker regression and QTL interval mapping of RAD‐Tag markers identified one or two large‐effect QTL in each of the families and these QTL tended to be shared between multiple mapping families (Figure 2a,b). A large‐effect QTL on CHR18 was shared among all four resistant populations. This CHR18 QTL accounted for a large fraction of resistance variation in each family (11% in NBH, 27% in BRP, 40% in NEW, and 34% in ELR). Another large‐effect QTL on CHR2 was shared among the three northern resistant populations. This CHR2 QTL accounted for a large fraction of resistance variation in each northern family (40% in NBH, 22% in BRP, and 18% in NEW). Thus, families from the three northern populations shared the same set of two major‐effect QTLs on chromosomes 2 and 18 (CHR2 and CHR18, respectively), despite a different breeding design and phenotypic distribution for the BRP family. Together, this set of two QTLs accounted for a very large fraction of variation in sensitivity to PCB‐induced toxicity in northern families (70% in NBH, 53% in BRP, and 55% in NEW). We also found that, at least for the NBH mapping family, the QTL on CHR2 appears to be recessive, and the few heterozygous individuals that are phenotypically resistant carry the resistant genotype on CHR18 (Figure 2c). Parameter estimates for the large‐effect QTL models are included in Table S6 and File S1. We detected support for medium to small effect QTL that do not overlap between families. When included in the full model, these additional QTL explain a relatively small proportion of variation (variance explained <6%, drop‐one‐qtl analysis estimates). We included parameter estimates for large effect loci with the more inclusive models as shown in Table S6 and File S1.

3.3. QTL and signatures of selection

We tested for an overlap between our QTL peaks and regions of the genome showing signatures of recent strong natural selection from genome‐wide scans (Reid et al., 2016), such that overlaps between QTL and selection signatures could narrow the interval for inferring candidate genes. We found that both major QTLs contain genomic regions encoding core components of the AHR signaling pathway (Appendix S3). The two large‐effect QTL intervals found on CHR2 and CHR18 (Figure 2) include AIP and AHR genes, respectively. As defined by high allele frequency differentiation between resistant and reference populations (F _ST) and reduced nucleotide diversity in the resistant population (delta‐p _i), the AIP genomic region of CHR2 ranks among the strongest signatures of natural selection in all resistant populations (Reid et al., 2016) (Figure 3). This region overlaps with the large‐effect QTL for the three northern‐population crosses (NBH, BRP, NEW; Figure 3), but not for the ELR mapping family (Figures 2 and 3). The large‐effect QTL for all four resistant populations on CHR18 (Figure 2) contains genes AHR1b and AHR2b. This region included the 20th ranked signature of selection in ELR, 164th in NEW, but did not include highly ranked signatures of selection for NBH or BRP (Figure 3).

In addition to the large‐effect QTL, we detected support for smaller effect QTL. Markers on chromosomes 8,11,13, and 24 exceeded the permuted threshold in the single QTL scans. However, we found nominal support for the effect and position of these QTL in the full‐model with drop‐one‐qtl analysis. Of the moderate/low effect QTL, the CHR8 QTL had the highest individual LOD score in NBH, but only explained <2.3% variation when included in the full model. None of these smaller‐effect QTL coincide with genomic regions showing strong signatures of selection in the wild‐resistant populations from which mapping families were derived.

We performed targeted genotyping of the mapping families to verify the presence or the absence of putatively adaptive haplotypes that were discovered by the whole‐genome scans for selection. In ELR, a highly ranking signature of recent natural selection on CHR1 was previously discovered to contain a 99.7‐kb deletion affecting coding sequences of AHR1a and AHR2a (Figure S7A), and this deletion had swept from low frequency in a nearby reference population (not detected) to high frequency (81%) in the resistant ELR population (Reid et al., 2016). PCR tests, using primers inside and outside of the deletion region (Figure S7B), verified that the AHR1a/AHR2a deletion on CHR1 is segregating within our ELR mapping family (Table S3). However, in ELR we detect no evidence for a QTL associated with variable sensitivity to PCB‐126 toxicity anywhere on CHR1 (Figure S8). The PCR genotypes of 40 resistant individuals from our ELR mapping family indicated that eight are homozygous for the deletion, 22 are heterozygous, and 10 are homozygous intact. Genotyping of 47 sensitive individuals from our mapping family indicated six homozygous for the deletion, 25 heterozygous, and 16 homozygous intact (Table S3). The phenotype association LOD score of the PCR deletion genotype and RAD markers near the deletion haplotype do not support a QTL on CHR1 for ELR or any of the mapping families (LOD of 0.051 for the deletion genotype in ELR).

Additional targeted genotyping of nonsynonymous SNPs in the coding sequence for AIP (CHR2) was used to confirm whether the adaptive genotype at that locus (inferred from results of selection scans from Reid et al., 2016) was segregating within our mapping families. We were particularly keen to determine why the ELR mapping family was the only one where the genomic region containing the AIP locus was not a QTL, even though the AIP locus is within one of the highest ranked signatures of natural selection in all populations including ELR (the AIP locus was within the third highest ranked signature of natural selection detected in ELR; Reid et al., 2016). We examined two nonsynonymous SNPs in AIP that have high allele frequency differences between resistant populations and their paired sensitive reference populations (Tables S3 and S4). These variants fall within the region containing the adaptive haplotype on CHR2 (Figure 3), and we hereby refer to these variants as “adaptive” insofar as they are indicative of large haplotype frequency differences (and they are in a candidate gene). The Trp → Leu variant in AIP amino acid position 224 (aa224) shows large frequency differences between northern resistant populations (where the Leu residue is between 70% and 74% frequency) and their paired reference populations (where the Leu residue is between 0% and 4% frequency) (Table S4). This suggests that the variant is linked to a variant(s) that promotes fitness advantage in polluted estuaries. The Northern resistant founders (NBH, BRP, NEW) of mapping families are all homozygous for the aa224 Leu variant, which is consistent with the support of a QTL at AIP in the Northern mapping families (NBH, BRP, NEW). However, the Trp residue at aa224 is fixed in both the ELR population and its paired reference population (Table S4). In contrast, a different amino acid variant (Ile → Val at aa252) shows large allele frequency differentiation at the AIP locus between ELR and its paired reference population, where the Val residue is absent in the reference population, but at 48% frequency in ELR (Table S5). We also found that both the Ile and Val residues were segregating within our ELR mapping family at equal frequency, but variation at this locus was not predictive of resistance (Table S5). This indicates that putatively adaptive variation at the AIP locus was captured in all of our mapping families but was only predictive of resistance to PCB‐126 in the three northern populations and not predictive of resistance in the ELR population.

4.1. Genetic complexity of evolved PCB resistance

Two major large‐effect QTL were identified among the four killifish families, with the full model that included these two QTL accounting for a very large proportion of phenotypic variation in developmental sensitivity to our model toxicant (34.2%–70.1%) (Figure 2). One major QTL was shared among the four mapping families: the QTL on CHR18, which includes candidate genes AHR1b/2b paralogs, accounted for 11%–43% of the variation in DLC sensitivity. The other major QTL was shared among the Northern resistant killifish but was not identified as a QTL in ELR: the QTL on CHR2, which includes candidate gene AIP, accounted for 18.1%–45.4% of the variance in DLC sensitivity. A couple of other loci, private to some of the families, accounted for a small fraction of phenotypic variation in some families. For example, minor QTL were noted on CHR8 in the NBH family and on CHR11 in the NEW family. These results are largely congruent with the prior NBH QTL study, where a small number of QTL accounted for 69% of phenotypic variation, with the largest effect loci associated with AHR pathway genes and a region on CHR8 (Nacci et al., 2016). We conclude that variants in a small number of key genes that regulate AHR signaling largely achieve the desensitization of AHR signaling that is protective of DLC‐induced development toxicity in urban killifish populations. These results indicate that the genetic basis of evolved variation in developmental sensitivity to PCB‐126 is not complex or highly polygenic. Resistance evolved rapidly through increased frequency of large‐effect loci, where the adaptive alleles are found at very low frequency (or are undetectable) in sensitive populations from clean habitats (Reid et al., 2016). This core finding is consistent with expectations when adaptation requires a rapid and large shift in the phenotypic optimum (Stetter et al., 2018).

It is possible that variants in AHR and/or AIP may result in similar phenotypic outcomes, by effecting adaptive desensitization of PCB‐induced AHR signaling. Variation at the DLC‐binding site of the AHR is highly predictive of DLC sensitivity in bird species (e.g., Farmahin et al., 2013). Similarly, evolved adaptive DLC resistance in Hudson River (NY, USA) tomcod is largely explained by variants in AHR localized near the site where AIP complexes with AHR (Wirgin et al., 2011). Moreover, AIP is an evolutionarily conserved AHR chaperone that influences the stability of the AHR complex and eventual nuclear translocation of AHR (Bell & Poland, 2000; Lecoq et al., 2016; Meyer et al., 2000). AIP loss‐of‐function mutations in mice protect against some forms of DLC toxicity (Nukaya et al., 2010). However, the role of AHR–AIP interactions in fish, and whether there are ligand‐specific differences in AHR–AIP interactions, are not well understood. Convergent evolution that is underpinned by convergent genetic changes, as we observe with these two AHR pathway genes, is more likely when selection is strong and genes have large phenotypic effects (MacPherson & Nuismer, 2017). Our results are consistent with these expectations, insofar as selection has been strong, and AHR and AIP mutations have large effect sizes. The likelihood of genetic convergence also increases when functional constraints on relevant pathways are high (Feldman et al., 2012; Weinreich et al., 2006). We hypothesize that changes to either AHR or AIP locus may result in an altered AHR‐AIP‐agonist complex that modifies the molecular events that would normally initiate AHR signal transduction. Whether this convergence emerges from pleiotropic constraints, bias in the likelihood of certain types of mutations (e.g., Martin & Orgogozo, 2013), shared selective environments, or shared genetic background, or some combination of these, remains to be experimentally determined.

Few studies have revealed the genomic basis of evolved pollution resistance (but see, e.g., Wright et al., 2013, 2015). However, many studies have examined adaptations to environmental pollutants that are toxic by design, such as pesticides, which may offer relevant lessons. Often few loci of major effect are implicated in evolved resistance to pesticides, for example, loci that are the biochemical targets of toxicity or loci that more broadly govern metabolic detoxification (Ffrench‐Constant, 2013). Although the literature on the genetic basis of pesticide resistance is relatively large, most studies focus on candidate genes; genome‐wide approaches are rare (e.g., Kreiner et al., 2019, 2021; Pélissié et al., 2022; Van Etten et al., 2020). Using knowledge from pesticide resistance to make predictions about the genomic basis of pollution resistance may be fraught, partly because genome‐wide understanding of the relevant genetics is thin, and partly because the nature of selection in polluted environments may be different from that in agricultural settings. For example, pesticides and antibiotics are typically designed to target very specific biochemical processes, and they are usually encountered in high concentrations as single chemicals or simple mixtures. In contrast, pollution can present fitness challenges that are more complex, because these chemicals may perturb multiple biochemical pathways (Gerstein et al., 2012, 2015), and multiple pollutants often co‐occur in the environment (Whitehead et al., 2017). Toxic chemicals are often encountered in human‐modified environments where multiple co‐occurring habitat alterations may contribute additional selective pressures (e.g., altered competitive interactions, thermal and flow regimes, hypoxia), which further increases the dimensionality of fitness challenges and the required adaptations (Iriart et al., 2021; Rivkin et al., 2019).

Many more loci than those few QTL reported here appear to be important for fitness in urban populations of killifish (Reid et al., 2016). It is, therefore, important to consider whether more broadly defined resistance to PCB toxicity, and fitness in polluted urban estuaries even more generally, has a simple or complex genetic basis. In this study, we only examined resistance to a limited suite of toxic effects caused by a single dose of a single PCB congener during one life stage. Different sets of mutations may underlie sensitivity to different doses (Wang & Kruglyak, 2014), or at different life stages (Everman et al., 2021). PCB‐126 is only one representative of the many DLCs and other pollutant classes that may affect fitness of urban fish. We currently do not know whether other aspects of DLC toxicity in embryos or later life stages (e.g., immune dysfunction, endocrine dysfunction, neurological dysfunction, or cancer [White & Birnbaum, 2009]) are resolved in resistant fish by the same or different QTL as those identified here. Similarly, our model toxicant PCB‐126 is only one of a complex mixture of chemicals that pollute these urban estuaries and threaten fitness. Populations inhabiting these sites have likely evolved resistance to many of those chemicals (e.g., Bugel et al., 2014; Celander et al., 2021; Grans et al., 2015; Greytak et al., 2010). Resistance to some of those chemicals may be underpinned by the same loci as we identified here, but other loci are also likely to contribute, especially for chemicals with unique mechanisms of toxicity. Furthermore, it is plausible that initial costs of adaptive AHR desensitization caused by mutations of large effect may have prompted compensatory adaptive fine‐tuning (MacLean et al., 2010) to restore some function to AHR signaling and connected pathways. By extension of the above reasoning, the many environmental changes that distinguish urban from natural estuaries include diverse nonchemical perturbations (e.g., biotic interactions, flow regimes, pathogens, etc.) that may drive adaptation through entirely unique suites of traits, pathways, and genes. We conclude that evolved resistance to the developmental toxicity induced by PCB‐126 has a relatively simple genetic basis in resistant individuals, while acknowledging that the genetic changes required for adaptation to urban estuaries are likely much more complex.

4.2. Variation in QTL among resistant populations

Our mapping families derived from four focal populations share some similarities and also some differences in the loci that contribute to resistance. Mapping families from the three northern populations share each of the two major‐effect QTL (on CHR2 and CH18), whereas the mapping family from the southern population (ELR) is most different and shares just one QTL (on CHR18) with those from the northern populations. A lack of parallelism at the genetic level following parallel resistance to toxicants has been observed in many systems including pesticide resistance in weeds (Van Etten et al., 2020) and beetles (Pélissié et al., 2022), and heavy metal resistance in sea campion (Papadopulos et al., 2021). It is often concluded that this is consistent with a complex genetic architecture for resistance, but this is rarely tested experimentally (e.g., with quantitative genetics). However, for copper resistance in fruit flies, an extensive overlap was observed between QTL and selection scan outliers, although QTL were clearly a small subset of the selection scan outliers (Everman et al., 2021, 2023). For killifish, where the genetic architecture is not complex, similarities and differences between QTL for different populations may be driven by any of three factors: (1) the presence of functionally redundant variants, (2) differences in genetic background, or (3) differences between selective environments; we consider these in sequence.

First, for many adaptations, the genetic basis is polygenic. This could manifest as many variants, each of small effect, contributing to phenotypic differences between two individuals (one resistant, one sensitive). However, this does not appear to be the case for the data presented here. Alternatively, polygenic adaptation could manifest from multiple redundant large‐effect variants contributing to adaptive fitness among different pairs of individuals (e.g., Goldstein & Holsinger, 1992). Indeed, there was limited overlap in PCB resistance QTL among the three mapping families created with different NBH progenitors (Nacci et al., 2016). This type of polygenic adaptation, perhaps better termed “oligogenic,” is likely relevant for toxicant‐resistant killifish populations, insofar as the resistance phenotype is fixed but almost no SNPs are fixed in those populations (Reid et al., 2016). Therefore, individuals selected as QTL mapping family progenitors may harbor different combinations of adaptive loci (e.g., as observed in Nacci et al., 2016). If adaptation is polygenic or oligogenic, then QTL studies (because they include only one or very few progenitors) may not reveal the complete sets of loci that contribute to evolved phenotypes. However, we do not consider this a good explanation for the QTL differences on CHR2 presented here. Variants on CHR2 have been subject to strong recent natural selection in all four populations (Figure 3; Osterberg et al., 2018; Reid et al., 2016), yet association with resistance was found only in NBH, BRP, and NEW families but not the ELR family. Perhaps the adaptive haplotype was missing from the ELR individual that we selected as the mapping family progenitor? We consider this unlikely since both putatively adaptive and sensitive genotypes at this locus were segregating in the ELR mapping family (Table S5). However, it is still possible that these variants are not closely linked with the adaptive variant(s). Ongoing high‐density genetic mapping should resolve this.

The second factor that could contribute to inter‐family variation in QTL is genetic background. Convergent phenotypic evolution is more likely to be underpinned by convergent genetic evolution when genetic variation is ancestrally shared and selection is strong (e.g., Alves et al., 2019; Jones et al., 2012; Therkildsen et al., 2019). The three northern populations are genetically quite similar to each other (average pairwise F _ST = 0.09), whereas the southern population is genetically distinct (average pairwise F _ST between ELR and each of the northern populations = 0.19) (Reid et al., 2016). Accordingly, the genetic variation available for contemporary natural selection was different between the northern populations and the southern population, such that the genotypes underlying parallel adaptation were shared in the north but different in the south. However, we do not consider this sufficient to explain the QTL differences presented here. Although the adaptive haplotypes were different between northern and southern populations, the genomic region surrounding the CHR2 QTL was among the very top‐ranked signatures of natural selection in all four resistant populations (Figure 3; Osterberg et al., 2018; Reid et al., 2016). If this adaptive variation supported resistance to PCBs, then the region on CHR2 should be a QTL in all four populations.

The third factor that could contribute to inter‐family variation in QTLs is differences between the environments that select for variants at different adaptive loci. The northern populations are all predominantly polluted by halogenated DLCs (PCBs and dioxins), whereas the southern ELR site is primarily polluted by PAHs. Although DLCs and PAHs may exert some aspects of developmental toxicity in part through a common mechanism (AHR signaling; Clark et al., 2010), the detailed mechanisms of action of these two chemical classes differ (Billiard et al., 2008; Denison et al., 2011; Incardona, 2017). One important difference is that PAHs are readily metabolized by AHR‐regulated proteins (potentially creating many reactive metabolites), whereas PCBs are not. Another difference is that dioxins and dioxin‐like PCBs mediate cardiac toxicity largely through the AHR, but many PAHs exert their cardiac toxicity through AHR and AHR‐independent mechanisms (Brette et al., 2014; Incardona et al., 2005). Clearly, variation at CHR2 has been subject to strong selection in all four populations (Figure 3). Perhaps PCB‐induced selection in northern estuaries favored variation at CHR2 that enables PCB‐resistant phenotypes, whereas PAH‐induced selection in ELR also favored variation at CHR2 but those variants promote aspects of adaptation (e.g., PAH resistance) that are independent of PCB resistance. We consider this a plausible explanation for our data. Ongoing experiments are testing for QTL for multiple classes of chemical resistance (including PAHs) in multiple killifish individuals and populations.

4.3. Coincidence between QTL and signatures of natural selection

Some QTLs for PCB resistance showed signatures of selection in wild populations, but some did not. The QTL on CHR2 in mapping families from NBH, BRP, and NEW overlaps with a top‐ranked signature of selection in genome scans of those same populations. In contrast, the QTL on CHR18 found in all four mapping families was a ranked signature of selection in only two of the four wild populations; this locus was ranked 20th and 164th in ELR and NEW populations, respectively, but no evidence for recent strong selection was detected in genome scans of BRP and NBH populations (Figure 3). One possible explanation is that multiple functionally redundant adaptive variants at this locus contributed to resistance in wild killifish populations, or a single adaptive variant embedded in multiple genetic backgrounds contributed to resistance. The soft sweeps or sets of partial sweeps that would emerge to underlie adaptation are notoriously hard to detect because of diminished signal of linked selection (Pritchard et al., 2010; Yeaman, 2015). Additional studies that, for example, take advantage of linkage‐disequilibrium based tests (Pennings & Hermisson, 2006) would be necessary to detect a soft sweep on CHR18 and test this hypothesis.

Many loci that showed strong signatures of selection in wild populations (Reid et al., 2016) were not QTL in this study. Genome scans may identify loci that affect fitness in wild populations, but do not directly implicate the agent of selection. Urban estuaries likely pose not only many fitness challenges for killifish, including chemical pollution, but also changes in biotic interactions and in nutrient, temperature, hypoxia, and flow regimes. The full suite of adaptive solutions to these complex and multifarious fitness challenges can imprint in genome scans. It is, therefore, not surprising that loci that are QTLs for PCB resistance are only a small subset of loci that are subject to selection in urban estuaries (e.g., as observed in evolved copper resistance in fruit flies [Everman et al., 2021, 2023]). However, we were surprised that a large deletion that spans AHR2a and AHR1a on CHR1 (Figure S7A) was not associated with PCB resistance in the ELR mapping family. This is surprising because AHR desensitization through AHR knockout/knockdown is protective of PCB (or dioxin) toxicity in rodents (Fernandez‐Salguero et al., 1996), zebrafish (Garcia et al., 2018; Goodale et al., 2012; Prasch et al., 2003; Souder & Gorelick, 2019), and killifish (Clark et al., 2010), and an AHR2a deletion had swept to high frequency in resistant populations of Atlantic tomcod (Wirgin et al., 2011), Gulf killifish (Oziolor et al., 2019), and the ELR population of Atlantic killifish (Reid et al., 2016). Since the deletion was segregating in the ELR mapping family (Table S3), we conclude that it must be associated with some other aspect of evolved pollution sensitivity.

One locus that showed among the strongest signatures of selection in all wild populations was also strongly associated with PCB resistance but not in all populations. The locus on CHR2 was the top ranked signature of selection in all four wild resistant populations. Adaptive variants at this locus were segregating in all four mapping families (Tables S3 and S4). However, this locus was a major‐effect QTL in only the three northern populations but not in the southern ELR family (Figure 2). We speculate that strong selection at AIP confers resistance to PCB toxicity in northern populations, but selection at this locus contributes to other aspects of adaptation in the ELR population. Since PAH pollution distinguishes the ELR site from the others, and AHR signaling is important for PAH toxicity, perhaps variation at AIP is important for evolved resistance to PAH toxicity that is independent of the mechanisms that drive resistance to PCB toxicity.

Conspicuously, two tandem AHR paralogs show strong signatures of selection in the wild, whereas two other AHR paralogs are QTL. The AHR paralogs AHR1a and AHR2a are found on CHR1. This locus is under strong selection in all four urban killifish populations (Proestou et al., 2014; Reid et al., 2016; Reitzel et al., 2014). However, variation on CHR1 was not associated with PCB resistance in any of our mapping families. This was unexpected because AHR2a was shown to mediate developmental toxicity of PCBs and PAHs in killifish (Clark et al., 2010). However, in contrast, the QTL region on CHR18 harbors the other two AHR paralogs AHR1b and AHR2b. This region is a QTL in all populations, but shows only a moderately ranked signature of selection in only two of the wild populations (Reid et al., 2016). In zebrafish, AHR2b mediates the developmental toxicity of PCB‐like chemicals (Garcia et al., 2018; Goodale et al., 2012; Prasch et al., 2003; Souder & Gorelick, 2019) and AHR1b protein binds with DLCs (Karchner et al., 2005). It is, therefore, plausible that the AHR1b/2b locus plays a role in PCB toxicity, but the relative roles of each of the four killifish AHR paralogs remain to be resolved. All four AHR proteins exhibit high‐affinity binding to [³H]TCDD in vitro and TCDD‐inducible transcriptional activation in cell culture (Karchner, Franks, Hahn, unpublished results), suggesting that DLC effects in killifish could involve multiple AHRs. Evidence for strong selection at AHR1a/2a but no QTL for PCB resistance at that locus, and QTL for PCB resistance at AHR1b/2b with evidence for weaker and perhaps polygenic selection at that locus, suggests that evolved pollution resistance is underlain by complex interactions among AHR paralogs and their protein partners (e.g., AIP). In any case, it is evident from both genotype–phenotype mapping and scans for signatures of recent natural selection that genetic variation in core regulatory components of the AHR signaling pathways contributes to evolved resistance to pollution in urban estuaries.

4.4. Caveats and considerations

Selective genotyping may have limited power to detect epistatic or small effect loci. We chose selective genotyping because of the power and efficiency (costs savings) that are afforded (Lander & Botstein, 1989). However, modeling studies suggest that if there are loci of large effect then the effectiveness of selective genotyping becomes unpredictable (Sen et al., 2009). Clearly, our results reveal loci of large effect, but small effect loci and nonadditive or epistatic loci may also be present (e.g., as in Nacci et al., 2016). Indeed, our data suggest that the QTL on CHR2 is recessive (Figure 2c). Ongoing studies by our group involve genotyping and phenotyping entire families that span the entire range of sensitive to resistant phenotypes, which should resolve this potential issue.

Although our results indicate that resistance to PCB‐126‐induced developmental toxicity is underlain by few variants of large effect, one should be cautious in concluding that this evolved resistance is not more polygenic. This is because our experiments included only one resistant individual as the progenitor for each mapping family (except for BRP which had two resistant progenitors), so sampling of within‐population variation in resistance was far from thorough. By sampling many more individuals, GWAS studies could reveal more loci that contribute to resistance, especially if resistance is underlain by multiple redundant variants that are not fixed within resistant populations. For example, GWAS for evolved herbicide resistance has revealed many loci, even when large‐effect loci were present, and these additional loci were segregating across a range of allele frequencies, suggesting a complex and polygenic architecture (Kreiner et al., 2021). For killifish, some QTL varied between three NBH mapping families in a previous study (Nacci et al., 2016). However, our three mapping families derived from the three resistant northern populations (NBH, BRP, and NEW) all implicated the same major‐effect QTL (Figure 2). It is plausible that minor‐effect and epistatic QTL, which we had little power to detect in our study, may be variable among individuals, so that the complexity of resistance to PCB‐126‐induced developmental toxicity during development may be greater than it currently appears.

One should also be cautious about concluding that all aspects of resistance to PCB‐126, or adaptation to pollution more broadly, or adaptation to urban estuaries even more broadly, is not highly polygenic. This is simply because multiple traits, each with a simple or complex genetic basis, likely contribute to fitness in these highly altered urban environments (as we discuss in more detail in the section above “Genetic complexity of evolved PCB resistance”). Additional quantitative genetic studies are needed to draw genotype–phenotype associations for the diversity of phenotypes induced by exposure to the array of stressors that distinguish urban estuaries.