Anthony Bertucci
email: anthony.bertucci@u-bordeaux.fr

Fabien Pierron

Tao Ye

Patrice Gonzalez

Patrice Couture

Magalie Baudrimont

Identification and expression of microRNAs in european eels Anguilla anguilla from two natural sites with different pollution levels

Keywords: MicroRNA, Anguilla anguilla, Aquatic ecotoxicology, Epigenetics, Biomarker

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

MicroRNAs (miRNAs) are a class of small non-coding RNA [START_REF] Großhans | Molecular biology: The expanding world of small RNAs[END_REF][START_REF] Yasuda | The RNA continent[END_REF] . These 20-24 nucleotides-long sequences associate with the 3'-untranslated region (3'-UTR) of target messenger RNAs (mRNAs), and post-transcriptionally regulate the expression of numerous genes by mediating translational repression or mRNA degradation [START_REF] Bartel | MicroRNAs: Genomics, Biogenesis, Mechanism, and Function[END_REF][START_REF] Chekulaeva | Mechanisms of miRNA-mediated post-transcriptional regulation in animal cells[END_REF] . Owing to this ability to regulate gene expression without involving DNA sequence changes, miRNAs are hence considered as epigenetic factors like histone modification and DNA methylation. In mammals, more than 50% of mRNAs are predicted to be the subject of miRNA-mediated control [START_REF] Friedman | Most mammalian mRNAs are conserved targets of microRNAs[END_REF] . One miRNA may regulate hundreds of target mRNAs, and one mRNA may contain multiple binding sites for multiples miRNAs [START_REF] Lim | Microarray analysis shows that some microRNAs downregulate large numbers of target mRNAs[END_REF] , thus resulting in a complex regulatory network. Although miRNAs are involved in regulation of almost all cellular processes, such as development, growth, apoptosis, immunity and maintenance of tissue-specific function [START_REF] Bushati | microRNA functions[END_REF] , mechanistic aspects of this regulation are not fully understood. In mammals, the aberrant expression of miRNAs has been linked to various diseases and cancer [START_REF] Huang | MicroRNA Machinery Genes as Novel Biomarkers for Cancer[END_REF][START_REF] Cho | MicroRNAs in cancer -from research to therapy[END_REF] and to the inheritance of stress-induced phenotypes across generations [START_REF] Gapp | Implication of sperm RNAs in transgenerational inheritance of the effects of early trauma in mice[END_REF] . Toxic environmental factors such as nanoparticles, organic pollutants and metals can alter miRNA expression [START_REF] Yu | The role of microRNAs in toxicology[END_REF][START_REF] Hou | Environmental chemicals and microRNAs[END_REF] . These changes in miRNAs expression may be part of the adaptive responses of organisms to pollutants and of the mechanisms of toxicity and disease aetiology.

Advances and accessibility of sequencing technologies, coupled with the development of microRNAs dedicated tools are improving the sensitivity of analyses and the ability to detect lowly abundant small RNAs [START_REF] Creighton | Expression profiling of microRNAs by deep sequencing[END_REF][START_REF] Friedländer | Discovering microRNAs from deep sequencing data using miRDeep[END_REF][START_REF] Li | Performance comparison and evaluation of software tools for microRNA deepsequencing data analysis[END_REF][START_REF] Eminaga | Quantification of microRNA Expression with Next-Generation Sequencing[END_REF] even in non-model species. The latest release of the miRNA database (miRBase 21) contains 28,645 hairpin precursors from 223 species [START_REF] Kozomara | miRBase: annotating high confidence microRNAs using deep sequencing data[END_REF] . In aquatic organisms, the role of miRNAs has been investigated in response to various environmental stresses, such as immune challenge, metals, salinity, hypoxia or acidic water [START_REF] Andreassen | Identification of differentially expressed Atlantic salmon miRNAs responding to salmonid alphavirus (SAV) infection[END_REF][START_REF] Wang | Copper-induced deregulation of microRNA expression in the zebrafish olfactory system[END_REF][START_REF] Wang | MicroRNA-Sequence Profiling Reveals Novel Osmoregulatory MicroRNA Expression Patterns in Catadromous Eel Anguilla marmorata[END_REF][START_REF] Zhang | Integrated analysis of mRNA-seq and miRNA-seq in the liver of Pelteobagrus vachelli in response to hypoxia[END_REF][START_REF] Zhou | The identification and characteristics of immune-related microRNAs in haemocytes of oyster Crassostrea gigas[END_REF][START_REF] Lau | Identification and expression profiling of microRNAs in the brain, liver and gonads of marine medaka (Oryzias melastigma) and in response to hypoxia[END_REF][START_REF] Kure | Molecular responses to toxicological stressors: profiling microRNAs in wild Atlantic salmon (Salmo salar) exposed to acidic aluminum-rich water[END_REF] .

