The landscape of extreme genomic variation in the highly adaptable Atlantic killifish

Impacts on Atlantic Killifish from Neurotoxicants: Genes, Behavior, and Population-Relevant Outcomes

Molecular, cellular, and organismal alterations are important descriptors of toxic effects, but our ability to extrapolate and predict ecological risks is limited by the availability of studies that link measurable end points to adverse population relevant outcomes such as cohort survival and growth. In this study, we used laboratory gene expression and behavior data from two populations of Atlantic killifish Fundulus heteroclitus [one reference site (SCOKF) and one PCB-contaminated site (NBHKF)] to inform individual-based models simulating cohort growth and survival from embryonic exposures to environmentally relevant concentrations of neurotoxicants. Methylmercury exposed SCOKF exhibited brain gene expression changes in the si:ch211-186j3.6, si:dkey-21c1.4, scamp1, and klhl6 genes, which coincided with changes in feeding and swimming behaviors, but our models simulated no growth or survival effects of exposures. PCB126-exposed SCOKF had lower physical activity levels coinciding with a general upregulation in nucleic and cellular brain gene sets (BGS) and downregulation in signaling, nucleic, and cellular BGS. The NBHKF, known to be tolerant to PCBs, had altered swimming behaviors that coincided with 98% fewer altered BGS. Our models simulated PCB126 decreased growth in SCOKF and survival in SCOKF and NBHKF. Overall, our study provides a unique demonstration linking molecular and behavioral data to develop quantitative, testable predictions of ecological risk.

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

The genetic architecture of phenotypic traits can affect the mode and tempo of trait evolution. Human-altered environments can impose strong natural selection, where successful evolutionary adaptation requires swift and large phenotypic shifts. In these scenarios, theory predicts that adaptation is due to a few adaptive variants of large effect, but empirical studies that have revealed the genetic architecture of rapidly evolved phenotypes are rare, especially for populations inhabiting polluted environments. Fundulus killifish have repeatedly evolved adaptive resistance to extreme pollution in urban estuaries. Prior studies, including genome scans for signatures of natural selection, have revealed some of the genes and pathways important for evolved pollution resistance, and provide context for the genotype-phenotype association studies reported here. We created multiple quantitative trait locus (QTL) mapping families using progenitors from four different resistant populations, and using RAD-seq genetically mapped variation in sensitivity (developmental perturbations) following embryonic exposure to a model toxicant PCB-126. We found that one to two large-effect QTL loci accounted for resistance to PCB-mediated developmental toxicity. QTLs harbored candidate genes that govern the regulation of aryl hydrocarbon receptor (AHR) signaling. One QTL locus was shared across all populations and another was shared across three populations. One QTL locus showed strong signatures of recent natural selection in the corresponding wild population but another QTL locus did not. Some candidate genes for PCB resistance inferred from genome scans in wild populations were identified as QTL, but some key candidate genes were not. We conclude that rapidly evolved resistance to the developmental defects normally caused by PCB-126 is governed by few genes of large effect. However, other aspects of resistance beyond developmental phenotypes may be governed by additional loci, such that comprehensive resistance to PCB-126, and to the mixtures of chemicals that distinguish urban estuaries more broadly, may be more genetically complex.

The Ah Receptor from Toxicity to Therapeutics: Report from the 5th AHR Meeting at Penn State University, USA, June 2022

The aryl hydrocarbon receptor (AHR) is a sensor of low-molecular-weight molecule signals that originate from environmental exposures, the microbiome, and host metabolism. Building upon initial studies examining anthropogenic chemical exposures, the list of AHR ligands of microbial, diet, and host metabolism origin continues to grow and has provided important clues as to the function of this enigmatic receptor. The AHR has now been shown to be directly involved in numerous biochemical pathways that influence host homeostasis, chronic disease development, and responses to toxic insults. As this field of study has continued to grow, it has become apparent that the AHR is an important novel target for cancer, metabolic diseases, skin conditions, and autoimmune disease. This meeting attempted to cover the scope of basic and applied research being performed to address possible applications of our basic knowledge of this receptor on therapeutic outcomes.

Genome resequencing clarifies phylogeny and reveals patterns of selection in the toxicogenomics model Pimephales promelas

Background

The fathead minnow (Pimephales promelas) is a model species for toxicological research. A high-quality genome reference sequence is available, and genomic methods are increasingly used in toxicological studies of the species. However, phylogenetic relationships within the genus remain incompletely known and little population-genomic data are available for fathead minnow despite the potential effects of genetic background on toxicological responses. On the other hand, a wealth of extant samples is stored in museum collections that in principle allow fine-scale analysis of contemporary and historical genetic variation.

Methods

Here we use short-read shotgun resequencing to investigate sequence variation among and within Pimephales species. At the genus level, our objectives were to resolve phylogenetic relationships and identify genes with signatures of positive diversifying selection. At the species level, our objective was to evaluate the utility of archived-sample resequencing for detecting selective sweeps within fathead minnow, applied to a population introduced to the San Juan River of the southwestern United States sometime prior to 1950.

Results

