Transgenerational inheritance of chemical-induced signature: A case study with simvastatin

Contents lists available at ScienceDirect Environment International journal homepage: www.elsevier.com/locate/envint Transgenerational inheritance of chemical-induced signature: A case study with simvastatin T. Neuparth a,⁎ , A.M. Machado a,b , R. Montes c , R. Rodil c , S. Barros a , N. Alves a , R. Ruivo a , L. Filipe C. Castro a,b , J.B. Quintana c , M.M. Santos a,b,⁎ a CIMAR/CIIMAR—Interdisciplinary Centre of Marine and Environmental Research, University of Porto, Avenida General Norton de Matos, S/N, 4450-208 Matosinhos, Portugal b FCUP –Department of Biology, Faculty of Sciences, University of Porto, Porto, Portugal c Department of Analytical Chemistry, Nutrition and Food Sciences, IAQBUS –Institute of Research on Chemical and Biological Analysis, Universidade de Santiago de Compostela, R. Constantino Candeira S/N, 15782 Santiago de Compostela, Spain ARTICLE INFO Handling Editor: Marti Nadal Keywords: Simvastatin Transgenerational effects Reproduction Transcriptomic Epigenetics Regulatory agencies ABSTRACT The hypothesis that exposure to certain environmental chemicals during early life stages may disrupt reproduction across multiple non-exposed generations has significant implications for understanding disease etiology and adverse outcomes. We demonstrate here reproductive multi and transgenerational effects, at environmentally relevant levels, of one of the most prescribed human pharmaceuticals, simvastatin, in a keystone species, the amphipod Gammarus locusta. The transgenerational findings has major implications for hazard and risk assessment of pharmaceuticals and other contaminants of emerging concern given that transgenerational effects of environmental chemicals are not addressed in current hazard and risk assessment schemes. Considering that the mevalonate synthesis, one of the key metabolic pathways targeted by simvastatin, is highly conserved among metazoans, these results may also shed light on the potential transgenerational effects of simvastatin on other animals, including humans. 1. Introduction A growing body of evidence suggests that the exposure to certain stressors during embryonic development may induce permanent endocrine disturbance, which can be transgenerationally inherited through multiple non-exposed generations, even in the absence of the insult (Chamorro-Garcia et al., 2017;Fraz et al., 2019;Klosin et al., 2017;Vera-Chang et al., 2018;Wang et al., 2016). This may have significant implications for understanding disease etiology and adverse outcomes that can be passed across generations (Anway et al., 2005; Vera-Chang et al., 2018). In the last years, an increasing scientific interest on transgenerational inheritance has emerged. However, despite the growing number of studies claiming transgenerational impacts, many present gaps in the transgenerational experimental design by investigating only direct exposure (F0, F1 and/or F2), excluding the truly non-exposed generations (F3 or above), or testing chemical concentrations far above ecological relevance (Jeremias et al., 2020;Shaw et al., 2017). Here, we exposed the keystone marine amphipod species Gammarus locusta (F0) (Fig. 1A) to environmentally relevant concentrations of the hypocholesterolaemic drug, simvastatin (SIM), and evaluated the transgenerational effect on F3. Statins are one of the most prescribed human pharmaceuticals in western European countries and in the United States, being ubiquitous in aquatic environments in the vicinity https://doi.org/10.1016/j.envint.2020.106020 Received 20 April 2020; Received in revised form 26 June 2020; Accepted 28 July 2020 Abbreviations: AACT, acetoacetyl-CoA thiolase; chiA, chitinase; CYP clan 2, steroidogenic CYP (CYP 2L1); ELOVL6, elongation of very long chain fatty acids protein 6; ELOVL7, elongation of very long chain fatty acids protein 7; FALD, farnesal dehydrogenase; FAMET, farnesoic acid O-methyltransferase; FOLD, farnosol dehydrogenase; FPPP, farnesyl diphosphate pyrophosphatase; FPPS, fanesyl diphosphate synthase; HEXA-B, hexosaminidase; HMGR, hydroxymethylglutaryl-CoA reductase; HMGS, hydroxymethylglutaryl-CoA synthase; IDI1, isopentenyl diphosphate isomerase; JHAMT, juvenile hormone acid methyltransferase; JHEH, juvenile hormone epoxide hydrolase; LSS, lanosterol synthase; MDC, mevalonate-5-decarboxylase; MFE, methyl farnesoate epoxidase; MUFA, monounsaturated fatty acids; MVK, mevalonate kinase; PFK, 6-phosphofructokinase 1; PUFA, Polyunsaturated fatty acids; PVK, phosphomevalonate kinase; SCD, stearoyl-CoA desaturase; SFA, Saturated fatty acids; SQLE, squalene monooxygenase; SQS, squalene synthase; SORD, L-iditol 2-dehydrogenase; UGT, glucuronosyltransferase; xylA, xylose isomerase ⁎ Corresponding authors at: CIMAR/CIIMAR—Interdisciplinary Centre of Marine and Environmental Research, University of Porto, Avenida General Norton de Matos, S/N, 4450-208 Matosinhos, Portugal (T. Neuparth) and FCUP –Department of Biology, Faculty of Sciences, University of Porto, Porto, Portugal (M.M. Santos). E-mail addresses: [email protected] (T. Neuparth), [email protected] (M.M. Santos). Environment International 144 (2020) 106020 Available online 26 August 2020 0160-4120/ © 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). T of urban areas, with concentrations above 100 ng/L being detected in surface waters and up to 2650 and 11,700 ng/L in waste water treatment plant effluents and influents, respectively (Pereira et al., 2015; Santos et al., 2016;Tete et al., 2020). This therapeutic class disrupts cholesterol synthesis by inhibiting the enzyme 3-hydroxy-3-methyl- glutaryl-CoA reductase (HMGR), responsible for the rate limiting step in the mevalonate synthesis, a pathway highly conserved among metazoans, including humans (Fig. 1B). Given that SIM was previously found to impact G. locusta reproduction in the range of concentrations reported in surface waters (Neuparth et al., 2014), and it is a known modulator of epigenetic markers in mammalian cells (Karlic et al., 2015), we investigated the multigenerational and transgenerational (F0 to F3) effects of environmentally relevant concentrations of SIM, using as model species the amphipod G. locusta (Fig. 1C). 2. Materials and methods 2.1. Test species The amphipod Gammarus locusta was the organism selected to investigate the transgenerational effects of Simvastatin (SIM) considering its advantages in comparison with other animal models, i.e., its short life-cycle (35 days at 20 °C), reduced logistical support for testing and breeding and ethical constraints. They are excellent models for toxicity testing and ecotoxicological risk assessment due to their high sensitivity to anthropogenic stressors and wide geographic distribution, but also because they display a key role in the aquatic ecosystem (Costa and Costa, 2000;Neuparth et al., 2002). A deep knowledge on the biology and ecology of the species has been obtained over the past two decades (Costa and Costa, 1999, 2000;Neuparth et al., 2002, 2005, 2014). 2.2. Experimental design All experiments conducted in this study were carried out at Biotério de Organismos Aquáticos (BOGA, CIIMAR) aquatic animal facilities and have been approved by the CIIMAR ethical committee and by CIIMAR Managing Animal Welfare Body (ORBEA) according to the European Union Directive 2010/63/EU on the protection of animals used for scientific purposes. To investigate the multigenerational and the transgenerational (F0 to F3) effects of environmentally relevant concentrations of SIM (CAS n°. 79902-63-9; Sigma Aldrich), G. locusta was continuously maintained under SIM exposure during four successive generations (F0E, F1E, F2E and F3E) –exposed group (Fig. 1C). Additionally, in order to evaluate the putative SIM transgenerational effects, in parallel to the continuous SIM exposure, a sub-set of offspring produced in F0E were raised under SIM-free water (control water) for the subsequent three generations e.g., F1T, F2T, F3T –transgenerational group (Fig. 1C). To conduct the experiment, 7L aquaria in a semi-static system were set-up in controlled laboratory conditions (20 °C, 16 h light/8 h dark). Each aquarium was filled with 5L of natural filtered seawater (33-35‰) dosed with SIM –exposed group –or free of SIM –transgenerational group. A 1 cm-deep layer of natural clean sediment and small stones were supplied to provide shelter and aeration was provided with plastic tips. Each generation (F0 to F3) was initiated with exactly fifty offspring that were allocated to each treatment aquaria, and lasted for 55 to Fig. 1. Rationale of the study. (A) Gammarus locusta in amplexus. (B) Simplified scheme of the Mevalonate pathway. Statins inhibit the enzyme 3-hydroxy-3-methyl- glutaryl-CoA reductase (HMGR), responsible for the rate limiting step in the mevalonate pathway. The dimethylallyl pyrophosphate synthesis from acetyl-CoA is conserved for all metazoan (gray box). The farnesyl pyrophosphatase synthesis (green box) is established in vertebrates and arthropods. Downstream farnesyl pyrophosphatase, there are group-specific variations: cholesterol biosynthesis in vertebrates (blue box) and the sesquisterpenoid biosynthesis in arthropods with the production of methyl farnesoate (yellow box). Beyond cholesterol and methyl farnesoate, other isoprenoids are synthesized in the mevalonate pathway and farnesyl and geranyl moieties also participate in protein prenylation (green Box). (C) Experimental Design: Multigenerational (orange) and transgenerational (blue) studies to investigate the effects of Simvastatin (SIM) on the amphipod G. locusta during 4 consecutive generation (F0 to F3). Multigenerational study (F0E-F3E), transgenerational study (F1T-F3T). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) T. Neuparth, et al. Environment International 144 (2020) 106020 2 65 days. Throughout the experiment amphipods were fed ad libitum with Ulva sp. collected from a coastal area in the Aguda beach, Portugal −41°02′55.2″N 8°39′16.6″W–a site devoid of direct contamination sources. The aquaria were inspected daily for feeding requirements and aeration and the pH, conductivity, dissolved oxygen, ammonia level, salinity and temperature were strictly monitored and controlled at each water change. Two weeks before the beginning of the assay, twenty couples of sexually mature males and females with about 12 week old, from the main culturing system, were separated and used to produce the F0 generation. Each female at this age produce in average 48 offspring (Neuparth et al., 2002). These sub set of adults were maintained at the same temperature and salinity of the bioassay, with unlimited food (Ulva sp.). The bioassay started with fifty, one-week old, offspring per aquarium, randomly isolated from the culturing sub-set and were continually exposed up to adulthood to three environmental concentrations of SIM (32, 64 and 320 ng/L) plus control (filtered natural seawater) and solvent control (0.0005% acetone), four replicate aquaria for each condition (F0E). The concentrations of SIM used in the present study were based upon our previous research (Neuparth et al., 2014) where SIM was found to impact G. locusta reproduction in the range of concentrations reported in surface waters (ng/L range). The SIM solutions for the different treatment groups, 32, 64 and 320 ng/L, were prepared by serial dilutions of a 0.64 mg/mL stock solution prepared in acetone (0.0005%). These solutions were stored in the dark at −20 °C. At each water renewal of the experimental aquaria, that occurred every three days, SIM solutions were applied directly in the middle of each aquarium and the final SIM concentrations were reestablished by properly stirred. This approach was previously validated (Neuparth et al., 2014), thus avoiding the stress of daily water change. Given that in previous experiments SIM was found to be stable in water under a similar experimental design (Neuparth et al., 2014,Barros et al., 2018), the actual concentration of SIM in each treatment/generation was quantified using LC-MS/MS after one of the water change (0 h, SIM addiction) and before the next water change (72 h post SIM addiction) (see Neuparth et al., 2020, Section 2.2). Given that no significant effects in growth and reproduction were observed at the lowest SIM exposure concentration (32 ng/L) (Fig. 2), the treatments at 64 and 320 ng/L were selected to continue the multi and transgenerational studies (F1, F2 and F3 exposed group –64E and 320E; and transgenerational group –64T and 320T). At maturity and after three consecutive reproductions, males and females were sampled (see next section) and the offspring of each treatment were divided into two branches: one was maintained under continuous SIM exposure and raised for more three consecutive generations (F1E, F2E and F3E) and the other was maintained in SIM-free water (natural filtered seawater and SIM free sediments, similar to the control group) for three consecutive generations (F1T, F2T and F3T). 2.3. Sample collection, ecological endpoints In each generation, after the second/third reproduction, the overlying water of each treatment replicate (exposed or transgenerational groups) was sieved through 1000 and 250 µm screens, to collect the adults and their offspring, respectively, following the protocol of Neuparth et al. (2014). Sediments were washed five times to assure that all organisms were collected. Fifty offspring from each treatment replicate were randomly selected to start the subsequent generation. The adults’(males and females) survival, body length, as well as female fecundity, were quantified in each of the four replicate aquaria per condition. Survival and body length were determined separately for each sex and expressed as percentage relative to control groups. Similar to our previous study (Neuparth et al., 2014), survival did not differ among SIM treatments and control groups. After being anesthetized in ice water, the amphipodśbody length was measured to the nearest 0.1 mm using a stereomicroscope (Leica EZ4) by the distance between the anterior end of the rostrum and the posterior end of the last metasomatic segment –metasomatic length (Neuparth et al., 2014). After being measured, adults were individually preserved in RNA later or liquid nitrogen and stored at −80 °C until further use in RNA-seq analyses and methyl farnesoate quantification, respectively. The fecundity (average number of offspring produced by female) was quantified after each reproduction. All offspring present in each aquarium were collected and preserved in 70% ethanol with Bengal Rose for later quantification of reproductive success. 2.4. Quantification of methyl farnesoate levels in Gammarus locusta The concentration of methyl farnesoate (MF) in G. locusta was determined for each treatment and sex of exposed and transgenerational groups in the four generations by matrix-solid phase dispersion (MSPD) and gas chromatography-chemical ionization mass spectrometry (GCCI-MS). This analytical methodology was developed as an alternative to hemolymph extraction and is described in detail in Montes et al. (2017). The MF levels (ng/g) were determined in four male and four female samples, each composed by three pooled individuals belonging to different replicates of each treatment of the exposed and transgenerational groups from the four generations. The sample weight was recorded, and the results were expressed as wet weight concentration. Deuterated methyl heptadecanoate (Sigma-Aldrich) was used as internal standard (IS). Each sample (average weight 0.302 g and 0.147 g for males and females, respectively) was dispersed in a mortar with 0.5 g of PSA (primary-secondary amine) and loaded into a MSPD cartridge previously filled with 1.5 g of Florisil (co-sorbent) and 0.5 g of sodium sulfate anhydrous (desiccant). Then, MF and IS were eluted from the MSPD cartridges using 1.5 mL of ethyl acetate. The final extracts were concentrated to 20 µL and injected in the GC-CI-MS system (Agilent) using isobutane as reagent gas. Both MF and the IS were measured in the single ion monitoring mode using the protonated molecular ion for quantification (Q) and two additional ions as qualifiers (q), per compound. These ions (m/z) were 251 (Q) and 219, 191 (q) for MF and 317.5 (Q) and 267.5, 282.5 (q) for IS. The concentration of the extracts was calculated by internal standard calibration in a concentration range between 5 and 250 ng/mL and then expressed as wet weight concentration (ng/g) for each sample. The MQLs were 3.1 ng/g for males 2 ng/g for females. The recoveries were evaluated at two different spiked levels (2 and 10 ng for females and 1 and 5 ng for males, over the native concentration) and turned out 96 and 103% for females and 81 and 85% for males with a maximum relative standard deviation (RSD) of 16%. 2.5. Statistical analysis For each generation (F0 to F3), the ecological data (males and femaleślength and fecundity) and the methyl farnesoate levels from the exposed and transgenerational groups were first checked for normality (Kolmogorov-Smirnov test) and homogeneity of variances (Levene's test). If the ANOVA assumptions were fulfilled, each endpoint data from the exposed and transgenerational groups were analyzed together by one-way analysis of variance (ANOVA). When significant differences were found by ANOVA, Fisher’s least significant difference (LSD) test was used to discriminate statistically significant differences between exposed or transgenerational treatments and the control. For males and females length, as the normality of the data were not fulfilled, even after data transformation, a Kruskal-Wallis ANOVA by ranks was used. Significant differences were established at p < 0.05 for every analysis and all statistics were computed with Statistica 12.5 (Stat-soft, USA). Similar to our previous study (Neuparth et al., 2014), the solvent used, acetone at 0.0005%, rendered no significant effect (P > 0.05) in comparison to the water control at any of the generations/sex. Therefore, following Green and Wheeler (2013) proposal, the data from these groups were pooled to increase statistical power and referred to as T. Neuparth, et al. Environment International 144 (2020) 106020 3 control group. 2.6. RNA extraction, library preparation and Illumina sequencing of G. locusta samples RNA-seq was performed independently in three individual G. locusta females, per group, over the four generations: solvent control groups (F0.C, F1.C, F2.C and F3.C), F0 320 exposed group (F0.320E) and F1 to F3 320 transgenerational group (F1.320T; F2.320T and F3.320T) (Detailed information displayed in Neuparth et al., 2020, Annex 1). Total RNA from three females per group was extracted using the Illustra RNAspin Min RNA Isolation Kit (GE Healthcare), according to the manufacturer’s protocol. RNA sample integrity and concentration were initially checked by agarose gel electrophoresis and then by an Agilent Bioanalyzer 2100 (Agilent Technologies, USA), with RIN values varying from 6.3 to 7.7. RNA-seq libraries and sequencing, Illumina HiSeq2500 paired-end (2x150), were obtained commercially at Novogene (Hong Kong). 2.7. In silico clean-up of raw datasets and de novo transcriptome assembly The 563 M of PE sequenced reads were analysed using several gold standard methods. Initially, to have a quality overview of the datasets, the FastQC (v.0.11.8) software was applied (http://www. bioinformatics.babraham.ac.uk/projects/fastqc/). After that, the Rcorrector (v.1.0.3) (Song and Florea, 2015) was applied for error correction, with the defaults, and Centrifuge (v.1.0.3-beta) (Kim et al., 2016) to taxonomically filter each read dataset. To perform the reference transcriptome assembly of G. locusta the de novo assembly method Trinity (v.2.8.4) (Grabherr et al., 2011) was applied (Detailed information displayed in Neuparth et al., 2020, section 2.3). 2.8. Decontamination, open reading frame (ORF) prediction, and assessment of transcriptome assembly The transcriptome assembly decontamination and assessments were done using three distinct approaches, blast searches against NCBI Fig. 2. Multigenerational and transgenerational impairment of simvastatin across multiple generations (F0 to F3). (A) Growth, expressed as percentage relative to control group, n = 52–58. (B) Fecundity, expressed as percentage relative to control group, n = 4 replicate tanks per treatment. (C) methyl farnesoate (MF) levels (ng/g) expressed as percentage relative to control, n = 4 pools for males and females. Data are presented as the means ± standard error and the asterisks represent significant difference in the exposed group (E) or transgenerational group (T) compared with the control (CTR): *p < 0.05 and **p ≤0.01. T. Neuparth, et al. Environment International 144 (2020) 106020 4 databases, calculation of basic statistics with Trinity and Transrate scripts (v.1.0.3) (Smith-Unna et al., 2016) and comparisons with databases of conserved orthologs genes, through the Benchmarking Universal Single –Copy Orthologs software (BUSCO v.3.0.2) (Simão et al., 2015). The open reading frames were predicted using the TranDecoder software (v.5.3.0) (Haas et al., 2013). This process generated three versions of the transcriptome assembly; the 1st version –Raw transcriptome assembly –Decontaminated; 2nd version –Protein coding assembly –All transcripts with ORF; 3rd version –Protein coding assembly –Unigenes with ORF (Detailed information displayed in Neuparth et al., 2020, section 2.4). 2.9. Transcriptome annotation To perform the transcriptome annotation, the second version of the transcriptome assembly and Trinotate tool (v.3.1.1) were used (Bryant et al., 2017). Trinotate uses several methods to functionally annotate transcriptome assemblies. Here, a bunch of tools were used to annotate the transcriptome of G. locusta: blast searches, identification of protein domains (PFAM Database) (39), classification of orthologous groups of genes (eggNOG) (Powell et al., 2012) and attribution of Gene Ontology (GO) terms (Ashburner et al., 2000) (Detailed information displayed in Neuparth et al., 2020, section 2.5). 2.10. Differential gene expression analyses For the differential gene expression (DEGs) analyses the second version of transcriptome assembly and the clean reads were used. The second version of the transcriptome was used to maximize mapping of the raw reads against the protein coding transcriptome. The read mapping and transcript quantification were done with the abundance_estimates_to_matrix.pl script of Trinity the pipeline, under the defaults (Haas et al., 2013). Here, the Bowtie2 (v.2.3.5) (Langmead and Salzberg, 2012) software was applied to map the reads and the RSEM (v.1.3.0) (Li and Dewey, 2011) tool to estimate the transcript abundance. The differential gene expression analyses were carried out in Degust (v.4.1.1) (Powell et al., 2015) platform (http://degust.erc. monash.edu/) being the edgeR (v.3.26.8) package (Robinson et al., 2010) of R (v.3.6.1) used to calculate differential gene expression. To determine DEG, the exposed or transgenerational samples of each generation were compared against the respective control (F0.C vrs F0.320E; F1.C vrs F1.320T; F2.C vrs F2.320T; F3.C vrs F3.320T), and all genes with False Discovery Rate –corrected (FDR) p-value < 0.05 and log2|fold change| ≥2 were considered differentially expressed (Detailed information displayed in Neuparth et al., 2020, section 2.6). 2.11. KEGG pathways analyses The main metabolic pathways affected by the transgenerational effects of SIM across all generations were identified using the KAAS webserver (Moriya et al., 2007). The DGEs genes were blasted onto KEGG pathways database using the single-directional best hit (SBH) method, the proteins corresponding to the differential expressed genes (collected from the third version of the transcriptome), and 710,890 reference sequences of 40 manually selected species (Detailed information displayed in Neuparth et al., 2020, section 2.7). 3. Results and discussion 3.1. Multigenerational and transgenerational impairment of simvastatin across generations: effects on G. locustagrowth and reproduction and methyl farnesoate (MF) levels We first exposed G. locusta juveniles (one-week old) from our laboratory culture to environmentally relevant concentrations of SIM (32, 64 and 320 ng/L) (Neuparth et al., 2020, Table A) up to adulthood (F0) and evaluated ecologically relevant endpoints (female fecundity and growth). At the F0 generation, we found significant effects at 320 ng/L, i.e., decrease of reproductive output and growth, and reduction of female growth at 64 ng/L (Fig. 2A and B). We then evaluated these endpoints in the three subsequent generations (F1 to F3), both in animals continuously exposed to SIM (exposed group, 64E and 320E) and in the offspring of animals exposed in F0 thereafter maintained in control water (transgenerational group, 64T and 320T). We found a significant impact in fecundity, in several treatments in comparison with control, up to F3. A significant decrease in reproduction was observed for 64E, 320E and 320T (p < 0.05) in the F1 and for 320E (p < 0.01) in the F2. In the F3, the transgenerational and exposed groups (64T and 320T and 320 E) also presented a significant decrease in the reproductive output i.e., 40% decreased fecundity for 64T, 50% for 320T and 48% for 320E (Fig. 2B). Growth was also impacted in the exposed group (320E) and in the transgenerational group (320T), but in the latter only in F2 and F3 (Fig. 2A). Although 64T displayed a significant impact in the F3 fecundity, animals from this treatment did not present, in any of the generations or sexes, growth impairment. This supports the hypothesis that the fecundity decrease observed for 64T in F3 is independent from female length. The final steps of the mevalonate pathway display group-specific variations, i.e., cholesterol synthesis in vertebrates and synthesis of the isoprenoid methyl farnesoate (MF) in arthropods (Fig. 1B; Santos et al., 2016). Given that in crustaceans MF has a central role in molting and reproduction (Laufer and Biggers, 2001), we next developed a highly sensitive analytical method (matrix solid-phase dispersion with gas chromatography coupled to mass spectrometry) to detect MF in arthropod tissues, including G. locusta (Montes et al., 2017). Both males and females from all generations (F0 to F3) were screened for MF tissue levels, with this isoprenoid being detected in all samples. Levels of MF in males were significantly increased at 32E and 320E in F0 and 64T in F1 (Fig. 2C). For F2 and F3, both in the exposed and the transgenerational groups, male MF levels did not differ significantly from control (P > 0.05) (Fig. 2C). In contrast, a significant decrease of MF levels was recorded in females for most of the treatments and generations (Fig. 2C). With exception for the females from 320E in F0 and F3 that did not differ significantly from control, female MF levels were significantly reduced in all F0, F1 and F3 groups with a decrease in comparison to control ranging between 27 and 45%. In F2, the exception was the females from 64T and 320T. In fact, in F2, similar to females’fecundity, levels of female MF in the transgenerational groups (64T and 320T) did not differ significantly from control. In F3, females from both transgenerational groups (64T and 320T) presented a significant reduction on MF levels in comparison with control, with the decrease ranging from 25% (64T, p < 0.05) to 35% (320T, p < 0.01). These findings suggest that the transgenerational inheritance of the reproduction impairment in the 64T and 320T groups could be associated with a decrease of the MF levels in females, a key regulator of reproduction in crustaceans. In our previous study to address the chronic life-cycle effect of SIM in G. locusta reproduction (Neuparth et al., 2014), we observed a similar impact to that reported here in F0 females exposed to SIM at 320 ng/L. In order to further understand the effects in reproduction, Neuparth et al., 2014 performed a detailed histological analysis of male and female gonads. Whereas female gonad life-cycle exposed to 320 ng/L displayed delayed maturation revealing predominance of primary oocytes, the majority of control females were fully mature. In contrast, spermatogenesis was not affected at SIM concentrations ranging from 320 ng/L to 8000 ng/L. Hence, taken together, our previous findings (Neuparth et al., 2014) and the results reported here for MF levels and fecundity support the hypothesis that SIM effects on reproduction seem to be associated with females rather than males. T. Neuparth, et al. Environment International 144 (2020) 106020 5 3.2. Transcriptomic analysis of females across generations To gain additional insights into the biological functions and canonical pathways transgenerationally disrupted by SIM a comprehensive transcriptome assembly was produced and the transcriptomic profiles of females from the control group (F0 to F3), the 320E group (F0) and the 320T group (F1 to F3) were analyzed (Fig. 3,Fig. 4 and Neuparth et al., 2020, Annex 1 and Annex 2). In order to mine the large set of multigenerational data, we focused on the interconnections and the network of genes sharing common biological processes and functions rather than on single genes. We addressed in particular the common gene network pathways affected in females from F0 (320E) and F3 (320T), that can contribute to the understanding of the transgenerational impaired fecundity and growth on F3 (Fig. 3 and Fig. 4). The RNA sequencing (RNA-seq) analysis (Neuparth et al., 2020, Annexes 1 to 8, 12 to 14 and 16) identified a total of 3785 differentially expressed unigenes in females across all generations (F0 to F3), with a False Discovery Rate –corrected (FDR) p-value < 0.05 and log2|fold change| ≥2(Neuparth et al., 2020, Annexes 4 to 7). Of these, 1482 and 2017 genes were altered in F0 and F3, respectively (Fig. 3C and Neuparth et al., 2020, Annex 3 and Annex 11). Hierarchical clustering heatmap (Fig. 3A and Neuparth et al., Annex 9) revealed a distinct pattern of gene transcription between control and 320E or 320T, in both F0 and F3, which is consistent with the multidimensional scaling analysis (Fig. 3B and Neuparth et al., 2020, Annex 10) that grouped control and exposed females in different axes. The Venn diagram corroborates this pattern identifying over 800 DEGs common to both F0 and F3 (Fig. 3C and Neuparth et al., 2020, Annex 11). We further investigated the most prominent overrepresented canonical pathways using KEGG Mapper –reconstruct pathway. Twenty four pathways contained DEGs in F3 generation (320T versus control): metabolism of terpenoids and polyketides (mevalonate pathway/sesquiterprenoid biosynthesis and degradation, ecdysteroid biosynthesis), carbohydrate metabolism, lipid metabolism, energy metabolism, amino acid metabolism, glycan biosynthesis and metabolism, metabolism of cofactors and vitamins, xenobiotics biodegradation and metabolism, transcription, folding sorting and degradation, signal transduction, signaling molecules and interaction, transport and catabolism, cell growth and dead, cell motility, immune system, endocrine system, digestive system, circulatory system, nervous system, sensory system, development and regeneration, aging and environmental adaptation (Fig. 4 and Neuparth et al., 2020, Annex 7 and Annex 15). Importantly, all F3 DEGs found in these pathways, between 320T and control, were also differentially expressed in F0 females exposed to SIM (320E) (Neuparth et al., 2020, Annexes 4, 8 and 15). We next evaluated whether impaired MF levels and reproduction caused by exposure to SIM in F0 were associated with shifts in the transcriptional profiles of mevalonate biosynthesis modules, as well as the downstream sesquiterpenoid biosynthesis and degradation and/or the ecdysteroid biosynthesis in F3 (Fig. 4A and Neuparth et al., 2020, Annex 16). Three genes were differentially expressed in F3 320T in comparison with control females: methyl farnesoate epoxidase (MFE), juvenile hormone epoxide hydrolase (JHEH) and CYP-2L1 (Fig. 4A). Interestingly JHEH that catalyzes an important step of sesquiterpenoid degradation is modulated in F0 and F3 (Fig. 4A). In insects, JHEH is involved in the regulation of juvenile hormone (JH) titer and metamorphosis/molt (Sin et al., 2015). Although JH has not been identified in non-insect arthropods, several gene orthologues involved in JH production and degradation, as JHEH, previously thought to be insect-specific, were identified in water flea, shrimp and centipede (Guittard et al., 2011;Lee et al., 2011;Sin et al., 2015). In crustaceans, similarly to insects, the degradation of MF appears to occur through ester hydrolysis by specific carboxylesterases (Lee et al., 2011). Therefore, as crustaceans lack JH, the potential substrate for the hydrolysis reaction of JHEH could be the MF (Lee et al., 2011). The CYP-2L1, a steroidogenic enzyme from the CYP clan 2, was also found to be differentially expressed in both F0 and F3 in the ecdysteroid biosynthesis pathway. Although this enzyme has not been functionally characterized in crustaceans, the function of CYP18A1, belonging also to the CYP clan 2, has been characterized in Drosophila (Guittard et al., 2011), displaying a central role in the maintenance of ecdysteroids within physiological levels, regulating Fig. 3. Transcriptomic analysis. (A) Hierarchical clustering heatmap depicting the patterns of gene transcription between F0 control (F0.C.R), 320E (F0.320E.R), F3 control (F3.C.R) and 320T (F3.320T.R), n = 3 per treatment; p-value < 0.05 and fold change ≥2. (B) multidimensional scaling analysis between F0 and F3 (C) Venn diagram of the overlapping significantly (p < 0.05) modulated common genes to F0 and F3. T. Neuparth, et al. Environment International 144 (2020) 106020 6 ecdysteroids inactivation, thus affecting molting and metamorphosis (Guittard et al., 2011). Additionally, the ecdysone inducible gene (E75), a primary target of ecdysone receptor with a critical role in the molting process of arthropods (Priya et al., 2010), by activating the degradation of old cuticle and the formation of a new cuticle, was found to be downregulated in the F3 (Neuparth et al., 2020, Annex 7). In amphipods, females’reproductive cycles are correlated with the molt cycle and reproduction is preceded by a molt (Hyne, 2011). The molt cycle is regulated by several hormones (including MF), under direct control of the steroid hormone, 20-hydroxyecdysone, the putative substrate of CYP-2L1. Thus, disruption of CYP-2L1 transcription could also affect this process. The transcription of genes coding for other enzymes, i.e., phosphofructokinase (PFK), chitinase (CHIT1) and elongation of very long Fig. 4. Schematic representation of selected metabolic pathways differentially expressed in females from F0 (320E versus control) and F3 (320T versus control). (A) Carbohydrate metabolism (Glycolysis/gluconeogenesis, Pentose phosphate pathway, Pentose and glucoronate interconversions, Fructose and mannose metabolism, Galactose metabolism and Amino sugar and nucleotide sugar metabolism), Metabolism of terpenoids and polyketides (Mevalonate pathway/Sesquiterpenoid biosynthesis and degradation, and Ecdysteroid biosynthesis). (B) Lipid metabolism (Biosynthesis of unsaturated fatty acids). The genes modulated by SIM according to KEGG annotation are represented in rectangles. T. Neuparth, et al. Environment International 144 (2020) 106020 7 chain fatty acids protein 7 (ELOVL 7), with a key function in carbohydrate and lipid metabolism, were found to be disrupted in the F3 between 320T and control (Fig. 4A and B and Neuparth et al., 2020, Annex 7). PFK, ELOVL 7 and several others genes were also impacted in the F0 (320E versus control) (Fig. 4A and B). The findings of changes on carbohydrate and lipid metabolism integrate well with our previous histological observation (Neuparth et al., 2014): females’life-cycle exposed to SIM concentrations ranging from 320 to 8000 ng/L yielded reduced glycogen and lipid storage. Interestingly, a recent study with short-term exposure of mussel (Mytilus edulis) to atorvastatin at 1200 ng/L, led to depletion of lipids and carbohydrates (Falfushynska et al., 2019). Therefore, reproductive impairment and/or growth may also be linked to metabolic dysfunction. Additional pathways modulated by SIM relate to the cuticle metabolism. The transcription of twenty genes were down-regulated in F3 (Neuparth et al., 2020, Annex 12). Interestingly, most of these genes altered in F3 were also down-regulated in F0 (Neuparth et al., 2020, Annex 12). This differential gene expression profile suggests that SIM transgenerational exposure modulates exoskeleton maintenance and molting. Since the cuticle degradation process is closely linked to the molt cycle and synchronized with amphipod’s growth and reproduction (Trapp et al., 2014), the differential expression of genes coding for cuticle proteins could also be linked to the observed effects on growth and reproduction of the F3 exposed and transgenerational groups (Fig. 2). SIM exposure also induced the transcription modulation of several processes involved in epigenetic regulation, notably, histone posttranslational modifications such as methylation, acetylation, crotonylation and sumonylation, methylation index modification and nucleosome assembly (Table 1 and Neuparth et al., Annex 13). Interestingly, the protein FAM98A (FAM98A) and YEATS domain-containing protein 2-like (YEATS2) were down-regulated in 320T females from F0, F1 and F3 generations (Table 1 and Neuparth et al., 2020, Annex 13). Since these pathways are implicated in epigenetic modifications, their modulation may be associated with the observed transgenerational disruption of the reproduction and growth in F3 (Aristizaba et al., 2019). Previous studies with mammalian cells demonstrated that SIM modulates DNA methylation, histone modifications and miRNAs (Karlic et al., 2015), thus supporting the transcriptomic responses observed here. Given the complex network of pathways involved in statin-in- duced epigenetic modifications (Karlic et al., 2015), future experiments are needed to unravel the molecular mechanisms underlying SIM modulation of epigenetic markers and how these relate with the observed fingerprint between F0 and F3. 3.3. Implications for hazard and risk assessment Overall, this study provides the first evidence that one of the most prescribed pharmaceuticals in western European countries, at environmental relevant concentrations, induces transgenerational effects in the reproduction of a keystone marine amphipod species. This is very relevant, since the effects of the SIM exposure experienced by G. locusta during its lifetime in the parental generation (F0E) were truly inherited to non-exposed generations (F3T) and were not a simply result of direct exposure. These results have broad implications for our understanding of the hazard of transgenerational impairment and how environmental pollutants may contribute to such disruption. The current chemical risk assessment paradigm includes steps of hazard identification, dose-re- sponse assessment, exposure assessment, and risk characterization (Shaw et al., 2017). Available test guidelines already address multigenerational effects, but, the potential transgenerational effects of environmental pollutants are still not integrated in hazard and risk assessment frameworks. Thus, the results reported here highlight the importance of integrating chemical-induced transgenerational effects in toxicity testing strategies and risk assessment. CRediT authorship contribution statement T. Neuparth: Conceptualization, Methodology, Validation, Formal analysis, Investigation. :Resources, Data curation, Writing - original Fig. 4. (continued) T. Neuparth, et al. Environment International 144 (2020) 106020 8 draft, Writing - review & editing, Supervision, Project administration. A.M. Machado: Methodology, Validation, Formal analysis, Investigation. :Data curation, Writing - review & editing. R. Montes: Methodology, Formal analysis, Investigation, Data curation. R. Rodil: Formal analysis, Investigation, Data curation. S. Barros: Formal analysis, Investigation, Data curation. N. Alves: Formal analysis, Investigation, Data curation. R. Ruivo: Formal analysis, Investigation, Data curation, Writing - review & editing. L. Filipe C. Castro: Validation, Formal analysis, Investigation, Resources, Data curation, Writing - review & editing. J.B. Quintana: Formal analysis, Investigation, Data curation, Writing - review & editing, Funding acquisition. M.M. Santos: Conceptualization, Methodology, Validation, Formal analysis, Investigation, Resources, Data curation, Writing - original draft, Writing - review & editing, Supervision, Project administration, Funding acquisition. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments This research was funded by COMPETE 2020, Portugal 2020, the European Union through the ERDF and the Portuguese Foundation for Science and Technology –FCT (Transobesogen project –Trans-phyletic obesogenic responses: from epigenetic modules to transgenerational environmental impacts, reference: PTDC/CTA-AMB/31544/2017 – NORTE-01-0145-FEDER-031544). This research was also supported by the National Funds through FCT under the projects (UIDB/04423/2020; UIDP/04423/2020), by the Spanish Agencia Estatal de Investigación (CTM2017-84763-C3-2-R), and by the Galician Council of Culture, Education and Universities (ED431C2017/36), cofounded by ERDF. A PhD grant awarded to Susana Barros (PD/BD/143090/2018) was funded by the FCT. Data availability All versions of the transcriptome were submitted to the online Figshare repository and can be consulted in the following link: https:// figshare.com/s/1110e0d14fcc6a275acb (available if requested). All the clean read datasets were submitted to Sequence Read Archive (SRA) database of NCBI, and can be consulted under the Bio project number PRJNA600472. References Anway, M.D., Cupp, A.S., Uzumcu, M., Skinner, M.K., 2005. Epigenetic transgenerational actions of endocrine disruptors and male fertility. Science 308, 1466–1469. https:// doi.org/10.1126/science.1108190. Aristizaba, M.J., Anreiter, I., Halldorsdottir, T., Odgers, C.L., McDade, T.W., Goldenberg, A., et al., 2019. Biological embedding of experience: a primer on epigenetics. PNAS 201820838. https://doi.org/10.1073/pnas.1820838116. Ashburner, M., Ball, C.A., Blake, J.A., Botstein, D., Butler, H., Cherry, J.M., et al., 2000. Gene ontology: tool for the unification of biology. Gene Ontol. Consortium. Nat. Genet. 25, 25–29. https://doi.org/10.1038/75556. Barros, S., Montes, R., Quintana, J.B., Rodil, R., André, A., Capitão, A., Soares, J., Santos, M.M., Neuparth, T., 2018. Chronic environmentally relevant levels of simvastatin disrupt embryonic development, biochemical and molecular responses in zebrafish (Danio rerio). Aquat. Toxicol. 201, 47–57. https://doi.org/10.1016/j.aquatox.2018. 05.014. Bryant, D.M., Johnson, K., DiTommaso, T., Tickle, T., Couger, M.B., Payzin-Dogru, D., et al., 2017. A Tissue-Mapped axolotl de novo transcriptome enables identification of limb regeneration factors. Cell. Rep. 18, 762–776. https://doi.org/10.1016/j.celrep. 2016.12.063. Chamorro-Garcia, R., Diaz-Castillo, C., Shoucri, B.M., Käch, H., Leavitt, R., Shioda, T., Blumberg, B., 2017. Ancestral perinatal obesogen exposure results in a transgenerational thrifty phenotype in mice. Nat. Commun. 8, 2012. https://doi.org/10. 1038/s41467-017-01944-z. Costa, F.O., Costa, M.H., 1999. Life history of the amphipod Gammarus locusta fromthe Sado estuary (Portugal). Acta Oecol. 20, 305–314. https://doi.org/10.1016/S1146- 609X(99)00136-8. Costa, F.O., Costa, H.H., 2000. Review of the ecology of Gammarus locusta (L.). Pol. Arch. Hydrobiol. 48, 541–559. Falfushynska, H., Sokolov, E.P., Haider, F., Oppermann, C., Kragl, U., Rut, W., et al., Table 1 Overview of the genes involved in epigenetic regulation altered in F0 to F3. F0/F3 indicates genes common to both F0 and F3; Red indicates upregulation, blue indicates downregulation; (-) indicates no changes in gene expression; DDB – putative methyltransferase DDB; FAM98A –protein FAM98A; YEATS2 –YEATS domain-containing protein 2-like; PRMT1 –protein arginine N-methyltransferase 1- like; NAA40 –N-alpha-acetyltransferase 40-like isoform X1; AHCY –adenosylhomocysteinase; INO80C –INO80 complex subunit C; HN1N2 –histone-binding protein N1/N2-like isoform X1; SAE1 –SUMO-activating enzyme subunit 1-like; SAM –S-adenosyl methionine; SAH –S-adenosyl homocysteine; HCY –homocysteine; Ado – adenosine. T. Neuparth, et al. Environment International 144 (2020) 106020 9