As a catadromous fish species, the European eel Anguilla anguilla (Anguillidae; Teleostei) exhibits a complex life cycle including marine (spawning, larval phase and sexual maturation) and continental (somatic growth) environments. Moreover, A. anguilla is a panmictic species and its unique life cycle makes it particularly vulnerable to pollution. After spawning in the Sargasso Sea, eel larvae drift back toward European and North African continental waters following oceanic currents such as the Gulf Stream and the North Atlantic Drift. There, they first metamorphose into glass eel (unpigmented) before reaching the juvenile growth phase (yellow eel) in continental water bodies where they remain up to 15-20 years, often moving across different habitat types. This stage ends with a second metamorphosis and maturing eels (silver eels) will migrate back to the Sargasso Sea to reproduce and die [START_REF] Feunteun | Management and restoration of European eel population (Anguilla anguilla): An impossible bargain[END_REF][START_REF] Guimarães | Yellow eel (Anguilla anguilla) development in NW Portuguese estuaries with different contamination levels[END_REF] . For their unusual life cycle and vulnerability to pollution, eels are considered as sentinel species in ecotoxicology [START_REF] Belpaire | Pollution fingerprints in eels as models for the chemical status of rivers[END_REF][START_REF] Geeraerts | The effects of contaminants in European eel: a review[END_REF] . Moreover, the species has been considered as critically endangered since 2008 (IUCN red list) [START_REF] Iucn | The IUCN Red List of Threatened Species 2014[END_REF] . Its catadromous migratory behaviour, long life, serious habitat reduction, migration barriers, pollution, human-introduced diseases, overfishing, as well as climatic events may be amongst the causes of the catastrophic collapse of the European eel population observed over the past decades [START_REF] Knights | A review of the possible impacts of long-term oceanic and climate changes and fishing mortality on recruitment of anguillid eels of the Northern Hemisphere[END_REF][START_REF]Spawning Migration of the European Eel[END_REF][START_REF] Bevacqua | A global viability assessment of the European eel[END_REF] . A better understanding of the biology of this species may come from the use of high-throughput (next generation) sequencing technologies that have recently provided transcriptomic [START_REF] Coppe | Sequencing, de novo annotation and analysis of the first Anguilla anguilla transcriptome: EeelBase opens new perspectives for the study of the critically endangered European eel[END_REF][START_REF] Baillon | Transcriptome profile analysis reveals specific signatures of pollutants in Atlantic eels[END_REF][START_REF] Pujolar | Surviving in a toxic world: transcriptomics and gene expression profiling in response to environmental pollution in the critically endangered European eel[END_REF][START_REF] Podgorniak | Differences in brain gene transcription profiles advocate for an important role of cognitive function in upstream migration and water obstacles crossing in European eel[END_REF] and genomic [START_REF] Henkel | Primitive Duplicate Hox Clusters in the European Eel's Genome[END_REF][START_REF] Minegishi | Genomics in Eels -Towards Aquaculture and Biology[END_REF][START_REF] Callol | An Enriched European Eel Transcriptome Sheds Light upon Host-Pathogen Interactions with Vibrio vulnificus[END_REF] information. However, like in many non-model species, the non-coding portion of these data has been overlooked in spite of the emerging importance of their role in gene expression regulation. So far, in A. anguilla, microRNAs were only predicted from de novo assembled transcriptome [START_REF] Coppe | Sequencing, de novo annotation and analysis of the first Anguilla anguilla transcriptome: EeelBase opens new perspectives for the study of the critically endangered European eel[END_REF] while in the related species A. marmorata (marbled eel), microRNAs were specifically sequenced and their expression pattern assessed in response to salinity [START_REF] Wang | MicroRNA-Sequence Profiling Reveals Novel Osmoregulatory MicroRNA Expression Patterns in Catadromous Eel Anguilla marmorata[END_REF] .

In previous studies, the responses of wild immature yellow eels to pollution was investigated in 8 locations presenting a broad contamination gradient in France and Canada [START_REF] Baillon | Transcriptome profile analysis reveals specific signatures of pollutants in Atlantic eels[END_REF][START_REF] Laporte | RAD sequencing reveals within-generation polygenic selection in response to anthropogenic organic and metal contamination in North Atlantic Eels[END_REF][START_REF] Pannetier | A comparison of metal concentrations in the tissues of yellow American eel (Anguilla rostrata) and European eel (Anguilla anguilla)[END_REF][START_REF] Pierron | Retrotransposon methylation and activity in wild fish (A. anguilla): A matter of size[END_REF] . These studies identified Arcachon Bay and the Gironde Estuary as the cleanest and the most contaminated French sites, respectively. We used samples from these two highly contrasted sites in order i) to identify evolutionarily conserved as well as novel miRNAs in the European eel Anguilla anguilla by using next generation sequencing, ii) to compare their expression between a pristine and a polluted site and iii) to link miRNA and mRNA levels. Comprehension of gene regulation by miRNAs may give new insights into a neglected mechanism of response to environmental pollution in aquatic organisms.

Material and Methods

Sample collection -European eels Anguilla anguilla (Linnaeus, 1758) were collected from two sampling sites with a contrasted contamination profile. The Gironde estuary (GE, 45° 12.110'N, 0° 43.579'W) is considered as a highly impacted site while Arcachon Bay (AB, 44° 41.300'N, 1° 1.650'W) represents a pristine environment. Six individuals were collected on each site. The liver was dissected, transferred in RNA-later solution and placed at 4°C overnight before storage at -20°C until RNA extraction. These samples were part of a previous study carried out by our group [START_REF] Baillon | Transcriptome profile analysis reveals specific signatures of pollutants in Atlantic eels[END_REF] .

Condition indices, tissue composition and trace metal analyses -Information concerning the 12

individuals used in the present study (biometry, tissue composition and trace metals) were extracted from previously published data [START_REF] Baillon | Transcriptome profile analysis reveals specific signatures of pollutants in Atlantic eels[END_REF][START_REF] Pannetier | A comparison of metal concentrations in the tissues of yellow American eel (Anguilla rostrata) and European eel (Anguilla anguilla)[END_REF][START_REF] Patey | Biotransformation, antioxidant and histopathological biomarker responses to contaminants in European and American yellow eels from the Gironde and St. Lawrence estuaries[END_REF][START_REF] Caron | Organic and inorganic contamination impacts on metabolic capacities in American and European yellow eels[END_REF] . Fulton index and HSI were calculated as (total weight / length 3 ) x 100 and as (liver weight / total weight) x 100, respectively (weight is in grams and length is in centimeters). The silvering status of the eels was assessed by the calculation of the ocular index (OI) according to the Pankhurst method [START_REF] Pankhurst | Relation of visual changes to the onset of sexual maturation in the European eel Anguilla anguilla (L.)[END_REF] . Statistical differences between the two sites were assessed either with Student t tests or Wilcoxon tests after checking assumptions of normality and homoscedasticity.

Small RNA isolation and deep sequencing -Total RNA were isolated from individual liver tissues following the Trizol™ Reagent (Invitrogen™) protocol. MicroRNA were then purified by using the PureLink™ MiRNA Isolation kit (Invitrogen™) according to the manufacturer's instructions. Total RNA and MicroRNA concentration and quality were assessed by measuring the absorbance at 260nm and 280nm, and migration on a TBE-urea 12% polyacrylamide gel. The NEBNext® Small RNA Library Prep Set for Illumina® (New Englan BioLabs®) was used to prepare the sequencing libraries.

Fragment size (around 150 bp) and quality of the individual libraries were checked on a 2100 Bioanalyzer® (Agilent ®) with a High Sensitivity chip. The 12 individually tagged libraries were then sequenced using a HiSeq 4000 sequencer at the GenomeEast microarray and sequencing platform (Illkirch, France).