We recovered well-supported but discordant phylogenetic topologies for nuclear and mitochondrial sequences that we hypothesize arose from mitochondrial transfer among species. The nuclear tree supported bluntnose minnow (P. notatus) as sister to fathead minnow, with the slim minnow (P. tenellus) and bullhead minnow (P. vigilax) more closely related to each other. Using multiple methods, we identified 11 genes that have diversified under positive selection within the genus. Within the San Juan River population, we identified selective-sweep regions overlapping several sets of related genes, including both genes that encode the giant sarcomere protein titin and the two genes encoding the MTORC1 complex, a key metabolic regulator. We also observed elevated polymorphism and reduced differentation among populations (F_ST) in genomic regions containing certain immune-gene clusters, similar to what has been reported in other taxa. Collectively, our data clarify evolutionary relationships and selective pressures within the genus and establish museum archives as a fruitful resource for characterizing genomic variation. We anticipate that large-scale resequencing will enable the detection of genetic variants associated with environmental toxicants such as heavy metals, high salinity, estrogens, and agrichemicals, which could be exploited as efficient biomarkers of exposure in natural populations.

Cardiac physiology and metabolic gene expression during late organogenesis among F. heteroclitus embryo families from crosses between pollution-sensitive and -resistant parents

BACKGROUND

The teleost fish Fundulus heteroclitus inhabit estuaries heavily polluted with persistent and bioaccumulative chemicals. While embryos of parents from polluted sites are remarkably resistant to toxic sediment and develop normally, embryos of parents from relatively clean estuaries, when treated with polluted sediment extracts, are developmentally delayed, displaying deformities characteristic of pollution-induced embryotoxicity. To gain insight into parental effects on sensitive and resistant phenotypes during late organogenesis, we established sensitive, resistant, and crossed embryo families using five female and five male parents from relatively clean and predominantly PAH-polluted estuaries each, measured heart rates, and quantified individual embryo expression of 179 metabolic genes.

RESULTS

Pollution-induced embryotoxicity manifested as morphological deformities, significant developmental delays, and altered cardiac physiology was evident among sensitive embryos resulting from crosses between females and males from relatively clean estuaries. Significantly different heart rates among several geographically unrelated populations of sensitive, resistant, and crossed embryo families during late organogenesis and pre-hatching suggest site-specific adaptive cardiac physiology phenotypes relative to pollution exposure. Metabolic gene expression patterns (32 genes, 17.9%, at p < 0.05; 11 genes, 6.1%, at p < 0.01) among the embryo families indicate maternal pollutant deposition in the eggs and parental effects on gene expression and metabolic alterations.

CONCLUSION

Heart rate differences among sensitive, resistant, and crossed embryos is a reliable phenotype for further explorations of adaptive mechanisms. While metabolic gene expression patterns among embryo families are suggestive of parental effects on several differentially expressed genes, a definitive adaptive signature and metabolic cost of resistant phenotypes is unclear and shows unexpected sensitive-resistant crossed embryo expression profiles. Our study highlights physiological and metabolic gene expression differences during a critical embryonic stage among pollution sensitive, resistant, and crossed embryo families, which may contribute to underlying resistance mechanisms observed in natural F. heteroclitus populations living in heavily contaminated estuaries.

Draft genome assemblies using sequencing reads from Oxford Nanopore Technology and Illumina platforms for four species of North American Fundulus killifish

BACKGROUND

Whole-genome sequencing data from wild-caught individuals of closely related North American killifish species (Fundulus xenicus, Fundulus catenatus, Fundulus nottii, and Fundulus olivaceus) were obtained using long-read Oxford Nanopore Technology (ONT) PromethION and short-read Illumina platforms.

FINDINGS

Draft de novo reference genome assemblies were generated using a combination of long and short sequencing reads. For each species, the PromethION platform was used to generate 30-45× sequence coverage, and the Illumina platform was used to generate 50-160× sequence coverage. Illumina-only assemblies were fragmented with high numbers of contigs, while ONT-only assemblies were error prone with low BUSCO scores. The highest N50 values, ranging from 0.4 to 2.7 Mb, were from assemblies generated using a combination of short- and long-read data. BUSCO scores were consistently >90% complete using the Eukaryota database.

CONCLUSIONS

High-quality genomes can be obtained from a combination of using short-read Illumina data to polish assemblies generated with long-read ONT data. Draft assemblies and raw sequencing data are available for public use. We encourage use and reuse of these data for assembly benchmarking and other analyses.

Association between Chromosome 4 and mercury accumulation in muscle of the three-spined stickleback (Gasterosteus aculeatus)

Anthropogenic stressors, such as pollutants, act as selective factors that can leave measurable changes in allele frequencies in the genome. Metals are of particular concern among pollutants, because of interference with vital biological pathways. We use the three-spined stickleback as a model for adaptation to mercury pollution in natural populations. We collected sticklebacks from 21 locations in Flanders (Belgium), measured the accumulated levels of mercury in the skeletal muscle tissue, and genotyped the fish by sequencing (GBS). The spread of muscle mercury content across locations was considerable, ranging from 21.5 to 327 ng/g dry weight (DW). We then conducted a genome-wide association study (GWAS) between 28,450 single nucleotide polymorphisms (SNPs) and the accumulated levels of mercury, using different approaches. Based on a linear mixed model analysis, the GWAS yielded multiple hits with a single top hit on Chromosome 4, with eight more SNPs suggestive of association. A second approach, a latent factor mixed model analysis, highlighted one single SNP on Chromosome 11. Finally, an outlier test identified one additional SNP on Chromosome 4 that appeared under selection. Out of all ten SNPs we identified as associated with mercury in muscle, three SNPs all located on Chromosome 4 and positioned within a 2.5 kb distance of an annotated gene. Based on these results and the genome coverage of our SNPs, we conclude that the selective effect of mercury pollution in Flanders causes a significant association with at least one locus on Chromosome 4 in three-spined stickleback.

