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

Transcriptomic and Methylomic Analyses Show Significant Shifts in Biosynthetic Processes and Reduced Intrapopulation Gene Expression Variance in PAH-Adapted Atlantic Killifish

Environmental contaminants pose a significant selection pressure across taxa, potentiating evolved resistance to chemicals. However, rapid evolution may alter molecular and physiological homeostasis leading to trade-offs. To elucidate molecular underpinnings of evolved chemical resistance, we compared liver gene expression and methylation profiles in polycyclic aromatic hydrocarbon (PAH)-adapted Atlantic killifish (Fundulus heteroclitus) in the Republic site (RP), Elizabeth River, Virginia with PAH-sensitive Kings Creek (KC) fish. We found 1607 differentially expressed and 2252 alternatively spliced genes between RP and KC, with highly enriched genes involving lipid and amino acid metabolism, respectively. While 308 genes had differentially methylated regions, only 13 of these genes were differentially expressed. The aryl hydrocarbon receptor 2b gene (ahr2b) was differentially methylated and expressed, as well as alternatively spliced signifying its critical role in mediating PAH tolerance. Notably, the intrapopulation coefficient of variation (CoV) was lower in 82% of 17,566 expressed genes in RP fish compared to KC fish. Among other pathways, these genes with low CoV were highly enriched in bioenergetic processes inferring reduced metabolic physiological variation as a population in RP fish. Altered metabolic gene expression and overall reduced gene expression variance in RP fish warrant further studies on fitness trade-offs including altered susceptibility to other stressors associated with rapid adaptation to anthropogenic pressures.

Limited impact of weathered residues from the Deepwater Horizon oil spill on the gut-microbiome and foraging behavior of sheepshead minnows (Cyprinodon variegatus)

The Deepwater Horizon disaster of April 2010 was the largest oil spill in U.S. history and exerted catastrophic effects on several ecologically important fish species in the Gulf of Mexico (GoM). Within fish, the microbiome plays a key symbiotic role in maintaining host health and aids in acquiring nutrients, supporting immune function, and modulating behavior. The aim of this study was to examine if exposure to weathered oil might produce significant shifts in fish gut-associated microbial communities as determined from taxa and genes known for hydrocarbon degradation, and whether foraging behavior was affected. The gut microbiome (16S rRNA and shotgun metagenomics) of sheepshead minnow (Cyprinodon variegatus) was characterized after fish were exposed to oil in High Energy Water Accommodated Fractions (HEWAF; tPAH = 81.1 ± 12.4 µg/L) for 7 days. A foraging behavioral assay was used to determine feeding efficiency before and after oil exposure. The fish gut microbiome was not significantly altered in alpha or beta diversity. None of the most abundant taxa produced any significant shifts as a result of oil exposure, with only rare taxa showing significant shifts in abundance between treatments. However, several bioindicator taxa known for hydrocarbon degradation were detected in the oil treatment, primarily Sphingomonas and Acinetobacter. Notably, the genus Stenotrophomonas was detected in high abundance in 16S data, which previously was not described as a core member of fish gut microbiomes. Data also demonstrated that behavior was not significantly affected by oil exposure. Potential low bioavailability of the oil may have been a factor in our observation of minor shifts in taxa and no behavioral effects. This study lays a foundation for understanding the microbiome of captive sheepshead minnows and indicates the need for further research to elucidate the responses of the fish gut-microbiome under oil spill conditions.

Contamination acts as a genotype-dependent barrier to gene flow, causing genetic erosion and fine-grained population subdivision in Mussels from the Strait of Istanbul