MicroRNA discovery analysis -The 50bp raw reads from each of the 12 samples were quality filtered (Q > 30) and adaptors (5' AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC 3') were removed. The Anguilla anguilla genome 37 (GenBank accession number: AZBK00000000.1) was used as a reference genome. High quality, adaptor processed reads were aligned against the eel reference genome using the miRDeep2 package [START_REF] Friedländer | Discovering microRNAs from deep sequencing data using miRDeep[END_REF][START_REF] Friedländer | miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades[END_REF] with default parameters. Our strategy was similar to the one used for the discovery and characterization of miRNA in the salmon Salmo salar [START_REF] Andreassen | Discovery and characterization of miRNA genes in Atlantic salmon (Salmo salar) by use of a deep sequencing approach[END_REF] . Shortly, miRDeep2 provides a list of putative miRNA precursors with their corresponding mature and star sequences. Each precursor is associated with a score that is a measurement of the posterior probability that a predicted sequence is a true miRNA gene. Each score corresponds to a signal-to-noise ratio that estimates total miRNA reported/mean estimated false positive miRNA over 100 rounds of permutated controls. We used a miRDeep2 score of 1 that yielded a signal-to-noise ratio of 10:1 as a cut-off threshold. MiRNAs over this score were further compared by BLAST searches to all known stem-loop sequences from all species deposited in miRBase, release 21 [START_REF] Kozomara | miRBase: annotating high confidence microRNAs using deep sequencing data[END_REF] . Only miRNA with reads matching both harms (5p and 3p) of their precursors were considered. MicroRNA with an E-value < 1.e-6 were annotated as a true Anguilla anguilla miRNA and named accordingly to its ortholog in miRBase.

Sequences with no significant match in miRBase were used as queries in a BLAST analysis against the nr/nt and refseqRNA databases in GenBank 48 , the functional small RNA database [START_REF] Mituyama | The Functional RNA Database 3.0: databases to support mining and annotation of functional RNAs[END_REF] and the small RNA families in Rfam 50 to exclude other kind of small RNA. Finally, the remaining miRNAs were aligned with the reference genome. Sequences with an E-value < 1e-6 against >5 loci were considered as repeats and discarded. Final sequences were then considered as novel miRNAs. The nomenclature rules used for the miRNA naming were in agreement with guidelines from miRBase [START_REF]A uniform system for microRNA annotation[END_REF][START_REF] Griffiths-Jones | miRBase: microRNA sequences, targets and gene nomenclature[END_REF] .

Transcriptome annotation -The transcriptome assembly was previously annotated based on similarity (E value ≤ 1e-10) with known proteins from the swissprot and nr protein databases using the BLASTx algorithm [START_REF] Baillon | Transcriptome profile analysis reveals specific signatures of pollutants in Atlantic eels[END_REF] . For the present work, we enriched this annotation with Gene Ontology (GO) assignment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis that were performed using the BLAST2GO software [START_REF] Conesa | Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research[END_REF][START_REF] Götz | High-throughput functional annotation and data mining with the Blast2GO suite[END_REF] and the GhostKOALA server [START_REF] Kanehisa | BlastKOALA and GhostKOALA: KEGG Tools for Functional Characterization of Genome and Metagenome Sequences[END_REF] , respectively.

Enrichment tests were performed with the Fisher exact test procedure and a correction for multiple testing using the Benjamini and Hochberg method. A Pvalue threshold of 0.01 was used. Redundancy in GO terms was reduced by the GO Trimming method 56 using a uniqueness threshold of 50%. For a matter of clarity, only GO and KEGG categories containing ≥5 sequences were considered.

Differential expression of transcripts -Differences in transcription levels between the two study

sites were assessed with the RNAseq data extracted from Baillon et al [START_REF] Baillon | Transcriptome profile analysis reveals specific signatures of pollutants in Atlantic eels[END_REF] for the exact same individual eels. Normalization and expression analysis were performed using the DESeq2 package [START_REF] Love | Moderated estimation of fold change and dispersion for RNAseq data with DESeq2[END_REF] . P-values for differential gene expression were corrected for multiple testing using the Benjamini and Hochberg method, and a false discovery rate (FDR) threshold of 0.05 was used.

Expression analysis of Eel miRNAs -

The copy number of known and novel mature miRNAs were counted with the quantifier module in miRDeep2 (additional files 1 and 3). Normalization and expression analysis were performed using the DESeq2 package [START_REF] Love | Moderated estimation of fold change and dispersion for RNAseq data with DESeq2[END_REF] . The false discovery rate (FDR) threshold was 0.05. In order to validate these results, the expression profiles of 10 mature miRNAs were assessed using qRT-PCR (additional file 4). Reverse transcriptions were performed by using the miScript II RT system (Qiagen®). Real Time PCR Analyses were performed on a Roche LightCycler ® 480 Instrument II by using the miScript PCR Starter Kit (Qiagen®). Cycling conditions were 95°C for 15 seconds followed by 40 cycles of 15 seconds at 94°C, 30 seconds at 55°C and 30 seconds at 70°C.

Target prediction -The sequences and annotation of protein-coding transcripts were retrieved from

Baillon et al [START_REF] Baillon | Transcriptome profile analysis reveals specific signatures of pollutants in Atlantic eels[END_REF] . The target genes of mature eel miRNAs were predicted using the miRanda software [START_REF] Enright | MicroRNA targets in Drosophila[END_REF] that is based on the local alignment between the query miRNA sequences and the reference transcriptome and on the thermodynamic stability of miRNA/mRNA duplexes (parameters used were score ≥ 160, free energy ≤ -25kcal/mol).

Results

Trace metals analysis -Differences between individuals collected in both sites are presented in Table 1. First, AB and GE sites differed in salinity and dissolved oxygen concentrations. We observed no difference in Fulton condition factor between the two sites. However, the Hepato-Somatic Index (HSI) of AB eels was higher than GE eels. No difference was detected in OI and with a value below the threshold value of 6.5 [START_REF] Pankhurst | Relation of visual changes to the onset of sexual maturation in the European eel Anguilla anguilla (L.)[END_REF] , all the individuals were indeed at the yellow eels stage. As expected from previous results [START_REF] Pannetier | A comparison of metal concentrations in the tissues of yellow American eel (Anguilla rostrata) and European eel (Anguilla anguilla)[END_REF] , trace metal concentrations were significantly higher in the liver of eels from GE, with the exceptions of Cr and Zn that showed no difference and Ni that was higher in AB. Cadmium showed the largest difference with a concentration more than 200-times higher in the liver of GE individuals in comparison to fish from AB.