Mummichog gill and operculum exhibit functionally consistent claudin-10 paralog profiles and Claudin-10c hypersaline response

Claudin (Cldn)-10 tight junction (TJ) proteins are hypothesized to form the paracellular Na⁺ secretion pathway of hyposmoregulating mummichog (Fundulus heteroclitus) branchial epithelia. Organ-specific expression profiles showed that only branchial organs [the gill and opercular epithelium (OE)] exhibited abundant cldn-10 paralog transcripts, which typically increased following seawater (SW) to hypersaline (2SW) challenge. Post-translational properties, protein abundance, and ionocyte localization of Cldn-10c, were then examined in gill and OE. Western blot analysis revealed two Cldn-10c immunoreactive bands in the mummichog gill and OE at ∼29 kDa and ∼40 kDa. The heavier protein could be eliminated by glycosidase treatment, demonstrating the novel presence of a glycosylated Cldn-10c. Protein abundance of Cldn-10c increased in gill and OE of 2SW-exposed fish. Cldn-10c localized to the sides of gill and OE ionocyte apical crypts and partially colocalized with cystic fibrosis transmembrane conductance regulator and F-actin, consistent with TJ complex localization. Cldn-10c immunofluorescent intensity increased but localization was unaltered by 2SW conditions. In support of our hypothesis, cldn-10/Cldn-10 TJ protein dynamics in gill and OE of mummichogs and TJ localization are functionally consistent with the creation and maintenance of salinity-responsive, cation-selective pores that facilitate Na⁺ secretion in hyperosmotic environments.

Summary: The role of claudin-10 tight junction proteins in paracellular salt secretion across fish branchial epithelia is indicated by organ-specific responses to hyperosmotic conditions and their association with salt secreting transcellular proteins

Chromosome-Level Assembly of the Atlantic Silverside Genome Reveals Extreme Levels of Sequence Diversity and Structural Genetic Variation

The levels and distribution of standing genetic variation in a genome can provide a wealth of insights about the adaptive potential, demographic history, and genome structure of a population or species. As structural variants are increasingly associated with traits important for adaptation and speciation, investigating both sequence and structural variation is essential for wholly tapping this potential. Using a combination of shotgun sequencing, 10x Genomics linked reads and proximity-ligation data (Chicago and Hi-C), we produced and annotated a chromosome-level genome assembly for the Atlantic silverside (Menidia menidia)-an established ecological model for studying the phenotypic effects of natural and artificial selection-and examined patterns of genomic variation across two individuals sampled from different populations with divergent local adaptations. Levels of diversity varied substantially across each chromosome, consistently being highly elevated near the ends (presumably near telomeric regions) and dipping to near zero around putative centromeres. Overall, our estimate of the genome-wide average heterozygosity in the Atlantic silverside is among the highest reported for a fish, or any vertebrate (1.32-1.76% depending on inference method and sample). Furthermore, we also found extreme levels of structural variation, affecting ∼23% of the total genome sequence, including multiple large inversions (> 1 Mb and up to 12.6 Mb) associated with previously identified haploblocks showing strong differentiation between locally adapted populations. These extreme levels of standing genetic variation are likely associated with large effective population sizes and may help explain the remarkable adaptive divergence among populations of the Atlantic silverside.

Extensive genome-wide duplications in the eastern oyster (Crassostrea virginica)

Genomic structural variation is an important source of genetic and phenotypic diversity, playing a critical role in evolution. The recent availability of a high-quality reference genome for the eastern oyster, Crassostrea virginica, and whole-genome sequence data of samples from across the species range in the USA, provides an opportunity to explore structural variation across the genome of this species. Our analysis shows significantly greater individual-level duplications of regions across the genome than that of most model vertebrate species. Duplications are widespread across all ten chromosomes with variation in frequency per chromosome. The eastern oyster shows a large interindividual variation in duplications as well as particular chromosomal regions with a higher density of duplications. A high percentage of duplications seen in C. virginica lie completely within genes and exons, suggesting the potential for impacts on gene function. These results support the hypothesis that structural changes may play a significant role in standing genetic variation in C. virginica, and potentially have a role in their adaptive and evolutionary success. Altogether, these results suggest that copy number variation plays an important role in the genomic variation of C. virginica. This article is part of the Theo Murphy meeting issue 'Molluscan genomics: broad insights and future directions for a neglected phylum'.

20th Pollutant Responses in Marine Organisms (PRIMO 20): Global issues and fundamental mechanisms caused by pollutant stress in marine and freshwater organisms