This study provides evidence of fine-grained genetic structuring in Mediterranean mussels (Mytilus galloprovincialis) from the Strait of Istanbul, caused by barriers to gene flow via contaminant-mediated selection. In this study, mitochondrial D-loop sequences were analyzed in mussels from 8 localities, all less than 30 kilometers apart, with differing contaminant loads. The results were: 1) Intra-population genetic differentiation (ΦST) between sites with high and low contaminant loads was high (up to 0.459), even at distances of only a few kilometers. 2) Genetic diversity was negatively correlated with the contaminant load ("genetic erosion"). 3) There was evidence of selection, based on haplotype frequencies and neutrality tests (Tajima's D), with purifying selection at the most contaminated site and balancing selection at the least contaminated. 4) Genetic distance was not correlated with geographic distance (no isolation-by-distance), but was correlated with contaminant load at each site. 5) Population dendrograms and Bayesian estimators of migration indicated that gene flow between sites was affected by contamination. For the dendrograms of the sampling sites, the clades clustered according to contaminant load more than geographic distance. Overall, these results suggest that 1) contamination may serve as a genotype-dependent dispersal barrier (i.e., contamination may not affect total number of migrants, just the relative proportions of the haplotypes in the established immigrants), leading strong population differentiation over short distances, and 2) genetic erosion may occur by a combination of selection and altered patterns of haplotype-specific gene flow. These effects may be more pronounced in the Strait of Istanbul than in other locations because of the riverine nature and strong, uni-directional current of the strait.

Differential DNA methylation and metabolite profiling of Atlantic killifish (Fundulus heteroclitus) from the New Bedford Harbor Superfund site

Atlantic killifish (Fundulus heteroclitus) is a valuable model in evolutionary toxicology to study how the interactions between genetic and environmental factors serve the adaptive ability of organisms to resist chemical pollution. Killifish populations inhabiting environmental toxicant-contaminated New Bedford Harbor (NBH) show phenotypes tolerant to polychlorinated biphenyls (PCBs) and differences at the transcriptional and genomic levels. However, limited research has explored epigenetic alterations and metabolic effects in NBH killifish. To identify the involvement of epigenetic and metabolic regulation in the adaptive response of killifish, we investigated tissue- and sex-specific differences in global DNA methylation and metabolomic profiles of NBH killifish populations, compared to sensitive populations from a non-polluted site, Scorton Creek (SC). The results revealed that liver-specific global DNA hypomethylation and differential metabolites were evident in fish from NBH compared with those from SC. The sex-specific differences were not greater than the tissue-specific differences. We demonstrated liver-specific enriched metabolic pathways (e.g., amino acid metabolic pathways converged into the urea cycle and glutathione metabolism), suggesting possible crosstalk between differential metabolites and DNA hypomethylation in the livers of NBH killifish. Additional investigation of methylated gene regions is necessary to understand the functional role of DNA hypomethylation in the regulation of enzyme-encoding genes associated with metabolic processes and physiological changes in NBH populations.

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.

Bioenergetic Effects of Polycyclic Aromatic Hydrocarbon Resistance Manifest Later in Life in Offspring of Fundulus heteroclitus from the Elizabeth River

Shifts in key physiological processes can confer resistance to chemical pollutants. However, these adaptations may come with certain trade-offs, such as altered energy metabolic processes, as evident in Atlantic killifish (Fundulus heteroclitus) in Virginia's Elizabeth River (ER) that have evolved resistance to polycyclic aromatic hydrocarbons (PAHs). We seek to understand the bioenergetic costs of PAH resistance among subpopulations of Atlantic killifish with differing contamination levels in order to examine how these changes manifest across multiple life stages and how these costs might be exacerbated by additional stressors. Bioenergetics data revealed differences in metabolic rates between offspring of PAH-resistant fish and reference fish were absent or minimal in both the embryo and larval stages but pronounced at the juvenile life stage, suggesting that bioenergetic changes in pollution-adapted killifish manifest later in life. We also provide evidence that killifish from remediated sites are more sensitive to PAH exposure than killifish from nonremediated sites, suggesting loss of PAH tolerance following relaxed selection. Collectively, our data suggest that the fitness consequences associated with evolved resistance to anthropogenic stressors may manifest differently over time and depend on the magnitude of the selection pressure. This information can be valuable in effective risk and remediation assessments as well as in broadening our understanding of species responses to environmental change.

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 the influence of 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 genetically mapped variation in sensitivity (developmental perturbations) following embryonic exposure to a model toxicant PCB-126. We found that a few 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, where some (but not all) of these QTL loci were shared across all populations, and some (but not all) of these loci showed signatures of recent natural selection in the corresponding wild population. Some strong 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.

Identification of Candidate Genes and Functional Pathways Associated with Body Size Traits in Chinese Holstein Cattle Based on GWAS Analysis

Body size is one of the most economically important traits of dairy cattle, as it is significantly associated with cow longevity, production, health, fertility, and environmental adaptation. The identification and application of genetic variants using a novel genetic approach, such as genome-wide association studies (GWASs), may give more insights into the genetic architecture of complex traits. The identification of genes, single nucleotide polymorphisms (SNPs), and pathways associated with the body size traits may offer a contribution to genomic selection and long-term planning for selection in dairy cows. In this study, we performed GWAS analysis to identify the genetic markers and genes associated with four body size traits (body height, body depth, chest width, and angularity) in 1000 Chinese Holstein cows. We performed SNPs genotyping in 1000 individuals, based on the GeneSeek Genomic Profiler Bovine 100 K. In total, we identified 11 significant SNPs in association with body size traits at the threshold of Bonferroni correction (5.90 × 10−7) using the fixed and random model circulating probability unification (FarmCPU) model. Several genes within 200 kb distances (upstream or downstream) of the significant SNPs were identified as candidate genes, including MYH15, KHDRBS3, AIP, DCC, SQOR, and UBAP1L. Moreover, genes within 200 kb of the identified SNPs were significantly enriched (p ≤ 0.05) in 25 Gene Ontology terms and five Kyoto Encyclopedia of Genes and Genomes pathways. We anticipate that these results provide a foundation for understanding the genetic architecture of body size traits. They will also contribute to breeding programs and genomic selection work on Chinese Holstein cattle.

The aryl hydrocarbon receptor: A predominant mediator for the toxicity of emerging dioxin-like compounds

The aryl hydrocarbon receptor (AHR) is a member of the basic helix-loop-helix/Per-ARNT-Sim (bHLH-PAS) family of transcription factors and has broad biological functions. Early after the identification of the AHR, most studies focused on its roles in regulating the expression of drug-metabolizing enzymes and mediating the toxicity of dioxins and dioxin-like compounds (DLCs). Currently, more diverse functions of AHR have been identified, indicating that AHR is not just a dioxin receptor. Dioxins and DLCs occur ubiquitously and have diverse health/ecological risks. Additional research is required to identify both shared and compound-specific mechanisms, especially for emerging DLCs such as polyhalogenated carbazoles (PHCZs), polychlorinated diphenyl sulfides (PCDPSs), and others, of which only a few investigations have been performed at present. Many of the toxic effects of emerging DLCs were observed to be predominantly mediated by the AHR because of their structural similarity as dioxins, and the in vitro TCDD-relative potencies of certain emerging DLC congeners are comparable to or even greater than the WHO-TEFs of OctaCDD, OctaCDF, and most coplanar PCBs. Due to the close relationship between AHR biology and environmental science, this review begins by providing novel insights into AHR signaling (canonical and non-canonical), AHR's biochemical properties (AHR structure, AHR-ligand interaction, AHR-DNA binding), and the variations during AHR transactivation. Then, AHR ligand classification and the corresponding mechanisms are discussed, especially the shared and compound-specific, AHR-mediated effects and mechanisms of emerging DLCs. Accordingly, a series of in vivo and in vitro toxicity evaluation methods based on the AHR signaling pathway are reviewed. In light of current advances, future research on traditional and emerging DLCs will enhance our understanding of their mechanisms, toxicity, potency, and ecological impacts.

The role of gut microbial community and metabolomic shifts in adaptive resistance of Atlantic killifish (Fundulus heteroclitus) to polycyclic aromatic hydrocarbons

Altered gut microbiomes may play a role in rapid evolution to anthropogenic change but remain poorly understood. Atlantic killifish (Fundulus heteroclitus) in the Elizabeth River, VA have evolved resistance to polycyclic aromatic hydrocarbons (PAHs) and provide a unique opportunity to examine the links between shifts in the commensal microbiome and organismal physiology associated with evolved resistance. Here, 16S rRNA sequence libraries derived from fish guts and sediments sampled from a highly PAH contaminated site revealed significant differences collected at similar samples from an uncontaminated site. Phylogenetic groups enriched in the libraries derived from PAH-resistant fish were dissimilar to their associated sediment libraries, suggesting the specific environment within the PAH-resistant fish intestine influence the gut microbiome composition. Gut metabolite analysis revealed shifts between PAH-resistant and non-resistant subpopulations. Notably, PAH-resistant fish exhibited reduced levels of tryptophan and increased levels of sphingolipids. Exposure to PAHs appears to impact several bacterial in the gut microbiome, particularly sphingolipid containing bacteria. Bacterial phylotypes known to include species containing sphingolipids were generally lower in the intestines of fish subpopulations exposed to high concentrations of PAHs, inferring a complex host-microbiome relationship. Overall, killifish microbial community shifts appear to be related to a suppression of overall metabolite level, indicating a potential role of the gut in organismal response to anthropogenic environmental change. These results on microbial and metabolomics shifts are potentially linked to altered bioenergetic phenotype observed in the same PAH-resistant killifish populations in other studies.

Evolutionary responses of marine organisms to urbanized seascapes

Many of the world's major cities are located in coastal zones, resulting in urban and industrial impacts on adjacent marine ecosystems. These pressures, which include pollutants, sewage, runoff and debris, temperature increases, hardened shorelines/structures, and light and acoustic pollution, have resulted in new evolutionary landscapes for coastal marine organisms. Marine environmental changes influenced by urbanization may create new selective regimes or may influence neutral evolution via impacts on gene flow or partitioning of genetic diversity across seascapes. While some urban selective pressures, such as hardened surfaces, are similar to those experienced by terrestrial species, others, such as oxidative stress, are specific to aquatic environments. Moreover, spatial and temporal scales of evolutionary responses may differ in the ocean due to the spatial extent of selective pressures and greater capacity for dispersal/gene flow. Here, we present a conceptual framework and synthesis of current research on evolutionary responses of marine organisms to urban pressures. We review urban impacts on genetic diversity and gene flow and examine evidence that marine species are adapting, or are predicted to adapt, to urbanization over rapid evolutionary time frames. Our findings indicate that in the majority of studies, urban stressors are correlated with reduced genetic diversity. Genetic structure is often increased in urbanized settings, but artificial structures can also act as stepping stones for some hard-surface specialists, promoting range expansion. Most evidence for rapid adaptation to urban stressors comes from studies of heritable tolerance to pollutants in a relatively small number of species; however, the majority of marine ecotoxicology studies do not test directly for heritability. Finally, we highlight current gaps in our understanding of evolutionary processes in marine urban environments and present a framework for future research to address these gaps.

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.

Replicated Landscape Genomics Identifies Evidence of Local Adaptation to Urbanization in Wood Frogs

Native species that persist in urban environments may benefit from local adaptation to novel selection factors. We used double-digest restriction-side associated DNA (RAD) sequencing to evaluate shifts in genome-wide genetic diversity and investigate the presence of parallel evolution associated with urban-specific selection factors in wood frogs (Lithobates sylvaticus). Our replicated paired study design involved 12 individuals from each of 4 rural and urban populations to improve our confidence that detected signals of selection are indeed associated with urbanization. Genetic diversity measures were less for urban populations; however, the effect size was small, suggesting little biological consequence. Using an FST outlier approach, we identified 37 of 8344 genotyped single nucleotide polymorphisms with consistent evidence of directional selection across replicates. A genome-wide association study analysis detected modest support for an association between environment type and 12 of the 37 FST outlier loci. Discriminant analysis of principal components using the 37 FST outlier loci produced correct reassignment for 87.5% of rural samples and 93.8% of urban samples. Eighteen of the 37 FST outlier loci mapped to the American bullfrog (Rana [Lithobates] catesbeiana) genome, although none were in coding regions. This evidence of parallel evolution to urban environments provides a powerful example of the ability of urban landscapes to direct evolutionary processes.

No additive genetic variance for tolerance to ethynylestradiol exposure in natural populations of brown trout (Salmo trutta)

One of the most common and potent pollutants of freshwater habitats is 17-alpha-ethynylestradiol (EE2), a synthetic component of oral contraceptives that is not completely eliminated during sewage treatment and that threatens natural populations of fish. Previous studies found additive genetic variance for the tolerance against EE2 in different salmonid fishes and concluded that rapid evolution to this type of pollution seems possible. However, these previous studies were done with fishes that are lake-dwelling and hence typically less exposed to EE2 than river-dwelling species. Here, we test whether there is additive genetic variance for the tolerance against EE2 also in river-dwelling salmonid populations that have been exposed to various concentrations of EE2 over the last decades. We sampled 287 adult brown trout (Salmo trutta) from seven populations that show much genetic diversity within populations, are genetically differentiated, and that vary in their exposure to sewage-treated effluent. In order to estimate their potential to evolve tolerance to EE2, we collected their gametes to produce 730 experimental families in blockwise full-factorial in vitro fertilizations. We then raised 7,302 embryos singly in 2-ml containers each and either exposed them to 1 ng/L EE2 (an ecologically relevant concentration, i.e., 2 pg per embryo added in a single spike to the water) or sham-treated them. Exposure to EE2 increased embryo mortality, delayed hatching time, and decreased hatchling length. We found no population differences and no additive genetic variance for tolerance to EE2. We conclude that EE2 has detrimental effects that may adversely affect population even at a very low concentration, but that our study populations lack the potential for rapid genetic adaptation to this type of pollution. One possible explanation for the latter is that continuous selection over the last decades has depleted genetic variance for tolerance to this synthetic stressor.

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.