Small RNA sequencing -The sequencing of the 12 sRNA libraries led to a total of 362,831,354 raw reads out of which 318,968,039 (87.91%) were kept for the following analyses after quality filtering (Q > 30) and adaptor removal (detailed statistics are given in Table 2). With the exception of the sample AB03 (59.35%), between 72.12% (AB06) and 89.40% (GE02) of these reads were 21-23 nucleotides (nt) in length (Figure 1). Another population of 28-30 nt was identified, representing on average 12.99% and 2.36% of the reads in AB and GE samples, respectively. These two short reads populations represented between 81.70% and 90.76% of the total high-quality reads. Up to 87.73% of these reads could be mapped against the reference genome for miRNA prediction (Table 2). The raw sequencing data from this study were submitted to the NCBI Gene Expression Omnibus 59 under GEO Accession GSE109689.

Identification of known and novel miRNAs in A.

anguilla -A total of 820 putative miRNA precursors were predicted by the miRDeep2 analysis (additional file 1) amongst which 474 had a score above 1 and reads on both harms, i.e. were associated with both a 3p and a 5p sequence (Figure 2). Amongst those, 230 pre-miRNA provided a match against a stem-loop sequence in miRBase (E-value < 1e-6) and were considered Anguilla anguilla orthologs of an evolutionary conserved miRNA gene and named in accordance with the nomenclature rules. Most of these conserved miRNAs (181 out of 230) matched with miRNAs from other fish species in miRBase. They corresponded to 158 unique mature miRNAs from 96 evolutionary conserved families. The remaining 244 putative precursors were further analyzed by BLAST searches against the reference genome and other RNA databases.

That allowed us to identify 145 potential novel pre-miRNAs that were named based on their genome location. Redundancy was lower since 141 unique novel mature miRNAs were identified. Taken together we thus identified a total of 375 miRNA precursors with their respective mature and star sequences. The total number of unique mature miRNAs was 299 (detailed list provided in additional file 2). In all libraries, aan-mir-122 was on average 14 times more expressed than the second most abundant mature miRNA (ranging from ~4 in AB06 to ~25 times in GE05). Other abundant miRNAs included aan-mir-21, aan-mir-22 and aan-let-7a (additional file 3).

Prediction of transcriptomic targets -

In order to determine the biological processes potentially regulated by the 299 mature microRNAs we identified, we predicted their target genes in the anguilla's transcriptome assembly [START_REF] Baillon | Transcriptome profile analysis reveals specific signatures of pollutants in Atlantic eels[END_REF] . The A. anguilla transcriptome contains 18,113 contigs with high homology with known sequences (BLASTx, Evalue ≤ 1e-10). Although 3'UTR regions are the most important regulatory sites for post-transcriptional regulation of gene expression by miRNAs, we considered whole length transcripts. 242 miRNAs showed a potential to hybridize within the transcriptome wherein 3,637 contigs were identified as potential targets (additional file 4). As expected, one miRNA could target several transcripts and one transcript could be targeted by several miRNAs. For instance, the novel miRNA aan-scaffold-2762 could hybridize with up to 271 targets and the transcript contig_15454, encoding a homologue of E3 ubiquitin-protein ligase UBR2 (Uniprot accession number Q8IWV8), was potentially the target of 8 different miRNAs. The entire set of predicted mRNA targets was further functionally categorized through GO annotation and KEGG pathway analysis. No GO term, nor KEGG pathway, was significantly enriched in the list of target genes compared to the transcriptome background, suggesting that miRNAs might affect any biological process.

Differentially expressed mature miRNAs between the two sites -One objective of this study was to evaluate the response of miRNAs to differing levels of contamination. A multidimensional scaling of the count data showed that despite inter-individual variability, samples are well separated according to their site of origin (additional file 5). We identified 41 mature miRNAs with differential expression (DE) in GE compared to the reference site (AB). Nineteen (19) miRNAs were up-regulated, while 22 were down-regulated (Table 3). The most up-regulated and down-regulated miRNAs were respectively mir-29a (2.77 times more expressed in GE) and mir-462 (3.97 times less expressed in GE). These results were validated by qRT-PCR (r 2 = 0.88, additional file 5). Amongst all the DE miRNAs, 35 were likely to target 490 sequences in A. Anguilla's transcriptome (Table 3). The targets of up-regulated miRNAs were functionally enriched in transcripts encoding proteins involved in ribosome biogenesis. This process is hence expected to decrease in GE. Biological processes affected by down-regulated miRNAs were related to the response to hormones (e.g. GO:0016568), response to chemicals (e.g. GO:0071417) and chromatin organisation (e.g. GO:0006325) (Figure 3). These processes are then expected to increase in GE. [START_REF] Baillon | Transcriptome profile analysis reveals specific signatures of pollutants in Atlantic eels[END_REF] , we identified 1,554 down-regulated and 1,343 up-regulated genes in the Gironde Estuary compared to the Arcachon Bay (additional file 7). Down-regulated genes contained a significant proportion of sequences related to lipid metabolism (GO:0006629 and GO:0044255), response to chemical (GO:0001101 and GO:0009410), primary metabolism, development (GO:0051216 and GO:0032836) and hormone metabolism (GO:0042445). Up-regulated transcripts were involved in processes such as transport (GO:0006810, GO:0051649 and GO:0032940), development (GO:0007423, GO:0007389, GO:0048840 and GO:0035889) and response to stimulus (GO:0010035 and GO:0071216) (Figure 3). 4).

Links between DE miRNAs and DE transcripts -Based on previous data

Discussion