The 20th Pollutant Responses in Marine Organisms (PRIMO 20) conference provided a forum for scientists from around the world to communicate novel toxicological research findings specifically focused on aquatic organisms, by combining applied and basic research at the intersection of environmental and mechanistic toxicology. The work highlighted in this special issue of Aquatic Toxicology, a special issue of Marine Environmental Research, and presented through posters and presentations, encompass important and emerging topics in freshwater and marine toxicology. This includes multiple types of emerging contaminants including microplastics and UV filtering chemicals. Other studies aimed to further our understanding of the effects of endocrine disrupting chemicals, pharmaceuticals, and personal care products. Further research presented in this virtual issue examined the interactive effects of chemicals and pathogens, while the final set of manuscripts demonstrates continuing efforts to combine traditional biomonitoring, data from -omic technologies, and modeling for use in risk assessment and management. An additional goal of PRIMO meetings is to address the link between environmental and human health. Several articles in this issue of Aquatic Toxicology describe the appropriateness of using aquatic organisms as models for human health, while the keynote speakers, as described in the editorial below, presented research that highlighted bioaccumulation of contaminants such as PFOS and mercury from fish to marine mammals and coastal human populations such as the Gullah/GeeChee near Charleston, South Carolina, USA.

Sex-Specific Differences of Steroid Receptors Following Exposure to Environmentally Relevant Concentrations of Phenothiazine in Fundulus heteroclitus

Phenothiazine (PTZ) is a heterocyclic thiazine compound used for industrial and medical purposes. Through environmental surveillance studies, PTZ was found being discharged into a local river in Connecticut. Phenothiazine has been shown to act similarly to endocrine disrupting chemicals. This study sought to identify sex specific hormone receptor changes in Fundulus heteroclitus in response to PTZ exposure. Fundulus heteroclitus, also known as mummichog, are small fish native to the Atlantic coast of the United States and Canada. They reside in brackish waters and can survive harsh toxic environments. This model organism is native to the polluted waters found in Connecticut. In this study, fish were exposed to PTZ concentrations of 0.5 ppm, 1.0 ppm, and 2.0 ppm for 1 week. Following exposure, brain, liver, and gonad tissues were harvested; cDNA was synthesized; and mRNA expression was assessed for 6 different hormone receptors. Compared with vehicle control (ethanol) differences in mRNA expression, levels of hormone receptors were observed in various tissues from male and female fish. Many of the tissues assessed showed changes in expression level, while only female liver and testis showed no change. These results implicate PTZ as a potential endocrine disrupting compound to mummichog at environmentally relevant concentrations.

Using asexual vertebrates to study genome evolution and animal physiology: Banded (Fundulus diaphanus) x Common Killifish (F. heteroclitus) hybrid lineages as a model system

Wild, asexual, vertebrate hybrids have many characteristics that make them good model systems for studying how genomes evolve and epigenetic modifications influence animal physiology. In particular, the formation of asexual hybrid lineages is a form of reproductive incompatibility, but we know little about the genetic and genomic mechanisms by which this mode of reproductive isolation proceeds in animals. Asexual lineages also provide researchers with the ability to produce genetically identical individuals, enabling the study of autonomous epigenetic modifications without the confounds of genetic variation. Here, we briefly review the cellular and molecular mechanisms leading to asexual reproduction in vertebrates and the known genetic and epigenetic consequences of the loss of sex. We then specifically discuss what is known about asexual lineages of Fundulus diaphanus x F. heteroclitus to highlight gaps in our knowledge of the biology of these clones. Our preliminary studies of F. diaphanus and F. heteroclitus karyotypes from Porter's Lake (Nova Scotia, Canada) agree with data from other populations, suggesting a conserved interspecific chromosomal arrangement. In addition, genetic analyses suggest that: (a) the same major clonal lineage (Clone A) of F. diaphanus x F. heteroclitus has remained dominant over the past decade, (b) some minor clones have also persisted, (c) new clones may have recently formed, and iv) wild clones still mainly descend from F. diaphanus ♀ x F. heteroclitus ♂ crosses (96% in 2017-2018). These data suggest that clone formation may be a relatively rare, but continuous process, and there are persistent environmental or genetic factors causing a bias in cross direction. We end by describing our current research on the genomic causes and consequences of a transition to asexuality and the potential physiological consequences of epigenetic variation.

Tributyltin disrupts fin development in Fundulus heteroclitus from both PCB-sensitive and resistant populations: Investigations of potential interactions between AHR and PPARγ

Tributyltin (TBT) and dioxin-like polychlorinated biphenyls (PCBs) are environmental contaminants that are highly toxic to fish and co-occur in New Bedford Harbor (NBH), an estuarine Superfund site located in Massachusetts, USA. Atlantic killifish (Fundulus heteroclitus) that reside in NBH (and other highly contaminated sites along the east coast of the United States) have developed resistance to activation of the aryl hydrocarbon receptor (AHR) pathway and the toxicity of dioxin-like chemicals, such as 3,3',4,4',5-pentachlorobiphenyl, PCB126. In many biological systems, TBT disregulates adipose and bone development via the PPARγ-RXR pathway; AHR activation also disrupts adipose and bone homeostasis, potentially through molecular crosstalk between AHR and PPARγ. However, little is known about how co-exposure and the interaction of these pathways modulate the toxicological effects of these contaminants. Here, we tested the hypotheses that TBT would induce teratogenesis in killifish via activation of PPARγ and that PCB126 co-exposure would suppress PPARγ pathway activation in PCB-sensitive killifish from a reference site (Scorton Creek, SC, PCB-sensitive) but not in PCB-tolerant NBH killifish. Killifish embryos from both populations exposed to TBT (50 and 100 nM) displayed caudal fin deformities. TBT did not change the expression of pparg or its target genes related to adipogenesis (fabp11a and fabp1b) in either population. However, expression of osx/sp7, an osteoblast marker gene, and col2a1b, a chondroblast marker gene, was significantly suppressed by TBT only in SC killifish. An RXR-specific agonist, but not a PPARγ-specific agonist, induced caudal fin deformities like those observed in TBT-treated embryos. PCB126 did not induce caudal fin deformities and did not exacerbate TBT-induced fin deformities. Further, PCB126 increased expression of pparg in SC embryos and not NBH embryos, but did not change the expression of fabp1b. Taken together, these results suggest that in killifish embryos the PPARγ pathway is regulated in part by AHR, but is minimally active at least in this early life stage. In killifish, RXR activation, rather than PPARγ activation, appears to be the mechanism by which TBT induces caudal fin teratogenicity, which is not modulated by AHR responsiveness.