The increasing evidence that the expression of miRNAs is affected by several known toxicants as well as oxidative and other forms of cellular stress certainly suggests an important role of miRNAs in toxicology, which could provide a link between environmental influences and gene expression [START_REF] Lema | MicroRNAs and their implications in toxicological research[END_REF] . The use of miRNAs in environmental toxicology is expanding, particularly in human health [START_REF] Huang | MicroRNA Machinery Genes as Novel Biomarkers for Cancer[END_REF][START_REF] Yu | The role of microRNAs in toxicology[END_REF][START_REF] Hou | Environmental chemicals and microRNAs[END_REF] . This interest mainly results from the strong potential of miRNAs to become robust, minimally invasive and non-expensive biomarkers for various diseases [START_REF] Schraml | MicroRNAs and toxicology: A love marriage[END_REF] . In the same manner, examples of miRNA-related work are more frequent in aquatic organisms exposed to environmental stressors. One can cite some recent publication in bivalves [START_REF] Zhou | The identification and characteristics of immune-related microRNAs in haemocytes of oyster Crassostrea gigas[END_REF][START_REF] Bao | Identification and comparative analysis of the Tegillarca granosa haemocytes microRNA transcriptome in response to Cd using a deep sequencing approach[END_REF][START_REF] Xu | Identification of conserved and novel microRNAs in the Pacific oyster Crassostrea gigas by deep sequencing[END_REF] , crustaceans [START_REF] Chen | Novel Cadmium Responsive MicroRNAs in Daphnia pulex[END_REF][START_REF] Chen | MicroRNAs are involved in cadmium tolerance in Daphnia pulex[END_REF] or teleost fish [START_REF] Wang | MicroRNA-Sequence Profiling Reveals Novel Osmoregulatory MicroRNA Expression Patterns in Catadromous Eel Anguilla marmorata[END_REF][START_REF] Zhang | Integrated analysis of mRNA-seq and miRNA-seq in the liver of Pelteobagrus vachelli in response to hypoxia[END_REF][START_REF] Lee | Exposure to bisphenol S alters the expression of microRNA in male zebrafish[END_REF][START_REF] Best | Epigenetics in teleost fish: From molecular mechanisms to physiological phenotypes[END_REF] for instance. To our knowledge, our work is the first in situ study using miRNA in an aquatic organism. MicroRNAs are being characterized in an increasing number of teleost fish species where they repress target mRNA abundance or translation in a highly conserved manner across most eukaryotes [START_REF] Bizuayehu | MicroRNA in teleost fish[END_REF] .

In this study, we conducted a comprehensive identification of the microRNA repertoire in the liver of the European eel A. anguilla and its response to environmental pollution. This organ was chosen owing to its role in pollutant metabolism, accumulation and toxicity. Only trace metal concentrations were reported here since they were measured in the liver. However, one should note that eels from the Gironde Estuary also showed higher levels of persistent organic pollutants (Poly-Chlorinated Biphenyls, Organo-Chlorine Pesticides, and Poly-Brominated Diphenyl Ethers) in their muscle compared to eels from the Arcachon Bay, as evidenced by previous analyses [START_REF] Baillon | Transcriptome profile analysis reveals specific signatures of pollutants in Atlantic eels[END_REF][START_REF] Patey | Biotransformation, antioxidant and histopathological biomarker responses to contaminants in European and American yellow eels from the Gironde and St. Lawrence estuaries[END_REF] . Owing to panmixia and proximity between the two sites, this contrast is most likely due to water quality rather than population differences.

The size distribution of the small RNA libraries revealed that 22-23 nt microRNAs were dominant, followed by small RNAs of 28-29 nt in length (Figure 1). This distribution is similar to what was found in the sister species A. marmorata [START_REF] Wang | MicroRNA-Sequence Profiling Reveals Novel Osmoregulatory MicroRNA Expression Patterns in Catadromous Eel Anguilla marmorata[END_REF] . It is interesting to note that, in both eel species, the relative abundance of 22-23 nt and 28-29 nt varies depending on environmental conditions, which could represent a new path of investigation in the future. These 28-29 nt sRNA might correspond to piRNA (piwi-interacting RNA that have a length of 24-32 nt 1,69 ), another smallRNA silencing system that plays a role in germline specification, gametogenesis, transposon silencing, genome integrity, and stem cell maintenance across the animal phylogeny [START_REF] Iwasaki | PIWI-Interacting RNA: Its Biogenesis and Functions[END_REF] .

We identified a total of 299 mature miRNA in eel liver, of which 158 were known and 141 were novel. This number is comparable to what is currently available in the miRBase database for the model species Danio rerio (375 mature sequences). This identification benefited from the availability of a sequenced genome of A. anguilla [START_REF] Henkel | Primitive Duplicate Hox Clusters in the European Eel's Genome[END_REF] . Only mature miRNAs were considered in our analysis since they eventually regulate gene expression at the post-transcriptional level. In all libraries, mir-122 was the most abundant mature sequence. This was not surprising since it is the most abundant liver miRNA in vertebrates (>70% of the total hepatic miRNA expression) with exquisite tissue specificity.

Mir-122 is mainly involved in lipid homeostasis [START_REF] Willeit | Liver microRNAs: potential mediators and biomarkers for metabolic and cardiovascular disease?[END_REF] . Its expression level increased by 2.34-fold in polluted water (Table 3). Such accumulation is associated to injured hepatocytes and is a biomarker of liver toxicity in zebrafish [START_REF] Vliegenthart | Retro-orbital blood acquisition facilitates circulating microRNA measurement in zebrafish with paracetamol hepatotoxicity[END_REF] and mammals [START_REF] Vliegenthart | Comprehensive microRNA profiling in acetaminophen toxicity identifies novel circulating biomarkers for human liver and kidney injury[END_REF] . The alteration of lipid metabolism was also evidenced by prevalence of related Gene Ontology terms amongst down-regulated transcripts (Figure 3). In farmed tilapia (GIFT, Oreochromis niloticus), the level of mir-122 was negatively correlated to Cd stress and promotes hepatic antioxidant defence [START_REF] Qiang | miR-122 promotes hepatic antioxidant defense of genetically improved farmed tilapia (GIFT, Oreochromis niloticus) exposed to cadmium by directly targeting a metallothionein gene[END_REF] . This discrepancy with the present work, where Cd was the metal showing the largest difference between the two sites, could either be due to the action of other stressors or to taxa specific mechanisms of miRNA regulation. The two other abundant miRNAs, mir-21 and mir-22, were also shown to play a role in the development of hepatic diseases [START_REF] Yang | MiR-22-silenced Cyclin A Expression in Colon and Liver Cancer Cells Is Regulated by Bile Acid Receptor[END_REF][START_REF] Caviglia | MicroRNA-21 and Dicer are dispensable for hepatic stellate cell activation and the development of liver fibrosis[END_REF] . Liver damage was confirmed by a series of transcriptomic changes described previously [START_REF] Baillon | Transcriptome profile analysis reveals specific signatures of pollutants in Atlantic eels[END_REF] . This first set of results shows the high tissue specificity of our approach and confirms the potential of miRNAs in aquatic ecotoxicology.

We chose to study two contrasted sites both in terms of organic and metallic contaminants [START_REF] Baillon | Transcriptome profile analysis reveals specific signatures of pollutants in Atlantic eels[END_REF][START_REF] Laporte | RAD sequencing reveals within-generation polygenic selection in response to anthropogenic organic and metal contamination in North Atlantic Eels[END_REF] . These sites also differed in salinity (19.9‰ and 7.78‰ for Arcachon Bay and Gironde Estuary, respectively) and oxygen percentage (60% vs. 86.5%) whereas the temperature was relatively similar (21.9°C vs. 1). With a dissolved oxygen concentration of 4.68 mg/l, site AB faces moderate hypoxia according to the European Water Framework Directive (WFD, 2000/60/EC). In such a complex context, it is obviously challenging to identify causing agents. However, since contamination of aquatic ecosystems typically includes many pollutants that present different mechanisms of action likely interacting together and under the influence of natural factors, this situation might provide more realistic information as well as a first insight into the role of miRNAs in a polluted environment. In fish, salinity and hypoxia were previously shown to affect the miRNA regulatory machinery [START_REF] Wang | MicroRNA-Sequence Profiling Reveals Novel Osmoregulatory MicroRNA Expression Patterns in Catadromous Eel Anguilla marmorata[END_REF][START_REF] Zhang | Integrated analysis of mRNA-seq and miRNA-seq in the liver of Pelteobagrus vachelli in response to hypoxia[END_REF][START_REF] Sun | Dynamic mRNA and miRNA expression analysis in response to hypoxia and reoxygenation in the blunt snout bream ( Megalobrama amblycephala )[END_REF] . In our study, only 5 miRNAs (miR-19b, miR-139, miR99a, miR-454 and miR-221) and 3 miRNAs (miR-125a, miR-181c and miR-30b) were found to respond to salinity and hypoxia, respectively, in the same way than in the literature cited above. Moreover, eels were shown to be able to cope with much more severe hypoxia, down to about 30% O2 saturation [START_REF] Cruz-Neto | The effects of acute hypoxia and hypercapnia on oxygen consumption of the freshwater European eel[END_REF] .Then, the differential expression of the majority of the miRNA we identified could most likely be linked to the presence of pollutants. The expression level of 19 miRNAs increased in the polluted site while the expression of 22 miRNAs was higher in AB. This relatively small number of DE miRNA is similar to what was found in other fish species exposed to environmental stressors [START_REF] Wang | Copper-induced deregulation of microRNA expression in the zebrafish olfactory system[END_REF][START_REF] Wang | MicroRNA-Sequence Profiling Reveals Novel Osmoregulatory MicroRNA Expression Patterns in Catadromous Eel Anguilla marmorata[END_REF][START_REF] Zhang | Integrated analysis of mRNA-seq and miRNA-seq in the liver of Pelteobagrus vachelli in response to hypoxia[END_REF][START_REF] Kure | Molecular responses to toxicological stressors: profiling microRNAs in wild Atlantic salmon (Salmo salar) exposed to acidic aluminum-rich water[END_REF] and could be due to the relatively low number of samples per site and inter-individual variability (additional file 5) that should be taken into account in future studies. The most expressed known miRNAs in AB (>3 fold) were miR-462, miR-15e and miR-23a, and the most up-regulated known miRNAs in GE (>2 fold) were miR-29a, miR-122 (discussed above) and miR-203. As previously, the fold change values measured here are in agreement with what can be found for Teleosts in the literature [START_REF] Wang | Copper-induced deregulation of microRNA expression in the zebrafish olfactory system[END_REF][START_REF] Wang | MicroRNA-Sequence Profiling Reveals Novel Osmoregulatory MicroRNA Expression Patterns in Catadromous Eel Anguilla marmorata[END_REF][START_REF] Zhang | Integrated analysis of mRNA-seq and miRNA-seq in the liver of Pelteobagrus vachelli in response to hypoxia[END_REF][START_REF] Kure | Molecular responses to toxicological stressors: profiling microRNAs in wild Atlantic salmon (Salmo salar) exposed to acidic aluminum-rich water[END_REF] . The role of these miRNAs could be inferred from previous studies.