Mitochondria, sex and variation in routine metabolic rate

Variation in the metabolic costs associated with organismal maintenance may play a key role in determining fitness, and thus these differences among individuals are likely to be subject to natural selection. Although the evolvability of maintenance metabolism depends on its underlying genetic architecture, relatively little is known about the nature of genetic variation that underlies this trait. To address this, we measured variation in routine metabolic rate (ṀO_{2 routine} ), an index of maintenance metabolism, within and among three populations of Atlantic killifish, Fundulus heteroclitus, including a population from a region of genetic admixture between two subspecies. Polygenic association tests among individuals from the admixed population identified 54 single nucleotide polymorphisms (SNPs) that were associated with ṀO_{2 routine} , and these SNPs accounted for 43% of interindividual variation in this trait. However, genetic associations with ṀO_{2 routine} involved different SNPs if females and males were analysed separately, and there was a sex-dependent effect of mitochondrial genotype on variation in routine metabolism. These results imply that there are sex-specific genetic mechanisms, and potential mitonuclear interactions, that underlie variation in ṀO_{2 routine} . Additionally, there was evidence for epistatic interactions between 17% of the possible pairs of trait-associated SNPs, suggesting that epistatic effects on ṀO_{2 routine} are common. These data demonstrate not only that phenotypic variation in this ecologically important trait has a polygenic basis with considerable epistasis among loci, but also that these underlying genetic mechanisms, and particularly the role of mitochondrial genotype, may be sex-specific.

Developing a High-Quality Linkage Map for the Atlantic Killifish Fundulus heteroclitus

Killifish (Fundulus heteroclitus) are widely distributed among different aquatic environments where they demonstrate an impressive range of highly-plastic and locally adaptive phenotypes. High-throughput sequencing has begun to unravel the mechanisms and evolutionary history of these interesting features by establishing relationships in the genotype-phenotype map. However, some genotype-phenotype analyses require a higher order of contiguity than what initial scaffolded (fragmented genome assembly where contigs have been assemble into scaffolds) genome assemblies can provide. Here, we used 5,685 high-quality RAD-Seq markers from a single mapping family to order 84% of the scaffolded genome assembly to 24 chromosomes. This serves to: 1) expand the killifish genomic toolkit, 2) estimate genome-wide recombination rates, and 3) compare genome synteny to humans and other fishes. After initially building our map, we found that the selection of thresholds for sequence data filtration highly impacted scaffold placement in the map. We outline each step of the approach that dramatically improved our map to help guide others toward more effective linkage mapping for genome assembly. Our final map supports strong conservation of genomic synteny among closely related fish species and reveals previously described chromosomal rearrangements between more distantly related clades. However, we also commonly found minor scaffold misorientations in F. heteroclitus and in other assemblies, suggesting that further mapping (such as optical mapping) is necessary for finer scale resolution of genome structure. Lastly, we discuss the problems that would be expected from misoriented/unplaced scaffolds and stress the importance of a quality mapped genome as a key feature for further investigating population and comparative genomic questions with F. heteroclitus and other taxa.

Genome-Wide Stress Responses to Copper and Arsenic in a Field Population of Daphnia

Over the past decade, significant advances have been made to unravel molecular mechanisms of stress response in different ecotoxicological model species. Within this study, we focus on population level transcriptomic responses of a natural population of Daphnia magna Straus, (1820), to heavy metals. We aim to characterize the population level transcriptomic responses, which include standing genetic variation, and improve our understanding on how populations respond to environmental stress at a molecular level. We studied population level responses to two heavy metals, copper and arsenic, and their binary mixture across time. Transcriptomic patterns identified significantly regulated gene families and genes at the population level including cuticle proteins and resilins. Furthermore, some of these differentially regulated gene families, such as cuticle proteins, were also significantly enriched for genetic variations including SNPs and MNPs. In general, genetic variation was observed in specific gene families, many of which are known to be involved in stress response. Overall, our results indicate that molecular stress responses can be identified within natural populations and that linking molecular mechanisms with genetic variation at the population level could contribute significantly to adverse outcome frameworks.

Evolutionary Toxicogenomics of the Striped Killifish (Fundulus majalis) in the New Bedford Harbor (Massachusetts, USA)