21.3°C) (Table

Despite the highly conserved mode of action of miRNAs, caution must be taken when comparing fish and mammals as only a small fraction of miRNA-mRNA interactions (~10%) is conserved [START_REF] Xu | The evolution of evolvability in microRNA target sites in vertebrates[END_REF] and identification of conserved targets between species remains challenging 80 . For instance, in fish, miR-462 was shown to be upregulated after hypoxia treatment in liver [START_REF] Huang | The zebrafish miR-462/miR-731 cluster is induced under hypoxic stress via hypoxia-inducible factor 1α and functions in cellular adaptations[END_REF][START_REF] Huang | Involvement of the miR-462/731 cluster in hypoxia response in Megalobrama amblycephala[END_REF] or virus infection [START_REF] Andreassen | Identification of differentially expressed Atlantic salmon miRNAs responding to salmonid alphavirus (SAV) infection[END_REF][START_REF] Schyth | Two Virus-Induced MicroRNAs Known Only from Teleost Fishes Are Orthologues of MicroRNAs Involved in Cell Cycle Control in Humans[END_REF] . MiR-462 is only described in teleost fish and is actually an orthologue of the human miR-191 that plays a role in cell cycle regulation, but not in viral infection [START_REF] Schyth | Two Virus-Induced MicroRNAs Known Only from Teleost Fishes Are Orthologues of MicroRNAs Involved in Cell Cycle Control in Humans[END_REF] . The role of miR-23a instead seems to be conserved.

It was identified as an essential regulator of adipocyte differentiation of which increased expression might reduce lipid accumulation and triglyceride content in adipocytes [START_REF] Shen | MicroRNA-23a regulates 3T3-L1 adipocyte differentiation[END_REF] . This role was confirmed in our study by the identification of a transcript encoding a very long-chain acyl-CoA synthetase (Uniprot accession O14975, additional file 5) as a putative target of miR-23a. In human, miR-29a is known to regulate cell proliferation, differentiation and apoptosis as several of its targets are oncogenes or anti-apoptotic genes [START_REF] Kriegel | The miR-29 family: genomics, cell biology, and relevance to renal and cardiovascular injury[END_REF] . Amongst the predicted mRNA targets of miR-29a, we identified a thyroid receptor-interacting protein homologue (Uniprot accession Q15643, additional file 5) of which degradation could reduce gene transcription in response to thyroid hormone while miR-122 predicted targets included transcripts encoding proteins involved in lipid synthesis, energy generation, platelet activation, apoptosis and transcription regulation. The role of miR-203 remains unclear as its overexpression might promote both apoptosis and tumor growth, as well as host-virus interactions [START_REF] Deng | MiRNA-203 suppresses cell proliferation, migration and invasion in colorectal cancer via targeting of EIF5A2[END_REF] . This difference in lipid metabolism is not confirmed by our observation as fish from the two sites showed no significant difference regarding their muscular lipid content (table 1), contrary to previous studies that showed that eels exposed to Cd 87 and organic pollutants [START_REF] Sancho | Effects of sublethal exposure to a pesticide on levels of energetic compounds in Anguilla anguilla[END_REF] increased their fat consumption and/or reduced their energy reserves. We can here hypothesize that miRNA regulation in the liver contributes to maintain a significant production of triglycerides and their export to the main site of lipid storage in eels, i.e. the muscle [START_REF] Pierron | Impairment of lipid storage by cadmium in the European eel (Anguilla anguilla)[END_REF] . At the opposite, the lower HSI values in fish from the estuary could indicate an impairment in the accumulation of energy reserves under the form of carbohydrate such as glycogen in hepatocytes.

More comprehensive information was gained from the functional analysis of miRNA targets in A.

Anguilla's transcriptome (Figure 3). Based on gene ontology of these targets, changes in miRNA expression might result in reduced ribosome biogenesis and increased response to hormone and chemical, and chromatin modification. The effect on ribosome biogenesis could be evidenced by a significantly reduced level of total proteins in the liver of GE eels (Table 1). One could expect to find some similar functions in the predicted targets of up-regulated miRNAs and down-regulated transcripts, and conversely. In spite, we found no relation between enriched gene ontology in miRNA targets and mRNA, as well as a small number of miRNA-mRNA pairs with negative correlation (Table 4). This situation was also observed in a comparable study carried out in Pelteobagrus vachelli exposed to hypoxia [START_REF] Zhang | Integrated analysis of mRNA-seq and miRNA-seq in the liver of Pelteobagrus vachelli in response to hypoxia[END_REF] in which the authors identified only 97 negative correlations out of 308 miRNA-mRNA pairs. In the same study, and in the same manner as our results, most of the DE mRNA (43 out of 60) were target by a single DE miRNA. This highlights that the identification of functional miRNAtargets interaction is still challenging and/or false positives might still be frequent with in silico analyses [START_REF] Tabas-Madrid | Improving miRNA-mRNA interaction predictions[END_REF] . MicroRNAs are obviously not the only mode of regulation of gene expression, however based on the lack of correspondence between miRNA targets and DE transcripts, we could hypothesize that the degradation of messenger RNA was not the major mode of action of miRNAs here and that they act more likely by inhibiting mRNA translation into proteins, instead. This inhibition might be reinforced by the targeting of mRNA involved in ribosome assembly and function.

This result is in accordance with those obtained in juvenile zebrafish exposed to the antibacterial agent triclosan where "structural constituent of ribosome" and "ribosome biogenesis" are amongst the significantly enriched annotation of differentially expressed miRNAs target genes [START_REF] Wang | Exposure to triclosan changes the expression of microRNA in male juvenile zebrafish (Danio rerio)[END_REF] . Which might suggest a conserved mode of action of miRNA in response to pollutants in fish. Another way by which miRNA could regulate gene expression is through chromatin modification, as suggested by the presence of 28 epigenetic regulators in miRNA targets (additional file 8). In human, miRNA and epigenetic regulation are not entirely separable and a strong interplay exists, but the mechanisms of cooperation are still unclear [START_REF] Osella | Interplay of microRNA and epigenetic regulation in the human regulatory network[END_REF][START_REF] Wang | Coordinated action of histone modification and microRNA regulations in human genome[END_REF] . Compared to mammals, the study of epigenetics in teleost fish has received less attention so far [START_REF] Best | Epigenetics in teleost fish: From molecular mechanisms to physiological phenotypes[END_REF] . Pollutant-induced changes in epigenetic marks were already documented [START_REF] Pierron | Retrotransposon methylation and activity in wild fish (A. anguilla): A matter of size[END_REF][START_REF] Head | Epigenetics for ecotoxicologists[END_REF][START_REF] Vandegehuchte | Epigenetics and its implications for ecotoxicology[END_REF] . For instance, transcriptomic changes in the liver of wild yellow perch (Perca flavescens) exposed to cadmium and copper may be driven through histone modifications [START_REF] Pierron | Effects of chronic metal exposure on wild fish populations revealed by highthroughput cDNA sequencing[END_REF] . The present work suggests that similar effects in eels could be mediated by miRNA regulation.

The European eel spends most of its lifetime (about 8-15 years) in European estuarine and freshwater basins in a stage known as yellow eels. During this period, eels accumulate a considerable amount of lipid reserves, which are fundamental for the success of the long oceanic migration since during this period they do not feed. Before initiating their oceanic migration to reproduce in the Sargasso Sea, they metamorphose into silver eel [START_REF] Feunteun | Management and restoration of European eel population (Anguilla anguilla): An impossible bargain[END_REF][START_REF] Guimarães | Yellow eel (Anguilla anguilla) development in NW Portuguese estuaries with different contamination levels[END_REF] . Laboratory and field studies have shown that eels can bioaccumulate high concentrations of environmental contaminants, which have been shown to cause mortality, growth delay, reproductive alterations, tumours, malformations, nervous and endocrine disturbances and immunological changes [START_REF] Geeraerts | The effects of contaminants in European eel: a review[END_REF] . In the present work, we evidenced a role of miRNAs in the regulation of genes involved in all these processes, reinforcing the potential of miRNAs as biomarkers of yellow eel's health status in a polluted environment. MiRNAs involved in lipid metabolism (like miR-122) and endocrine function might be of particular interest owing to the critical importance of these biological processes in the reproductive success of A. anguilla. Although lipid storage might not be absolutely necessary to initiate metamorphosis (silvering [START_REF] Svedäng | Low fat contents in female silver eels: indications of insufficient energetic stores for migration and gonadal development[END_REF] ), a sufficient energy reserve is needed to cover the migratory needs, ensure gonadal maturation and reproduction.

Moreover, morphological and physiological changes that take place during the silvering phase are initiated and regulated by hormonal changes [START_REF] Geeraerts | The effects of contaminants in European eel: a review[END_REF] .

The present study was carried out on wild organisms, in a multi-stress context. In such conditions, it is difficult to establish significant relationships between the expression level of one particular miRNA and one particular contaminant. Further experiments carried out with more individuals, in controlled conditions, different organs and/or at different developmental stages are thus needed to fully reveal the potential of miRNAs as biomarkers.

Conclusion

To the best of our knowledge, this study is the first example of high-throughput sequencing of microRNAs in wild aquatic organisms. We were able to identify both conserved and putative novel miRNAs and their targets at the transcriptomic level in the liver of European eels from two sites with contrasted pollution profiles. We showed that some miRNA were differentially expressed in organisms experiencing chronic exposure to pollution and that these changes have the potential to affect important biological processes such as lipid and hormone metabolism or protein synthesis. Moreover, our work also revealed promising aspects of epigenetics in teleost fish, as miRNAs might regulate gene expression by interacting with other epigenetic mechanisms such as chromatin modification. Part of the biological processes affected by miRNAs for which expression levels were altered were similar to previous transcriptome studies, confirming the potential role of miRNAs as early molecular markers for monitoring stress responses in exposed animals.

Arcachon

  Due to the action mode of miRNAs, their expression patterns generally showed an inverse relationship with those of their mRNA targets. We compared the list of DE genes with the lists of DE miRNAs and their targets. Out of the 35 DE miRNAs that had putative targets in the transcriptome, 20 of them were able to hybridize to 49 (out of 490) DE mRNAs targets, forming 52 miRNA-mRNA pairs with both positive and negative relationships in their expression values (log2FC). In detail, there were 29 (out of 52) negative miRNA-mRNA interactions (Spearman's ρ = -0.69, Pvalue = 3.46e-05) involving 18 DE miRNAs and 29 DE mRNAs. Whilst a single DE miRNA could target up to 4 different DE transcripts, each DE transcript was targeted by only one DE miRNA. Many of these DE transcripts were involved in biological processes cited above, such as immune response and lipid metabolism (Table

Figure 1 :

 1 Figure 1: Length distribution of high-quality reads in Illumina libraries. Left panel: Arcachon Bay (AB), Right panel: Gironde Estuary (GE). 663

Figure 2 :

 2 Figure 2: Pie chart showing the distribution of the 820 miRNA sequences predicted by miRDeep2. N=number of sequences in each category.

Figure 3 :

 3 Figure 3: GO enrichment analysis (p-value ≤ 0.01; ≥ 5 sequences) of the transcript targets of differentially expressed miRNAs (Left panel) and differentially expressed mRNAs (Right panel), A: targets of up-regulated miRNAs. B: down-regulated mRNAs, C: targets of down-regulated miRNAs, D: up-regulated mRNAs. X-axis, enrichment score (-log10 Pvalue); Y-axis, Biological processes. The number of sequences in each category is presented on the right of the corresponding bar.

Table 1 .

 1 Environmental variables, condition indices, tissue composition, and metal concentrations in the liver (mean ± SE, n=6) of European eels sampled in Arcachon Bay and Gironde Estuary. HSI: hepato-somatic index, OI: ocular index, ww: wet weight, dw: dry weight. (*, ** and *** symbols are Pvalue ≤0.05, ≤0.01 and ≤0.001, respectively)

	Bay	Gironde Estuary

Table 2 :

 2 Summary of samples sequenced for discovery of Anguilla anguilla miRNA genes. AB, Arcachon Bay; GE, Gironde Estuary.

	653
	654
	655

Table 3 :

 3 Relative expression and number of target transcripts of differentially expressed miRNA 657 (FDR≤0.05) in GE compared to AB groups. 658

Table 4 :

 4 Negative relationship between the relative expression of DE miRNA and their DE targets. Up-regulated and Down-regulated genes (FDR≤0.05) in GE vs. AB are indicated in red and green, respectively.

			contig_14686 0.77	Q9NYF5 Protein FAM13B	2e-135 Signal transduction
			contig_8899	1.42	P60467	Protein transport protein Sec61 subunit beta	5e-21	Secretion system
	aan-scaffold_951	-1.71	contig_16299 1.13	A5PN09 Ubiquitin carboxyl-terminal hydrolase 20	9e-128 Ubiquitin system
			contig_38623 0.76	A2AQ19	RNA polymerase-associated protein RTF1 homolog	6e-167 Transcription machinery

Funding

MicroRNA libraries' quality check was performed at the Genome Transcriptome Facility of Bordeaux (grants from the Conseil Régional d'Aquitaine 20030304002FA and 20040305003FA, from the European Union FEDER 2003227 and from Investissements d'Avenir ANR-10-EQPX-16-01). Sequencing was performed by the IGBMC Microarray and Sequencing platform, a member of the 'France Génomique' consortium (ANR-10-INBS-0009). This work was supported by the Regional Council of Aquitaine, the Agence Nationale de la Recherche of France (ANR Immorteel ANR-10-INTB-1704) and a strategic project grant of the Natural Sciences and Engineering Research Council of Canada to Patrice Couture, jointly funded by the Innovative Project grant of the UMR EPOC 5805.

Appendix S1: miRNAs prediction by miRDeep2. Appendix S2: List of the 375 conserved and novel miRNA in A. Anguilla