In this paper, we used a Genotyping-by-Sequencing (GBS) approach to find and genotype more than 4000 genome-wide SNPs (Single Nucleotide Polymorphisms) from striped killifish exposed to a variety of polychlorinated biphenyls (PCBs) and other aromatic pollutants in New Bedford Harbor (NBH, Massachusetts, USA). The aims of this study were to identify the genetic consequences of exposure to aquatic pollutants and detect genes that may be under selection. Low genetic diversity (H_E and π) was found in the site exposed to the highest pollution level, but the pattern of genetic diversity did not match the pollution levels. Extensive connectivity was detected among sampling sites, which suggests that balanced gene flow may explain the lack of genetic variation in response to pollution levels. Tests for selection identified 539 candidate outliers, but many of the candidate outliers were not shared among tests. Differences among test results likely reflect different test assumptions and the complex pollutant mixture. Potentially, selectively important loci are associated with 151 SNPs, and enrichment analysis suggests a likely involvement of these genes with pollutants that occur in NBH. This result suggests that selective processes at genes targeted by pollutants may be occurring, even at a small geographical scale, and may allow the local striped killifish to resist the high pollution levels.

Profiling microRNA expression in Atlantic killifish (Fundulus heteroclitus) gill and responses to arsenic and hyperosmotic stress

The Atlantic killifish (Fundulus heteroclitus), native to estuarine areas of the Atlantic coast of the United States, has become a valuable ecotoxicological model as a result of its ability to acclimate to rapid environmental changes and adapt to polluted habitats. MicroRNAs (miRNAs) are highly conserved small RNAs that regulate gene expression and play critical roles in stress responses in a variety of organisms. Global miRNA expression in killifish and the potential roles miRNA have in environmental acclimation have yet to be characterized. Accordingly, we profiled miRNA expression in killifish gill for the first time and identified a small group of highly expressed, well-conserved miRNAs as well as 16 novel miRNAs not yet identified in other organisms. Killifish respond to large fluctuations in salinity with rapid changes in gene expression and protein trafficking to maintain osmotic balance, followed by a secondary phase of gene and protein expression changes that enable remodeling of the gills. Arsenic, a major environmental toxicant, was previously shown to inhibit gene expression responses in killifish gill, as well the ability of killifish to acclimate to a rapid increase in salinity. Thus, we examined the individual and combined effects of salinity and arsenic on miRNA expression in killifish gill. Using small RNA sequencing, we identified 270 miRNAs expressed in killifish, and found that miR-135b was differentially expressed in response to arsenic and at 24 h following transfer to salt water. Predicted targets of miR-135b are involved in ion transport, cell motility and migration, GTPase mediated signal transduction and organelle assembly. Consistent with previous studies of these two environmental stressors, we found a significant interaction (i.e., arsenic dependent salinity effect), whereby killifish exposed to arsenic exhibited an opposite response in miR-135b expression at 24 h post hyperosmotic challenge compared to controls. By examining mRNA expression of predicted miRNA targets during salinity acclimation and arsenic exposure, we found that miR-135b targets were significantly more likely to decrease during salinity acclimation than non-targets. Our identification of a significant interaction effect of arsenic and salinity on miR-135b expression supports the hypothesis that arsenic alters upstream regulators of stress response networks, which may adversely affect the killifish response to osmotic stress. The characterization of miRNAs in this ecotoxicological model will be a valuable resource for future studies investigating the role of miRNAs in response to environmental stress.

Hepatic metabolite profiling of polychlorinated biphenyl (PCB)-resistant and sensitive populations of Atlantic killifish (Fundulus heteroclitus)

Atlantic killifish inhabiting polluted sites along the east coast of the U.S. have evolved resistance to toxic effects of contaminants. One such contaminated site is the Acushnet River estuary, near New Bedford Harbor (NBH), Massachusetts, which is characterized by very high PCB concentrations in the sediments and in the tissues of resident killifish. Though killifish at this site appear to be thriving, the metabolic costs of survival in a highly contaminated environment are not well understood. In this study we compared the hepatic metabolite profiles of resistant (NBH) and sensitive populations (Scorton Creek (SC), Sandwich, MA) using a targeted metabolomics approach in which polar metabolites were extracted from adult fish livers and quantified. Our results revealed differences in the levels of several metabolites between fish from the two sites. The majority of these metabolites are associated with one-carbon metabolism, an important pathway that supports multiple physiological processes including DNA and protein methylation, nucleic acid biosynthesis and amino acid metabolism. We measured the gene expression of DNA methylation (DNA methyltransferase 1, dnmt1) and demethylation genes (Ten-Eleven Translocation (TET) genes) in the two populations, and observed lower levels of dnmt1 and higher levels of TET gene expression in the NBH livers, suggesting possible differences in DNA methylation profiles. Consistent with this, the two populations differed significantly in the levels of 5-methylcytosine and 5-hydroxymethylcytosine nucleotides. Overall, our results suggest that the unique hepatic metabolite signatures observed in NBH and SC reflect the adaptive mechanisms for survival in their respective habitats.

Integrative Population and Physiological Genomics Reveals Mechanisms of Adaptation in Killifish

Adaptive divergence between marine and freshwater (FW) environments is important in generating phyletic diversity within fishes, but the genetic basis of this process remains poorly understood. Genome selection scans can identify adaptive loci, but incomplete knowledge of genotype-phenotype connections makes interpreting their significance difficult. In contrast, association mapping (genome-wide association mapping [GWAS], random forest [RF] analyses) links genotype to phenotype, but offer limited insight into the evolutionary forces shaping variation. Here, we combined GWAS, RF, and selection scans to identify loci important in adaptation to FW environments. We utilized FW-native and brackish water (BW)-native populations of Atlantic killifish (Fundulus heteroclitus) as well as a naturally admixed population between the two. We measured morphology and multiple physiological traits that differ between populations and may contribute to osmotic adaptation (salinity tolerance, hypoxia tolerance, metabolic rate, body shape) and used a reduced representation approach for genome-wide genotyping. Our results show patterns of population divergence in physiological capabilities that are consistent with local adaptation. Population genomic scans between BW-native and FW-native populations identified genomic regions evolving by natural selection, whereas association mapping revealed loci that contribute to variation for each trait. There was substantial overlap in the genomic regions putatively under selection and loci associated with phenotypic traits, particularly for salinity tolerance, suggesting that these regions and genes are important for adaptive divergence between BW and FW environments. Together, these data provide insight into the mechanisms that enable diversification of fishes across osmotic boundaries.

Arsenic Reduces Gene Expression Response to Changing Salinity in Killifish

Toxicogenomic approaches can detect and classify adverse interactions between environmental toxicants and other environmental stressors but require more complex experimental designs and analytical approaches. Here we use novel toxicogenomic techniques to analyze the effect of arsenic exposure in wild killifish populations acclimating to changing salinity. Fish from three populations were acclimated to full strength seawater and transferred to fresh water for 1 or 24 h. Linear models of gene expression in gill tissue identified 31 genes that responded to osmotic shock at 1 h and 178 genes that responded at 24 h. Arsenic exposure (100 μg/L) diminished the responses (reaction norms) of these genes by 22% at 1 h ( p = 1.0 × 10^-6) and by 10% at 24 h ( p = 3.0 × 10^-10). Arsenic also significantly reduced gene coregulation in gene regulatory networks ( p = 0.002, paired Levene's test), and interactions between arsenic and salinity acclimation were uniformly antagonistic at the biological pathway level ( p < 0.05, binomial test). Arsenic's systematic interference with gene expression reaction norms was validated in a mouse multistressor experiment, demonstrating the ability of these toxicogenomic approaches to identify biologically relevant adverse interactions between environmental toxicants and other environmental stressors.

Genome-wide scan reveals signatures of selection related to pollution adaptation in non-model estuarine Atlantic killifish (Fundulus heteroclitus)

In many human-altered ecosystems, organisms are increasingly faced with more diverse and complex environmental stressors and pollutant mixtures, to which the adaptations necessary to survive exposure are likely to be numerous and varied. Improving our understanding of the molecular mechanisms that underlie complex polygenic adaptations in natural settings requires significant toxicological, biochemical, physiological, and genomic data rarely available for non-model organisms. Here, we build upon two decades of study of adaptation to anthropogenic pollutants in a population of Atlantic killifish (Fundulus heteroclitus) that inhabits the creosote-contaminated Atlantic Wood Industries Superfund (AW) site on the Elizabeth River, Virginia in the United States. To better understand the genotypes that underlie previously characterized resistance to PCBs and PAHs, we performed Restriction site-Associated DNA sequencing (RADseq) on killifish from AW and two relatively clean reference sites (King's Creek-KC, and Mains Creek-MC). Across the genome, we analyzed over 83,000 loci and 12,000 single nucleotide polymorphisms (SNPs). Shared across both comparisons of killifish from polluted (AW) and relatively unpolluted (KC and MC) sites, we found eight genomic regions with smoothed FST values significantly (p < 0.001) elevated above background. Using the recently published F. heteroclitus reference genome, we identified candidate genes in these significant regions involved in the AHR pathway (e.g. AIP, ARNT1c), as well as genes relating to cardiac structure and function. These genes represent both previously characterized and potentially novel molecular adaptations involved with various aspects of resistance to these environmental toxins.

The Toxicogenome of Hyalella azteca: A Model for Sediment Ecotoxicology and Evolutionary Toxicology

Hyalella azteca is a cryptic species complex of epibenthic amphipods of interest to ecotoxicology and evolutionary biology. It is the primary crustacean used in North America for sediment toxicity testing and an emerging model for molecular ecotoxicology. To provide molecular resources for sediment quality assessments and evolutionary studies, we sequenced, assembled, and annotated the genome of the H. azteca U.S. Lab Strain. The genome quality and completeness is comparable with other ecotoxicological model species. Through targeted investigation and use of gene expression data sets of H. azteca exposed to pesticides, metals, and other emerging contaminants, we annotated and characterized the major gene families involved in sequestration, detoxification, oxidative stress, and toxicant response. Our results revealed gene loss related to light sensing, but a large expansion in chemoreceptors, likely underlying sensory shifts necessary in their low light habitats. Gene family expansions were also noted for cytochrome P450 genes, cuticle proteins, ion transporters, and include recent gene duplications in the metal sequestration protein, metallothionein. Mapping of differentially expressed transcripts to the genome significantly increased the ability to functionally annotate toxicant responsive genes. The H. azteca genome will greatly facilitate development of genomic tools for environmental assessments and promote an understanding of how evolution shapes toxicological pathways with implications for environmental and human health.

The Toxicogenome of Hyalella azteca: A Model for Sediment Ecotoxicology and Evolutionary Toxicology

Hyalella azteca is a cryptic species complex of epibenthic amphipods of interest to ecotoxicology and evolutionary biology. It is the primary crustacean used in North America for sediment toxicity testing and an emerging model for molecular ecotoxicology. To provide molecular resources for sediment quality assessments and evolutionary studies, we sequenced, assembled, and annotated the genome of the H. azteca U.S. Lab Strain. The genome quality and completeness is comparable with other ecotoxicological model species. Through targeted investigation and use of gene expression data sets of H. azteca exposed to pesticides, metals, and other emerging contaminants, we annotated and characterized the major gene families involved in sequestration, detoxification, oxidative stress, and toxicant response. Our results revealed gene loss related to light sensing, but a large expansion in chemoreceptors, likely underlying sensory shifts necessary in their low light habitats. Gene family expansions were also noted for cytochrome P450 genes, cuticle proteins, ion transporters, and include recent gene duplications in the metal sequestration protein, metallothionein. Mapping of differentially expressed transcripts to the genome significantly increased the ability to functionally annotate toxicant responsive genes. The H. azteca genome will greatly facilitate development of genomic tools for environmental assessments and promote an understanding of how evolution shapes toxicological pathways with implications for environmental and human health.

Turning natural adaptations to oncogenic factors into an ally in the war against cancer

Both field and experimental evolution studies have demonstrated that organisms naturally or artificially exposed to environmental oncogenic factors can, sometimes rapidly, evolve specific adaptations to cope with pollutants and their adverse effects on fitness. Although numerous pollutants are mutagenic and carcinogenic, little attention has been given to exploring the extent to which adaptations displayed by organisms living in oncogenic environments could inspire novel cancer treatments, through mimicking the processes allowing these organisms to prevent or limit malignant progression. Building on a substantial knowledge base from the literature, we here present and discuss this progressive and promising research direction, advocating closer collaboration between the fields of medicine, ecology, and evolution in the war against cancer.

Ryanodine receptor and FK506 binding protein 1 in the Atlantic killifish (Fundulus heteroclitus): A phylogenetic and population-based comparison

Non-dioxin-like polychlorinated biphenyls (NDL PCBs) activate ryanodine receptors (RyR), microsomal Ca²⁺ channels of broad significance. Teleost fish may be important models for NDL PCB neurotoxicity, and we used sequencing databases to characterize teleost RyR and FK506 binding protein 12 or 12.6kDa (genes FKBP1A; FKBP1B), which promote NDL PCB-triggered Ca²⁺ dysregulation. Particular focus was placed on describing genes in the Atlantic killifish (Fundulus heteroclitus) genome and searching available RNA-sequencing datasets for single nucleotide variants (SNV) between PCB tolerant killifish from New Bedford Harbor (NBH) versus sensitive killifish from Scorton Creek (SC), MA. Consistent with the teleost whole genome duplication (tWGD), killifish have six RyR genes, corresponding to a and b paralogs of mammalian RyR1, 2 and 3. The presence of six RyR genes was consistent in all teleosts investigated including zebrafish. Killifish have four FKBP1; one FKBP1b and three FKBP1a named FKBP1aa, FKBP1ab, likely from the tWGD and a single gene duplicate FKBP1a3 suggested to have arisen in Atherinomorphae. The RyR and FKBP1 genes displayed tissue and developmental stage-specific mRNA expression, and the previously uncharacterized RyR3, herein named RyR3b, and all FKBP1 genes were prominent in brain. We identified a SNV in RyR3b encoding missense mutation E1458D. In NBH killifish, 57% were heterozygous and 28% were homozygous for this SNV, whereas almost all SC killifish (94%) lacked the variant (n≥39 per population). The outlined sequence differences between mammalian and teleost RyR and FKBP1 together with outlined population differences in SNV frequency may contribute to our understanding of NDL PCB neurotoxicity.

Annotation of the Nuclear Receptors in an Estuarine Fish species, Fundulus heteroclitus

The nuclear receptors (NRs) are ligand-dependent transcription factors that respond to various internal as well as external cues such as nutrients, pheromones, and steroid hormones that play crucial roles in regulation and maintenance of homeostasis and orchestrating the physiological and stress responses of an organism. We annotated the Fundulus heteroclitus (mummichog; Atlantic killifish) nuclear receptors. Mummichog are a non-migratory, estuarine fish with a limited home range often used in environmental research as a field model for studying ecological and evolutionary responses to variable environmental conditions such as salinity, oxygen, temperature, pH, and toxic compounds because of their hardiness. F. heteroclitus have at least 74 NRs spanning all seven gene subfamilies. F. heteroclitus is unique in that no RXRα member was found within the genome. Interestingly, some of the NRs are highly conserved between species, while others show a higher degree of divergence such as PXR, SF1, and ARα. Fundulus like other fish species show expansion of the RAR (NR1B), Rev-erb (NR1D), ROR (NR1F), COUPTF (NR2F), ERR (NR3B), RXR (NR2B), and to a lesser extent the NGF (NR4A), and NR3C steroid receptors (GR/AR). Of particular interest is the co-expansion of opposing NRs, Reverb-ROR, and RAR/RXR-COUPTF.