domain: sdv info

text:
Taina Rocha De Almeida

Maud Alix

Aurélie Le Cam

Christophe C Klopp

Jérôme Montfort

Lola Toomey

Yannick Ledore

Julien Bobe

Dominique Chardard

Bérénice Schaerlinger
email: berenice.schaerlinger@univ-lorraine.fr

Tainá Rocha De Almeida

Yannick Ledoré

Pascal Fontaine

Domestication may affect the maternal mRNA profile in unfertilized eggs, potentially impacting the embryonic development of Eurasian perch (Perca fluviatilis)

Keywords: microarray, egg content, fish populations, reproductive cycle, gonadogenesis, egg quality

à la diffusion de documents scientifiques de 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

Domestication is an evolutionary and continuous process enabling wild animals to adapt to humans and anthropogenic environments [START_REF] Price | Behavioral Aspects of Animal Domestication[END_REF][START_REF] Price | Behavioral development in animals undergoing domestication[END_REF]. It involves the combination of several genetic and potentially epigenetic modifications driving gene expression and phenotypic changes. Differences in gene expression due to rapid adaptation to new environments were already reported from the first generation of reproduction in captivity in steelhead trout (Oncorhynchus mykiss) [START_REF] Christie | A single generation of domestication heritably alters the expression of hundreds of genes[END_REF]. Thus, the regulation of genes expression may change early during domestication. The main genetic mechanisms involved are inbreeding, genetic drift and natural, relaxed and artificial selections [START_REF] Price | Behavioral Aspects of Animal Domestication[END_REF][START_REF] Price | Behavioral development in animals undergoing domestication[END_REF]. The two first mechanisms have important consequences when founder populations are small because they rapidly lead to important changes in allelic frequencies [START_REF] Beaumont | Génétique et adaptation des animaux d'élevage : introduction[END_REF]. The relaxed selection consists of a reduction of selection pressure on traits which are not necessary anymore in captive conditions. The artificial selection results from the selection of breeders according to phenotypes chosen by humans. Finally, the natural selection occurs and usually eliminates animals which are not adapted to anthropogenic environments [START_REF] Boissy | Genetics of adaptation and domestication in livestock[END_REF]. Studies suggest a relationship between epigenetic modifications and phenotypic plasticity in response to the environment in some animal and plant species [START_REF] Vogt | Facilitation of environmental adaptation and evolution by epigenetic phenotype variation: insights from clonal, invasive, polyploid, and domesticated animals[END_REF]. In fish, the current knowledge shows a modification of the epigenetic signature in individuals reared in hatcheries in comparison to their wild counterparts in salmonids and the European sea bass (Dicentrarchus labrax) [START_REF] Luyer | Parallel epigenetic modifications induced by hatchery rearing in a Pacific salmon[END_REF][START_REF] Gavery | Characterization of Genetic and Epigenetic Variation in Sperm and Red Blood Cells from Adult Hatchery and Natural-Origin Steelhead, Oncorhynchus mykiss. G3 Genes Genomes[END_REF][START_REF] Barreto | DNA Methylation Changes in the Sperm of Captive-Reared Fish: A Route to Epigenetic Introgression in Wild Populations[END_REF][START_REF] Anastasiadi | Epimutations in developmental genes underlie the onset of domestication in farmed European sea bass[END_REF]. These modifications affect somatic and germinal cells [START_REF] Gavery | Characterization of Genetic and Epigenetic Variation in Sperm and Red Blood Cells from Adult Hatchery and Natural-Origin Steelhead, Oncorhynchus mykiss. G3 Genes Genomes[END_REF] and may play a role during the first steps of domestication [START_REF] Anastasiadi | Epimutations in developmental genes underlie the onset of domestication in farmed European sea bass[END_REF]. However, the relationships and timeline between genetic and epigenetic modifications remain unclear [START_REF] Vogt | Facilitation of environmental adaptation and evolution by epigenetic phenotype variation: insights from clonal, invasive, polyploid, and domesticated animals[END_REF].

All these mechanisms depend greatly on the breeding practices. Indeed, independent trials of domestication may lead potentially to various types of modifications that may have either beneficial (adaptation) or deleterious (maladaptation) effects on various biological functions of the new domesticated populations. The artificial selection, for one or several phenotypes is specific to the domestication process and has consequences that may not be predicted, since morphological, behavioral and physiological traits of animals are intrinsically related [START_REF] Dobney | Genetics and animal domestication: new windows on an elusive process[END_REF]. The consequences of modifications of the balance between these traits are not yet properly understood. Indeed, it appears that the artificial selection of specific phenotypes may have deleterious effects on other biological traits because most of resource intakes are dedicated to the selected traits. This imbalance often leads to a decrease of reproduction abilities [START_REF] Boissy | Genetics of adaptation and domestication in livestock[END_REF][START_REF] Waaij | A resource allocation model describing consequences of artificial selection under metabolic stress[END_REF], commonly seen in numerous terrestrial and aquatic species. A recent meta-analysis investigating the effect of birth-origin (captive vs. wild) on the reproductive success of animals reared in different anthropogenic environments was performed. For all of the 44 analyzed species, the offspring survival success was decreased in "captive-born" animals compared to their "wild-born" counterparts [START_REF] Farquharson | A meta-analysis of birth-origin effects on reproduction in diverse captive environments[END_REF]. It usually involves developmental failures characterized by fertilization issues, embryonic lethalities or the occurrence of deformities. They often ensue from defects of incorporation or synthesis of the eggs' molecular content. Indeed, the abundance of these molecules, controlling embryogenesis process after fertilization, can directly be affected by modifications of extrinsic or intrinsic factors faced by females during oogenesis [START_REF] Pelegri | Maternal factors in zebrafish development[END_REF]. Among them, the maternal mRNA expression profile may thus result from genetic and epigenetic changes in the breeders' ancestors during the domestication process. It could potentially help to make the link between mechanisms described above and the reproductive success of captive populations in comparison to their wild counterparts.

There are two main ways to study the domestication process [START_REF] Boissy | Genetics of adaptation and domestication in livestock[END_REF]. First, longitudinal studies allowing a continuous follow up of a population throughout the domestication process across generations. This method is the most efficient to understand phenotypic and molecular modifications occurring at each step of the domestication process. However, it is long and difficult to perform logistically. The second way, which is commonly used, corresponds to a comparison between wild and domesticated populations. Such method has been previously used and successfully highlighted differences in several fish species, such as steelhead trout (O. mykiss) [START_REF] Christie | A single generation of domestication heritably alters the expression of hundreds of genes[END_REF],

Atlantic salmon (Salmo salar) [START_REF] Bicskei | A comparison of gene transcription profiles of domesticated and wild Atlantic salmon (Salmo salar L.) at early life stages, reared under controlled conditions[END_REF][START_REF] Bicskei | Comparing the transcriptomes of embryos from domesticated and wild Atlantic salmon (Salmo salar L.) stocks and examining factors that influence heritability of gene expression[END_REF][START_REF] Yeates | Assessing risks of invasion through gamete performance: farm Atlantic salmon sperm and eggs show equivalence in function, fertility, compatibility and competitiveness to wild Atlantic salmon[END_REF], Atlantic cod (Gadus morhua) [START_REF] Lanes | Transcriptome of Atlantic Cod (Gadus morhua L.) Early Embryos from Farmed and Wild Broodstocks[END_REF][START_REF] Lanes | Biochemical composition and performance of Atlantic cod (Gadus morhua L.) eggs and larvae obtained from farmed and wild broodstocks[END_REF], greater amberjack (Seriola dumerili) [START_REF] Zupa | Comparative Study of Reproductive Development in Wild and Captive-Reared Greater Amberjack Seriola dumerili[END_REF] and Eurasian perch (Perca fluviatilis) [START_REF] Khendek | How does a domestication process modulate oogenesis and reproduction performance in Eurasian perch?[END_REF][START_REF] Křišt'an | Comparison of reproduction characteristics and broodstock mortality in farmed and wild eurasian perch (Perca fluviatilis L.) females during spawning season under controlled conditions[END_REF]. However, numerous pieces of information are often lacking (genetic relatedness between populations, rearing conditions and history of the domestication process) and thus prevent drawing accurate conclusions. Indeed, as previously explained, phenotypic modifications potentially leaded by genetic/epigenetic modifications, depend on rearing practices. Such information is not always tracked properly by farmers. For example, for several domesticated fish stocks, wild breeders are introduced to keep a sufficiently high genetic diversity [START_REF] Teletchea | Levels of domestication in fish: implications for the sustainable future of aquaculture[END_REF] without keeping track of these introductions leading to incomplete traceability [START_REF] Khadher | Genetic characterization and relatedness of wild and farmed Eurasian perch (Perca fluviatilis): Possible implications for aquaculture practices[END_REF]. Today, with the increasing knowledge accuracy and the development of molecular tools, such information becomes important to draw proper conclusions. One way to overcome this issue consists in investigating the genetic differentiation between wild and farmed studied populations. This preliminary step would help understanding differences between populations.

In the context of fish production diversification, numerous efforts are done to domesticate a large number of new species [START_REF] Teletchea | Levels of domestication in fish: implications for the sustainable future of aquaculture[END_REF]. However, the lack of knowledge on biological and physiological needs of some species may lead to inadequate domestication attempts with deleterious consequences on the biological traits described above. The Eurasian perch (P. fluviatilis) is a promising species in aquaculture for the production diversification. It is a freshwater fish species widely distributed in Europe and in the Northern part of Asia [START_REF] Stepien | Taxonomy, Distribution, and Evolution of the Percidae[END_REF]. It has a niche market with a traditional demand in several European countries [START_REF] Overton | Commercial Production: Factors for Success and Limitations in European Percid Fish Culture[END_REF][START_REF] Steenfeldt | Current Status of Eurasian Percid Fishes Aquaculture[END_REF]. The Eurasian perch is a synchronous early spring spawner and its oogenesis induction and progression are mostly controlled by temperature and photoperiod variations [START_REF] Fontaine | Broodstock Management and Control of the Reproductive Cycle[END_REF][START_REF] Gillet | A survey of the spawning of perch (Perca fluviatilis), pike (Esox lucius), and roach (Rutilus rutilus), using artificial spawning substrates in lakes[END_REF][START_REF] Wang | Photothermal control of the reproductive cycle in temperate fishes: Photothermal control of reproduction[END_REF][START_REF] Migaud | Gamete quality and broodstock management in temperate fish[END_REF][START_REF] Abdulfatah | Effects of photothermal kinetics and amplitude of photoperiod decrease on the induction of the reproduction cycle in female Eurasian perch Perca fluviatilis[END_REF][START_REF] Sulistyo | Reproductive cycle and plasma levels of sex steroids in female Eurasian perch Perca fluviatilis[END_REF]. Consequently, manipulation of these two factors allowed defining a photothermal program largely used in Eurasian perch farms for out-of season reproduction in recirculating aquaculture systems (RAS) [START_REF] Abdulfatah | Effects of photothermal kinetics and amplitude of photoperiod decrease on the induction of the reproduction cycle in female Eurasian perch Perca fluviatilis[END_REF]. Despite this successful control of its reproductive cycle allowing out-of-season spawning, the reproduction success remains variable even if the same rearing conditions are applied to the broodstock [START_REF] Schaerlinger | Evaluation and Improvements of Egg and Larval Quality in Percid Fishes[END_REF]. It is probably due to the lack of knowledge on potential intrinsic and/or extrinsic modulating factors, including the history and details of the domestication progression experienced by populations.

In the present study, we chose to compare two Eurasian perch populations, reared in the same conditions but with different histories of domestication. We hypothesized that the level of domestication may modulate the accumulation of maternal mRNA in eggs during oogenesis, potentially impacting the embryos early development after fertilization.

Material and Methods

Origin of fish and broodstock management

Fish were handled in accordance with national and international guidelines for animal welfare protection (Directive 2010/63/EU. Agreement number: APAFIS#1390-2018031516387833 v2 accepted by the Lorraine Ethic Committee for Animal Experimentation (CELMEA) and the French Ministry of Research). Two populations of three years old Eurasian perch originating from artificial reproductions in November 2011 were used. They correspond to (i) fish at an advanced stage of the domestication process coming from breeders reared in RAS for at least seven generations (F7+ population) and (ii) Eurasian perch originating from breeders collected in the Geneva Lake at the embryonic stage and reared in RAS conditions (F1 population).

All animals came from the fish farm "Lucas Perches" (Hampont, France), which provided us the information of the presumable number of generations of the F7+ population and that their ancestors had supposedly been caught in Geneva Lake. Fish were transferred to our indoor facilities after weaning period in February 2012 (mean weight of 3.78 ± 1.07 g). Both populations were reared separately in different tanks but under the same RAS conditions (constant photoperiod (L:D -16:8), 300 lux at the water surface during the lighting period and at 20-21°C) to avoid oogenesis stimulation [START_REF] Abdulfatah | Effects of photothermal kinetics and amplitude of photoperiod decrease on the induction of the reproduction cycle in female Eurasian perch Perca fluviatilis[END_REF][START_REF] Migaud | Influence of photoperiod on the onset of gonadogenesis in Eurasian perch Perca fluviatilis[END_REF] until they reached mean weight of 287 ± 89 g. About two months before the experiment began (May 2014), 654 Eurasian perch (313 from the F7+ and 341 from the F1 populations) were transferred to the Aquaculture Experimental Platform (AEP, registration number for animal experimentation C54-547-18) belonging to the URAPA lab and located at the Faculty of Sciences and Technologies of the University of Lorraine (France). They were divided into six independent groups (three per population with an equivalent number of fish in each group). Fish were put into six identical rooms consisting of independent RAS with 3000 liters tanks. Temperature, photoperiod and light intensity were accurately controlled in each room using dedicated software. Environmental conditions were the same during acclimation phase and growing period. In order to induce gonadogenesis, breeders were submitted to a photothermal program allowing effective induction and control of the reproduction cycle [START_REF] Fontaine | Broodstock Management and Control of the Reproductive Cycle[END_REF] from August 18 th 2014 to June 22 nd 2015 (day 1 to day 309, Fig 1). Water levels of dissolved oxygen (9.53 ± 0.07 mg/L and 9.73 ± 0.12 mg/L), pH (7.55 ± 0.10 and 7.47 ± 0.00), nitrite (0.29 ± 0.14 mg/L and 0.09 ± 0.01 mg/L) and ammonium (0.67 ± 0.22 mg/L and 0.27 ± 0.02 mg/L) were monitored twice a week and kept under the respective thresholds in breeders' tanks (for F1 and F7+ populations, respectively). At the beginning of the experiment, all animals were individually tagged with P.I.T. tags (Transponder ISO 2 x 12 mm, Biolog-id) to monitor individuals all along the experiment. All fish were fed twice a day to satiation five days a week. In alternation, they were fed three days with commercial pellets (Sturgeon Grower N°5, Le Gouessant) and two days with frozen squids and shrimps (Pomona). On Saturdays and Sundays, they were fed once with commercial pellets to satiety.

Experimental design, tissue sampling and morphometric measures

In order to follow the oogenesis process, five sampling points of females were performed throughout the reproduction cycle: T0 at the beginning of the experiment allowed us to determine the initial status of breeders, T31, T87 and T154 sampling points allowed us to check the oogenesis along its progression and T253 the final status of the gonads before the spawning season (Fig 1).

Four to five females per tank were collected at each sampling point. Firstly, fish were anesthetized by immersion into a Tricaine methanesulfonate solution (120 mg/L; Sigma-Aldrich) for five minutes to collect blood from the caudal vein. Blood was centrifuged at 8000 rpm for 10 minutes in previously heparinized (4.5 mg heparin sodium salt from porcine intestinal mucosa 100KU, Sigma-Aldrich) microtubes. Plasma aliquots were conserved at -80°C for further evaluation of sexual steroids and Vitellogenin concentrations measurements.

After blood sampling, fish were killed using an overdose of Tricaine methanesulfonate (240 mg/L; Sigma-Aldrich) in accordance to European Ethical guidelines (Directive 2010/63/UE).

Total weight was measured before collecting the gonads which were weighted to calculate the gonado-somatic index (GSI=gonad weight/total weight*100) and fixed as described below for histological studies.

Evaluation of Steroids and Vitellogenin concentrations in the plasma

The 17-β-estradiol (E2, ng/mL) and testosterone (T, ng/mL) assays were performed on 50 µL of plasma of each sampled female for all sampling points (around 15 females per population/sampling point) using the DIAsource E2-ELISA kit (KAP0621, DIAsource) and the DIAsource Testosterone ELISA kit (KAPD1559, DIAsource), respectively. Samples were diluted from 1/1 to 1/50 for E2 and from 1/1 to 1/10 for T measurements depending on the oogenesis developmental stage. The E2 assay sensitivity was 5±2 pg/mL and the range of use was from 0 to 880 pg/mL with an intra assays CV varying from 4 to 17% and an inter assay CV of about 27%.

Concerning T, the sensitivity was 0.083 ng/mL and the range of use was from 0 to 16 ng/mL. The intra assays CV varied from 7 to 19 % and the inter assay CV was about 16%.

The Vitellogenin plasmatic concentration was indirectly estimated in 80 µl of plasma by measuring concentrations of the alkali-labile phosphate level as described in [START_REF] Ammar | Continuous lighting inhibits the onset of reproductive cycle in pikeperch males and females[END_REF].

Gonads histology

Female gonads from all sampling points were fixed in Bouin's solution for one week before being washed in 70% alcohol. Samples were then dehydrated in ascending series of ethanol (70-100%) before being embedded in paraffin with an orientation allowing transversal cuts. Five micrometer sections were performed with a Leitz Wetzlar microtome and collected on glass slides.

Masson's trichrome staining was done according to a protocol adapted from [START_REF] Gabe | Techniques histologiques[END_REF] as follows:

Hematoxylin solution modified according to Gill III (Merck) was used from five to ten minutes; Phloxine (VWR) was diluted in water at 0.5% and used for five minutes; Light Green (Sigma) was diluted in water at 0.5% and used from three to five minutes. Stained sections were examined, photographed and analyzed using a light upright optical microscope (Nikon Eclipse Ni-U) associated with a DS-Fi1 digital camera and the software NIS BR (Nikon France, Champigny-sur-Marne, France) at low magnification (x2 and x4).

Oocytes stages were determined according to [START_REF] Rinchard | Comparative study of reproductive biology in single-and multiple-spawner cyprinid fish. I. Morphological and histological features[END_REF] and classified into six classes: primary growth (PG), early cortical alveoli stage (ECA), late cortical alveoli stage (LCA), early vitellogenesis (EV), late vitellogenesis (LV) and atresia (A). Oogonia (O) stages were also identified.

The gonadic maturation state was determined by counting all oocytes of each class on one complete and representative transversal stained sections of the ovary for T0-T154 and three representative transversal sections for T253 because the gonads were then too large to be laid on one slide.

Gamete collection and fertilization

Before the spawning season, all females were caught to take some oocytes using a catheter and determine their oocytes maturation according to [START_REF] Żarski | A new classification of a preovulatory oocyte maturation stage suitable for the synchronization of ovulation in controlled reproduction of Eurasian perch, Perca fluviatilis[END_REF]. On May 13 rd and 15 th 2015, females from all tanks were examined and were allocated to separate tanks for the spawning season, according to their oocyte maturation stage and regardless of their original population. Thus, one tank contained the females having oocytes at stages I and II, another tank contained the females at stages III and IV (which were all coming from the F1 population), and another tank contained females that could not be staged reliably. Males from F7+ and F1 populations were kept apart in two tanks depending on their origin. When one female spawned, all females of this tank were monitored daily to identify ovulation and collect the eggs by stripping them. This procedure was always performed between 4am and 5am and each spawn was treated individually. The first spawn observed in each tank was not considered for the experiment because once in the water, the eggs are activated and their ability to be fertilized decrease rapidly. The spawning season took place from May 29 th to June 22 nd 2015 and each female stripped was identified as to its original population by its P.I.T. tag. Eggs were fertilized as described in [START_REF] Żarski | Effect of different activating solutions on the fertilization ability of Eurasian perch, Perca fluviatilis L., eggs[END_REF] with sperm from three males (total volume of 100 µl sperm/g dry eggs). Eggs stripped in the same day were fertilized using the same pool of sperm and at the end of the spawning season no day effect was seen in the reproductive performance. All males used to perform fertilization came from the F1 population, that presented a higher quality (motility estimated under a light upright optical microscope (Nikon Eclipse Ni-U)). The other part of each spawn (18 individual samples of about 200mg) was immediately frozen in liquid nitrogen and kept at -80°C for subsequent transcriptomic analysis.

Study of reproductive performance

In total, 32 spawn, corresponding to 32 females (13 from F7+ population and 19 from the F1), were treated as described above. A previous study performed in our team showed that cell cleavage timetable can be highly variable between embryos even within the same spawn [START_REF] Alix | An alternative developmental table to describe non-model fish species embryogenesis: application to the description of the Eurasian perch (Perca fluviatilis L. 1758) development[END_REF].

Thus, we choose to perform our first evaluation of the embryonic survival at 24 hours post fertilization (hpf) because it was the earliest stage to establish a relevant estimation of survival.

We were thus not able to make a distinction between non-fertilized eggs and dead embryos.

Eurasian perch eggs are surrounded by a jelly coat that protects embryos from the outside environment [START_REF] Gillet | A survey of the spawning of perch (Perca fluviatilis), pike (Esox lucius), and roach (Rutilus rutilus), using artificial spawning substrates in lakes[END_REF]. This jelly coat prevents to sort the eggs depending on their survival at each stage as it is currently performed with other fish species for which eggs are not attached. Thus, we performed the following protocol to evaluate the survival in the same ribbon samples at several timepoints.

Just after fertilization, three ribbon fragments (samples) of around 500 eggs were randomly cut from each spawn/ribbon and incubated to determine the percentage of embryos alive at different timepoints. We considered that the borders of each fragment should be avoided since embryonic development in this region could be impacted by the cutting. At each specific timepoint, around 100 embryos were counted and, among them, the number of alive embryos (those reaching the proper developmental stage) was counted to evaluate the percentage of alive embryos. These evaluations were performed in the middle of each fragment at 24, 48, 72 and 120 hpf to avoid potential border effects described above. The percentage of embryos alive was estimated using the following formula: (number of alive embryos at one stage/total number of embryos studied at that stage) x 100. In addition, three other fertilized samples of around 100 eggs/sample were kept apart, without manipulation, until the hatching period. They were used to determine the percentage of embryos hatching using the following formula: (number of hatched embryos/total number of embryos in the fragment) x 100 and deformities rates according to [START_REF] Alix | Deformities in newly hatched embryos of Eurasian perch populations originating from two different rearing systems[END_REF]. The global deformities rate (Dr) and specific deformities rates allowed studying defects in the axis (Ad), yolk (Yd), cardiac (Cd), mouth (Md), eyes (Ed) and others (Od) as described in [START_REF] Alix | Deformities in newly hatched embryos of Eurasian perch populations originating from two different rearing systems[END_REF].

RNA extraction

Total RNA was extracted from unfertilized frozen eggs from the 32 spawn (mean weight of 100 mg, 10-15 eggs) using TRIzol reagent (Life Technologies) at a ratio of 100mg per mL of reagent and following the manufacturer's instructions with some modifications. Indeed, a milling step was added during the homogenization step to get rid more easily of the gelatinous envelope and chorion surrounding eggs. To do so, a bullet blender (Next Advance) and zirconium oxide beads 1.0 mm were used. In addition, a supplementary centrifugation (4°C, 30min. 13000 rpm) was performed before the addition of the chloroform to remove the lipid content in eggs. A NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies) was used to evaluate the quantity of total RNA and an Agilent 2100 Bioanalyzer (Agilent Technologies) was used to evaluate the integrity of the RNA extracted. Samples exhibited an integrity score higher than 7 and were used for the microarray analysis.

Microarray analysis

The Eurasian perch array (SurePrint G3 Custom Gene Expression Microarray, 8x60K -Agilent Technologies) contains 48,986 non-redundant probes previously identified and available from the PhyloFish Database [START_REF] Pasquier | Gene evolution and gene expression after whole genome duplication in fish: the PhyloFish database[END_REF]. The "One-Color Microarray-Based Gene Expression Analysis (Low Input Quick Amp Labeling) Protocol" was followed for samples preparation, hybridization, washing and scanning of slides and data extractions. Briefly, 150 ng of total RNA were used for the amplification/Cy3-labeling step. After this step, samples were purified (RNeasy mini kit, Qiagen) and quantified using a NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies). Samples exhibiting a yield higher than 1.03 µg of cRNA and a specific activity higher than 10.30 pmol of Cy3/µg of cRNA were fragmented and used to hybridize arrays (one sample failed the labeling step and was excluded from the experiment). Samples (600 ng of Cy3-cRNA) were randomly distributed onto four slides. After 17h of hybridization at 65°C, slides were washed, dried and scanned with an Agilent Technologies Scanner (G2505C). Scanned images were extracted with Agilent Feature Extraction software.

Data extracted from scanned images were normalized and log(2) transformed for statistical analyses (all data are available in the Gene Expression Omnibus database under the accession code GSE119802). Using the GeneSpring software, an unpaired t-test followed by a Benjamini-Hochberg correction was used to identify genes differentially expressed (DEG) between the two populations (p<0.05). Then, a hierarchical clustering analysis (unsupervised average linkage) was performed using Cluster 3.0 software (version1.52). Clusters were visualized using TreeView software (version 1.1.6r4).

Gene ontology analyses

In order to give an overview of gene ontology (GO) terms represented among differentially expressed genes between F1 and F7+ populations, they were functionally annotated and classified using Blast2GOv4.1.9 software [START_REF] Gotz | Highthroughput functional annotation and data mining with the Blast2GO suite[END_REF]. Default parameters were used for blastx and GO annotations.

In the following, we performed a GO overrepresentation analysis in which we compared the list of DEG to a reference list corresponding to all expressed genes in the microarray. To retrieve this reference list, we first filtered all genes present in the array. Genes were, considered as expressed when they presented a signal above background in at least 75% of the samples and in at least one of the populations. In a second time, as the Eurasian perch genome is not yet available in any databases allowing performing GO analyses, we chose to retrieve human orthologs identifiers for the reference and DEG lists. It allowed us to perform the analysis. In order to find these identifiers, we aligned probes designed for the array corresponding to each sequence in the lists against the Eurasian perch transcriptome extracted from the genome, provided by [START_REF] Ozerov | Highly Continuous Genome Assembly of Eurasian Perch (Perca fluviatilis) Using Linked-Read Sequencing[END_REF], with minimap2 (version 2.7 with -m 20 parameters). Because the gene prediction file had small UTRs, we extended the prediction on both transcript sides by 2kb. If a probe had a unique alignment, it was assigned to the corresponding transcript. If the probe had two corresponding transcripts located one after the other on the genome, we assigned it to the transcript having the match closest to its center. Probes with no match or over two matches were not assigned to a transcript. We then retrieved the human orthologs identifiers from the annotation file provided by [START_REF] Ozerov | Highly Continuous Genome Assembly of Eurasian Perch (Perca fluviatilis) Using Linked-Read Sequencing[END_REF]. The genes still missing identifiers in the DEG list were manually annotated using UniProt accessions. In total, among the 358 DEG, 265 human orthologs identifiers were found. The GO overrepresentation analysis was performed using Panther14.0 [START_REF] Mi | Protocol Update for largescale genome and gene function analysis with the PANTHER classification system (v.14.0)[END_REF]. Parameters were set at default, meaning that a Fisher's exact test and a Benjamini-Hochberg correction were applied. We asked for GO-Slim Biological Processes (BP) and Pathways and only corrected p-values <0.05 were considered as significant.

Real-time PCR analysis

Genes presenting a log(2) fold change (log(2)FC) > 4 in microarray were additionally studied by real-time qPCR in all samples previously used for the microarray. After RNA extraction, a DNase treatment (DNase I, RNase free -Thermo Scientific) was applied to 5 µg of all samples (n=32) following the manufacturer's protocol. The reverse transcription was performed in a final volume of 20 µl using a M-MLV Reverse Transcriptase (Sigma-Aldrich), 1 μg of RNA and random nonamers (2.5 µM -Sigma-Aldrich) and following the manufacturer's protocol.

Reverse transcript products were diluted 1:27 and 5 µl were used for the real-time PCR, using PerfeCTa SYBR Green SuperMix (Quanta Bioscience) and 5 pmol of each primer. Primers were designed using Primer3Plus or Primer design Tool-NCBI software (S1 Table ). The real-time qPCR was performed using a Step One Plus system (Applied Biosystems, Foster City, USA). The PCR program consists of a first step at 95°C for two minutes followed by 40 cycles consisting of a denaturation step at 94°C for 15s, an annealing step at 50-58°C depending on the primers pairs for 15s and an elongation step at 72°C for 30s. The amplification was followed by a melting curve stage, according to manufacturer's parameters in order to check the primers specificity. The abundance of the target cDNA in each sample was calculated using a serial dilution of a pool of all cDNA samples using the StepOne Software (Applied Biosystems, version 2.1). This dilution curve was used to certify the reaction efficiency (80-120%). All samples were analyzed in duplicate and the geometric mean of the expression level data of the adenosine kinase-like, RNA18s, TATA-box binding and ELAV-like protein 1-like (S1 Table ) was used as reference to normalize the data obtained. These genes were found as stable in microarray and preliminary realtime qPCR analyses in all samples. The values obtained after normalization were analyzed using a Mann-Whitney test (RStudio software version 1.0.143) to compare differences of gene expression between populations (p<0.05).

Genetic variability between populations

Genomic DNA was extracted from 42 fin samples representing the two populations (21 F1 and 21 F7+), using the universal salt-extraction method, according to [START_REF] Aljanabi | Universal and rapid salt-extraction of high quality genomic DNA for PCRbased techniques[END_REF]. Purity and quantity of genomic DNA were assessed using a NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies). Eight microsatellites previously used on P. fluviatilis [START_REF] Khadher | Patterns of genetic structure of Eurasian perch (Perca fluviatilis L.) in Lake Geneva at the end of the spawning season[END_REF] were selected: PflaL1, PflaL2, PflaL4, PflaL6 [START_REF] Leclerc | Isolation and characterization of microsatellite loci in the yellow perch (Perca flavescens), and cross-species amplification within the family Percidae[END_REF], SviL7 [START_REF] Wirth | Isolation and characterization of microsatellite loci in the walleye (Stizostedion vitreum), and cross-species amplification within the family Percidae[END_REF], Svi17 [START_REF] Borer | Microsatellites in walleye Stizostedion vitreum[END_REF], YP60 and YP111 [START_REF] Li | Isolation and characterization of microsatellites in yellow perch (Perca flavescens)[END_REF]. Two multiplex amplifications were done using fluorescently labeled primers. The first multiplex (A) contained PflaL2 (FAM), PflaL4 (PET), SviL7 (VIC), Svi17 (FAM) and YP111 (PET). The second one (B) contained YP60 (FAM), PflaL1 (VIC) and PflaL6 (FAM). Polymerase chain reaction was carried out using the Multiplex TEMPase 2X MasterMix (VWR), 10 pmol of fluorescent primer mix, genomic DNA and water for a final volume of 30 µl. PCR conditions for multiplex A were: 95°C for five minutes, 28 cycles at 95°C for 30 seconds, 55°C for 90 seconds and 72°C for 30 seconds, and a final extension of 45 seconds at 60°C. For multiplex B, PCR conditions were: 95°C for five minutes, six cycles at 95°C for 30s, 48°C for 90s and 72°C for 30s, 22 cycles at 95°C for 30s, 50°C for 90s and 72°C for 30s, and a final extension at 60°C for 45s. PCR products were diluted (1:151) with deionized water and added Hi-Di™ Formamide (Applied Biosystems) and GeneScan 600 LIZ Size Standard (Applied Biosystem). The fragment analysis was performed on a 3500 Genetic Analyzer (Applied Biosystems HITACHI) and alleles were scored with Geneious 11.0.2 [START_REF] Kearse | Geneious Basic: An integrated and extendable desktop software platform for the organization and analysis of sequence data[END_REF].

Genetic diversity was estimated through calculation of observed (Ho) and expected (He) heterozygosities in GENETIX [START_REF] Belkhir | GENETIX 4.05, logiciel sous Windows TM pour la génétique des populations[END_REF]. Population differentiation was assessed by estimating the "global" FST statistic on populations through an analysis of molecular variance (AMOVA when considering only one group of populations) performed in Arlequin with 10000 permutations.

Divergence between populations was estimated with a FST-pairwise test (10000 permutations).

Statistical Analysis

Differences in GSI, sexual steroids and Vitellogenin concentrations, oocytes stages abundance on histological cut and specific deformities were estimated using a non-parametric Wilcoxon-Mann-Whitney. Normality and homogeneity of variance were tested using respectively Shapiro-Wilk and Levene's test.

In order to investigate possible differences of embryonic development between populations, percentage of embryos alive at 24, 48, 72 and 120 hpf and at hatching were compared between populations using a one-way repeated-measures ANOVA. The statistical model included as fixed effect the percentage of alive embryos for F1 and F7+ populations. It has been controlled by the introduction of the females as covariate. This model was chosen after comparison with a second one without a covariate. They were ranked according to their Akaike information criterion and the one having the lowest criterion was chosen [START_REF] Cavanaugh | The Akaike information criterion: Background, derivation, properties, application, interpretation, and refinements[END_REF]. A TukeyHSD was performed to identify differences between populations at different times individually.

For all tests, a p-value ≤ 0.05 was considered statistically significant. All values given are represented as means ± standard error of mean (SEM). All tests were performed using R (v.1.1.423) [START_REF] Core | R: A language and environment for statistical computing[END_REF]. The package stats was used to perform Shapiro-Wilk, Wilcoxon-Mann-Whitney, ANOVA and TukeyHSD tests. While car was employed for Levene's test.

Results

Two Eurasian perch populations (F1 and F7+) were used in the present study. Analysis of their genetic differentiation using microsatellites revealed a FST of 0.1055 (p < 0.001). In addition, the F7+ population presented a higher observed heterozygosity (mean=0.440) compared to the F1 (mean=0.348), indicating a larger genetic diversity of the F7+ population.

Follow up of the gonadogenesis progression reveals few differences between populations.

Five sampling points of females were performed as shown in Fig 1 . The GSI of both populations increased progressively all along gonadogenesis to reach close to 13% one month before the spawning season (S1 FigA). The only difference between populations can be seen at T31 with a higher GSI (p = 0.02) for the F7+ females (0.91 ± 0.05%) compared to the F1 ones (0.69 ± 0.06%). A follow up of the hormonal status of females did not present any significant difference of the testosterone level between populations. At T253 a higher level was recorded for the F1 (50.35 ± 5.89 ng/mL) compared to the F7+ (35.81 ± 6.49 ng/mL, p = 0.06 -S1 FigB). Similar data were obtained for the 17-β-oestradiol (E2), except for the T253 for which the hormonal level was higher (p = 0.03) for F7+ (8.71 ± 0.34 ng/mL) than for F1 (7.39 ± 0.44 ng/mL), no other statistical difference was found (S1 FigC). The follow up of the Vitellogenin level in the blood did not present any significant difference between F7+ and F1 populations (S1 FigD). These data suggest that the oogenesis progression for both populations was similar.

However, the histological study of gonads revealed that F1 females presented higher percentage of late vitellogenesis oocytes (57 ± 6%) at the end of the oogenesis in comparison to the F7+ ones (41 ± 4%; T253, Fig 2A andB), suggesting that oogenesis was slightly more advanced in the F1 than in the F7+.

The embryonic survival is higher in F1 than in F7+ spawn

Following the observation of a slightly faster oocytes development in F1 population, the first spawning was more precocious for F1 than for F7+ fish. F1 females, coming from all three original replicate tanks, spawned between thirteen and seven days earlier than F7+ first spawner (Fig 3). The main effect for population yielded an F ratio of F(1,30) = 4.266, p = 0.0476, indicating a significance difference on the number of embryos alive depending on the population.

Since an overall population effect on the percentage of embryos alive was observed, TukeyHSD tests were performed between populations and revealed significant differences between them at 48 (p = 0.02), 72 (p = 0.03) and 120 hpf (p = 0.05) and at hatching (p = 0.03).

Consequently, embryonic survivorship was significantly higher in F1 than in F7+ population from 48 hpf (Fig 4A and S2 Table). Interestingly, more heterogeneity of survivorship is seen in F7+ population at all timepoints (coefficient of variation (CV) = 62%, 88%, 95%, 90% and 121% at 24, 48, 72 and 120 hpf and at hatching, respectively) in comparison to F1 in the same stages (CV = 51%, 56%, 56%, 59% and 64% at 24, 48, 72 and 120 hpf and at hatching, respectively). The overall occurrence of deformities and that of specific deformities were comparable between populations and did not present any statistical difference (Fig 4B and S3 Table). Bars with numbers correspond to the number of spawn obtained each day for F1 and F7+ females. preceded by a significant one-way ANOVA in repeated measures (p ≤ 0.05 were considered as significant). (B) Total deformities rates (Dr) and specific rates in various tissues at hatching in F1 and F7+ populations (Ad -Axis, Cd -Cardiac, Yd -Yolk, Md -Mouth, Od -Other, Ed -Eyes).

No significant difference has been identified between populations using non-parametric Wilcoxon-Mann-Whitney test (p<0.05; n=23).

Eggs transcriptomic analysis

A large scale analysis was performed on 31 spawn to compare the maternal transcriptomic profiles of non fertilized eggs. The statistical analysis revealed 358 differentially expressed genes (DEG) between populations (S4 Table ). An unsupervised average linkage clustering analysis was performed using the expression data of these 358 DEG. It allowed splitting apart both populations revealing that 172 genes were over-expressed in the F7+ population and 186 genes were overexpressed in the F1 population (Fig 5A and S5 Table ). A functional annotation of these 172 and 186 over-expressed genes allowed mapping 86 genes for each population, among which 73 and 76 genes containing specific GO annotations were identified for F7+ and F1 populations, respectively (S6 and S7 Tables).

Using Panther software, no Pathways was enriched and the term "immune system" was overrepresented (FDR p-value = 0.00251; S8 Table ). Four of the genes represented belong to the butyrophilin or the butyrophilin-like family but two of them were up-regulated in the F1 population while the others were down-regulated in this population. For three of these genes, the log(2)FC was >2 and were thus among the 16% most differentially expressed genes between the two populations.

Indeed, among the 358 DEG only 60 presented a log(2)FC > 2 between the two populations.

Fifty DEG showed a log(2)FC between 2 and 3-fold (S4 Table ) and 10 genes had a log(2)FC higher were more abundant in the F1 population. Expression level differences were confirmed for mex3b (log(2)FC = 7.26 for the RT-qPCR), the uncharacterized protein (log(2)FC = 1.65 for the RT-qPCR), hace1 (log(2)FC = 1.65 for the RT-qPCR) and per2 (log(2)FC = 3.13 for the RT-qPCR; Fig 6). However, for bloc1s1 and nibl1 the expression levels between both populations were not significantly different by RT-qPCR even if they followed the same profile than in the microarray (Fig 6). In addition, concerning iyd not only no significantly different expression was observed, but also the profile observed by RT-qPCR was in favor to a higher expression in the F7+ population, which is in contradiction with the microarray data.

Discussion

The Differentiation Index (FST) illustrates a genetic divergence between populations ranging from 0 (gene flow between populations) to 1 (isolated populations without any gene flow).

In the literature, a lower FST (0.002, p < 0.0023) was found between various samples of Eurasian perch all around the Geneva Lake [START_REF] Khadher | Patterns of genetic structure of Eurasian perch (Perca fluviatilis L.) in Lake Geneva at the end of the spawning season[END_REF] and authors considered that there was only one population in the lake. The same authors showed a minimum FST of 0.45 (p < 0.0001) between a wild and nine farmed populations of Eurasian perch that were all supposed to originate from the same geographic location [START_REF] Khadher | Genetic characterization and relatedness of wild and farmed Eurasian perch (Perca fluviatilis): Possible implications for aquaculture practices[END_REF]. They were considered as distinct populations. In comparison, our FST (0.1055, p < 0.0001) is in between the above mentioned examples. The interpretation of significant differences when using highly variable loci, such as microsatellites, has to be very carefully analyzed since its biological meaning can often be weak [START_REF] Hedrick | Perspective: Highly Variable Loci and Their Interpretation in Evolution and Conservation[END_REF]. In the present study, we selected the most variable microsatellites according to literature on Eurasian perch microsatellite differentiation. In addition, the higher heterozygosity index for the F7+ (0.440) compared to the F1 (0.348) shows a higher heterogeneity in the F7+ population compared to the F1. Our F7+ population came from a partner fish farmer but first originated from the transfer from another fish farm. The initial stock was supposedly from Geneva lake but it was difficult to obtain a reliable traceability. Consequently, it is possible that some individuals from different geographic locations may have been introduced in the F7+ stock, which is a common practice in aquaculture [START_REF] Khadher | Genetic characterization and relatedness of wild and farmed Eurasian perch (Perca fluviatilis): Possible implications for aquaculture practices[END_REF].

Alternatively, the random independent sampling process to create each population stock (F1 and F7+) could have led to different degrees of genetic drift in each stock. In our conditions, none of these hypotheses can be ruled out.

Control of the reproductive cycle

The oogenesis progression was similar between populations even if females from the F1 population seemed to respond faster to the photothermal stimulation and kept a slight nonsignificant advance all along the process. This advance became significant one month before the spawning with a higher proportion of late vitellogenesis oocytes and may contribute to explain the early onset of the spawning season of F1 females (thirteen to seven days before the F7+ ones). For technical reasons it has not been possible to determine the end of the spawning season of F7+ fish.

The short spawning season of F7+ population in comparison to F1 is thus not relevant. It is worthy to note that at the end of the experimental period all remaining females that did not spawn presented developed gonads, once they were all slaughtered at the end of the experiment.

In this study, sexual steroids and Vitellogenin levels were either slightly or not different during the reproduction cycle between populations. The same pattern has been observed for the GSI. Differences found may either be mediated by yet unknown mechanisms or very subtle variations of hormone levels controlling the oogenesis progression. From the current knowledge, it remains difficult to precisely point which of these parameters imposed these modulations between the two populations.

Concerning the developmental success, results showed that F1 population presented an embryonic survival higher than the F7+ one with differences ranging from 16 to 26% of survival depending on the timepoint. This difference is the lowest at 24 hpf (16%), leading to a nonsignificant result at this timepoint. In any case, differences between populations may be due to egg quality issues that led either to fertilization impairments or higher mortality occurrence. These results enable us to conclude that the overall egg quality was higher in F1 than in F7+ population.

Egg quality relies greatly upon its intrinsic content and some molecules may have very early and essential effect during embryogenesis, including pronuclear congregation and mitotic spindle assembly [START_REF] Dekens | The maternal-effect gene futile cycle is essential for pronuclear congression and mitotic spindle assembly in the zebrafish zygote[END_REF]. In addition, only few exchanges occur between the developing embryos and the environment [START_REF] Schaerlinger | Evaluation and Improvements of Egg and Larval Quality in Percid Fishes[END_REF] reinforcing the importance of the egg content during early embryogenesis. Moreover, the F7+ population shows more heterogeneous survivorship. It is potentially related to individuals' history in the tank during the oogenesis (e.g. behavioral difference). In any case this difference of egg quality may be explained by the resource allocation theory [START_REF] Waaij | A resource allocation model describing consequences of artificial selection under metabolic stress[END_REF]. Even if domestication of Eurasian perch begun 20 years ago, numerous questions remain to properly achieve their breeding in recirculating aquaculture systems. Since no selective program has begun yet for this species [START_REF] Teletchea | Levels of domestication in fish: implications for the sustainable future of aquaculture[END_REF], the artificial selection driven by farmers remains empirical. In addition, rearing conditions are not fully optimized and fish farmers continue to make changes in their practices. Consequently, an imbalance between biological functions may occur in this species because resource intakes are allocated toward certain traits at the expense of other functions, such as reproduction. No zootechnical practices may be undertaken to compensate the lack of intake that may occur at each generation. This data is in accordance with numerous studies showing a decrease of the reproductive success in domesticated populations. Indeed, lower fertilization and hatching rates in farmed Atlantic cod (G. morhua) [START_REF] Lanes | Transcriptome of Atlantic Cod (Gadus morhua L.) Early Embryos from Farmed and Wild Broodstocks[END_REF][START_REF] Lanes | Biochemical composition and performance of Atlantic cod (Gadus morhua L.) eggs and larvae obtained from farmed and wild broodstocks[END_REF], lower survival rate at the eyed stage in farmed Atlantic salmon (S. salar) [START_REF] Mcginnity | Fitness reduction and potential extinction of wild populations of Atlantic salmon, Salmo salar, as a result of interactions with escaped farm salmon[END_REF], lower hatching rate in cultured common sole (Solea solea) [START_REF] Lund | A comparison of fatty acid composition and quality aspects of eggs and larvae from cultured and wild broodstock of common sole ( Solea solea L . )[END_REF] and oogenesis impairments in captive-reared greater amberjack (S. dumerili) [START_REF] Zupa | Comparative Study of Reproductive Development in Wild and Captive-Reared Greater Amberjack Seriola dumerili[END_REF] were described in comparison with their respective wild counterparts. Therefore, this problem is widely observed in fish species and further studies are needed to test the hypotheses proposed in the resource allocation theory during fish domestication. If these hypotheses are confirmed, it means that a more accurate knowledge of metabolic needs for each biological function at each step of the life cycle is required before starting any selection program. Farmers should then take this information into account in order to implement a compensation program when the resource allocation is imbalanced. In this manner, all biological functions may beneficiate from an increase of the global resource intake. In addition, the choice of selected biological traits in addition to the follow up of the fitness of individuals should be assessed properly in order to detect early any deviation.

Concerning the Eurasian perch (P. fluviatilis), other studies showed that wild populations had higher reproductive performance than captive ones but with even more drastic differences than in the present work [START_REF] Khendek | How does a domestication process modulate oogenesis and reproduction performance in Eurasian perch?[END_REF][START_REF] Křišt'an | Comparison of reproduction characteristics and broodstock mortality in farmed and wild eurasian perch (Perca fluviatilis L.) females during spawning season under controlled conditions[END_REF]. Differences between our study and these ones may be due to accumulation of the domestication effect on oogenesis and spermatogenesis, since in those studies developmental performance of embryos was evaluated from pure crosses from each population. In the present study, our goal was to investigate exclusively the effect of the domestication on females' performance. Consequently, we chose to reduce the potential effect of sperm quality and fertilized all spawn with sperm from F1 males. Another explanation could be that rearing conditions using in other studies were different (tank size and colour, first spawn or not…) from ours and influenced differently wild populations in the different studies.

Variations of gene expression profiles between wild and "domesticated" populations have been demonstrated in whole larvae and embryos of Atlantic salmon (S. salar) [START_REF] Bicskei | A comparison of gene transcription profiles of domesticated and wild Atlantic salmon (Salmo salar L.) at early life stages, reared under controlled conditions[END_REF][START_REF] Bicskei | Comparing the transcriptomes of embryos from domesticated and wild Atlantic salmon (Salmo salar L.) stocks and examining factors that influence heritability of gene expression[END_REF], fertilized eggs of Atlantic cod (G. morhua) [START_REF] Lanes | Transcriptome of Atlantic Cod (Gadus morhua L.) Early Embryos from Farmed and Wild Broodstocks[END_REF] and whole larvae of steelhead trout (O. mykiss) [START_REF] Christie | A single generation of domestication heritably alters the expression of hundreds of genes[END_REF] but no study investigated the maternal mRNA profile. These observations strengthen the necessity to better understand the contribution of maternal mRNA to the embryonic early development in our conditions.

Eggs transcriptome

Our working hypothesis was that the domestication process may impact the maternal mRNA content of the eggs and thus the transcriptomic profile of unfertilized eggs coming from F1 and F7+ populations. In total, 358 genes presented a significant difference of expression between populations.

The GO analysis revealed one biological process term overrepresented among the DEG list that we analyzed, suggesting that this function may be affected by the domestication process. This biological process corresponds to the "immune system" grouping proteins usually involved in the regulation of the adaptative immune response such as activation and proliferation of effective T cells and cytokine production [START_REF] Huang | T Cell Receptor Signaling: Beyond Complex Complexes[END_REF], controlling inflammation. However, these genes can be expressed in several tissues and may thus be involved in other biological functions. In fish, the T cell receptor signaling pathway seems to be as complex as in mammals but remains yet to be understood [START_REF] Partula | Surface markers of fish T-cells[END_REF][START_REF] Laing | Fish T cells: Recent advances through genomics[END_REF]. Some transcriptomic studies in eggs, embryos and early larvae showed that genes related with the immune system are often differentially expressed depending on egg quality or the domestication level. Concerning studies on the egg quality, genes representing the immune system are mostly down-regulated in eggs presenting lower potential to properly develop [START_REF] Żarski | Transcriptomic Profiling of Egg Quality in Sea Bass (Dicentrarchus labrax) Sheds Light on Genes Involved in Ubiquitination and Translation[END_REF][START_REF] Mommens | Profiling of the embryonic Atlantic halibut (Hippoglossus hippoglossus L.) transcriptome reveals maternal transcripts as potential markers of embryo quality[END_REF].

Concerning domestication comparisons, several genes involved in the immune system were downregulated in a domesticated population of Atlantic salmon (S. salar) compared to their wild counterparts [START_REF] Bicskei | A comparison of gene transcription profiles of domesticated and wild Atlantic salmon (Salmo salar L.) at early life stages, reared under controlled conditions[END_REF][START_REF] Bicskei | Comparing the transcriptomes of embryos from domesticated and wild Atlantic salmon (Salmo salar L.) stocks and examining factors that influence heritability of gene expression[END_REF]. In the latter case, the populations used were selected for their growth abilities and authors proposed the existence of a trade-off between growth and immune response during the selection process. In another study performed in the steelhead trout (O. mykiss) [START_REF] Christie | A single generation of domestication heritably alters the expression of hundreds of genes[END_REF], authors proposed that it could simply reflect an adaptation of individuals to the captive environment, implying an up-regulation of genes involved in this function in domesticated fish. However, no information on the direction of perturbations (up-or down-regulation of genes in domesticated or wild populations) was given. In our study, the term "immune system" involved genes that were up-regulated either in the F1 population (five genes) or in the F7+ population (four genes).

Additionally, most of the highlighted genes code for Butyrophilin proteins involved in several biological functions, on top of their role in the immune system. More particularly, they regulate the oil droplets secretion in the milk produced by mammals [START_REF] Robenek | Butyrophilin controls milk fat globule secretion[END_REF]. Eurasian perch eggs contain a large oil droplet necessary for proper embryonic development since impairments in its formation have been correlated to embryogenesis defects [START_REF] Żarski | Oocyte quality indicators in Eurasian perch, Perca fluviatilis L., during reproduction under controlled conditions[END_REF]. It would be interesting to further investigate a potential role of these proteins in the formation of oil droplets in fish eggs.

In addition, some of the genes presenting the wider variations of expression between populations are known to be involved during the early embryogenesis and may partly explain developmental defects leading to early lethalities observed in the F7+ population. The period2 gene (per2) codes for a protein belonging to the basic-helix-loop-helix-PAS (bHLH-PAS) protein family and belongs to the clock genes controlling many developmental and physiological events [START_REF] Crews | Control of cell lineage-specific development and transcription by bHLH-PAS proteins[END_REF]. However, in several mammalian species, numerous clock genes, including per2, have been observed in the developing oocytes and eggs as maternal mRNA. These mRNA are not controlled by circadian rhythms and disappear after the MZT [START_REF] Amano | Expression analysis of circadian genes in oocytes and preimplantation embryos of cattle and rabbits[END_REF][START_REF] Amano | Expression and Functional Analyses of Circadian Genes in Mouse Oocytes and Preimplantation Embryos: Cry1 Is Involved in the Meiotic Process Independently of Circadian Clock Regulation1[END_REF]. Authors proposed that per2 could be involved in the control of meiosis but, up to now, no proof has corroborated this hypothesis and further investigations are needed to study their role in this tissue. Similar data have been observed in the developing spermatogonia in mice in which clock genes are expressed without any link with the circadian rhythm. In this study, authors made the hypothesis that this gene is mainly linked to the differentiation process [START_REF] Alvarez | Non-cyclic and Developmental Stage-Specific Expression of Circadian Clock Proteins During Murine Spermatogenesis1[END_REF]. Our study shows for the first time that a per2 gene is expressed in fish eggs, suggesting similarities with mammals. In addition, it is highly expressed in eggs of the F1 population. Our study does not allow making any hypothesis about the role of per2 or if it is dependent or independent of the circadian rhythm. However, it suggests that per2 role during the gonadogenesis may be evolutionary conserved.

The gene hace1 codes for the HECT domain and Ankyrin repeat Containing E3 ubiquitinprotein ligase 1. It was first identified as a potential tumor suppressor in humans [START_REF] Anglesio | Differential expression of a novel ankyrin containing E3 ubiquitin-protein ligase, Hace1, in sporadic Wilms' tumor versus normal kidney[END_REF]. Later, absence or mutation of this gene was related with developmental issues such as some neurodevelopmental syndromes in humans [START_REF] Hollstein | HACE1 deficiency causes an autosomal recessive neurodevelopmental syndrome[END_REF], a shortening of the body axis, an inhibition of eye pigments formation and a delay in neural tube closure in xenopus [START_REF] Iimura | The E3 ubiquitin ligase Hace1 is required for early embryonic development in Xenopus laevis[END_REF]. In this last study, authors compared their data with another work performed with zebrafish and showing no clear phenotype in this species. They argued that the zebrafish study revealed the role of HACE1 by using splicesites morpholino efficient to repress zygotic but not maternal RNA [START_REF] Daugaard | Hace1 controls ROS generation of vertebrate Rac1-dependent NADPH oxidase complexes[END_REF] while their study involved translational blocking morpholinos blocking both maternal and zygotic RNA. Thus they suggest that phenotypes in xenopus could be linked to the maternal pool of RNA. Recently, a study demonstrated that HACE1 is involved in the normal development and proper function of the heart in zebrafish [START_REF] Razaghi | hace1 Influences zebrafish cardiac development via ROS-dependent mechanisms[END_REF]. However, authors used splice-site morpholinos suggesting that this cardiac phenotype could be due to later zygotic expression of HACE1. Consequently, even if no maternal role of HACE1 has been studied in the zebrafish, studies performed on xenopus suggest that a maternal expression of these mRNA have consequences on the neural tube development. Up to now, two main targets of Hace1 have been identified: Rac1 involved in the cell cycle control and NADPH oxidase regulating the reactive oxygen species production. A previous study showed that some enzymes involved in the control of the oxidative stress present a variation of expression depending on oocytes quality in Eurasian perch [START_REF] Castets | Combined analysis of Perca fluviatilis reproductive performance and oocyte proteomic profile[END_REF]. In some marine invertebrate species, a modification of redox homeostasis may help to avoid the polyspermy during fertilization [START_REF] Schomer | Redox Changes during Fertilization and Maturation of Marine Invertebrate Eggs[END_REF]. In our study, the mRNA abundance of hace1 is higher in the F7+ population, potentially accelerating NADPH oxydase degradation and thus influencing fertilization. With the current knowledge, no transduction pathway could be favored.

Similarly, the mRNA abundance of mex3b mRNA was higher in the domesticated population. Mex3b, or muscle excess 3, codes for an evolutionary conserved RNA-binding protein involved in post-transcriptional regulations [START_REF] Buchet-Poyau | Identification and characterization of human Mex-3 proteins, a novel family of evolutionarily conserved RNA-binding proteins differentially localized to processing bodies[END_REF]. It is associated to proper embryonic development by establishing antero-posterior patterning in Tribolium [START_REF] Kimelman | Anterior-posterior patterning in early development: three strategies[END_REF], Caenorhabditis elegans [START_REF] Draper | MEX-3 Is a KH Domain Protein That Regulates Blastomere Identity in Early C. elegans Embryos[END_REF] and Xenopus laevis [START_REF] Takada | The RNA-binding protein Mex3b has a fine-tuning system for mRNA regulation in early Xenopus development[END_REF]. As a whole, embryonic patterning is regulated by expression and spatial distribution of many transcripts. In xenopus, mex3b mRNA presents a long conserved untranslated 3'UTR involved in its auto regulation [START_REF] Takada | The RNA-binding protein Mex3b has a fine-tuning system for mRNA regulation in early Xenopus development[END_REF]. In turn, the protein is involved in the mRNA stability and regulates the abundance of several genes, involved in diverse cellular functions [START_REF] Takada | The RNA-binding protein Mex3b has a fine-tuning system for mRNA regulation in early Xenopus development[END_REF]. In the present work, the maternal expression of mex3b mRNA presents the highest variation of expression level between the two populations (log2FC > 7 in the RT-qPCR). Consequently, the high expression of this transcript in the F7+ eggs may be linked to fine tuning of numerous molecular functions in the embryo and lead to diverse phenotypes.

Finally, a mRNA coding for an uncharacterized protein was found to be significantly more expressed in the eggs laid by F7+ females. BLAST analyses against the Uniprot and NCBI databases showed that they are highly homologous to other uncharacterized proteins in other fish species presenting homologies with some domains of ADP ribosylation enzymes. However, these sequences have not been identified yet.

Conclusions

Our study showed that reproductive performance of Eurasian perch females may be influenced by the domestication process which is probably closely related to the rearing practices potentially leading to several genetic and/or epigenetic modifications. This study revealed that even if the breeders of two Eurasian perch populations (F1 vs. F7+) were reared and induced in the same conditions, the F1 population started to spawn earlier than the F7+ during the spawning season. In other words, it shows that, in our conditions, the domestication process may influence the responsiveness of females to the reproductive environmental stimuli in captive environment.

The F1 population produced eggs having a better potential to develop properly until hatching compared to eggs from the F7+ population. These differences in egg quality may be linked with the important variation in the eggs transcriptomic content between populations. The identification of several genes presenting distinct expression between the two populations could open new paths of investigation to understand their role and mechanism of regulation during embryogenesis and depending on the domestication level.

Finally, the genetic differentiation analysis between studied populations did not allow us to isolate the domestication as the only factor explaining our data. It reinforces the necessity of studying populations presenting a clearer life history to further understand the dynamic of modifications occurring during the domestication process. It is particularly important in the perspective of a growing pressure toward fish farmers and scientists to initiate selection programs for several fish species.

Fig 1 .

 1 Fig 1. Graphic summary of the experimental design. The photothermal program was used to

Fig 2 .

 2 Fig 2. Histological follow up of the oogenesis progression. (A) Percentage of each cellular stage

Fig 3 .

 3 Fig 3. Timetable of the spawn obtained for both populations during the spawning season in

Fig 4 .

 4 Fig 4. Evaluation of the embryonic developmental success in F1 and F7+ populations. (A)

Fig 5 .

 5 Fig 5. Heatmaps showing differentially expressed genes between F1 and F7+ populations. (A)

than 4 -

 4 fold (5 over-expressed in the F7+ population and 5 in the F1; Fig 5B). We choose to check the expression level of the genes having a log(2)FC > 4 by real-time qPCR. Among them, three sequences could not be identified because the probes actually recognized contigs grouping numerous unidentified genes in the Eurasian perch transcriptomic database (Fig 5B). Among the remaining genes, four were more expressed in the F7+ population: mex3b (log(2)FC = 8.26), bloc1s1 (log(2)FC = 4.91), an uncharacterized protein (log(2)FC = 4.57) and hace1 (log(2)FC = 5.53). In the same way, per2 (log(2)FC = 5.68), nibl1 (log(2)FC = 4.04) and iyd (log(2)FC = 4.62)

Fig 6 .

 6 Fig 6. Compared analysis of seven genes among the most differentially expressed

Acknowledgements

We thank, Paul Georges and Franck Bertrand for their involvement in some experiments as trainees, Alain Iuretig for the care given to the breeders during the growing period and Yves Le Roux for the assistance with the survival statistical analysis. We are also grateful to Mikhail Ozerov and Freed Ahmad from University of Turku for kindly providing the annotated Eurasian perch genome used in GO analysis and the members of the UR AFPA DAC and LPGP Rennes for their support.

No significant differences have been identified between the two populations using non-parametric Wilcoxon-Mann-Whitney test (p<0.05; n=23). Population means ± SEM are presented.

S4 Table. Differentially expressed genes list. PhyloFish gene ID, gene description, False

Discovery Rate (p-value) and log(2)FC of the 358 differentially expressed genes on eggs from the F7+ and F1 populations.

S5 Table. Differentially expressed genes listed in the same order as they appear in


Mohamed Sfar

Slah Yaacoubi

Sofiene Dellagi

Mahjoub El Mountassir

Sadok Turki

-Agadir -Morocco "New advances and challenges for sustainable and smart industries" TOWARDS A PREDICTIVE MAINTENANCE BASED ON ULTRASONIC GUIDED WAVES NON-DESTRUCTIVE TESTING TECHNIQUE: A CASE STUDY

Keywords: Predictive maintenance, Ultrasonic Guided Waves, Prognostic, Residual Useful Life, ARIMA, EMD, Neural Network, Non-destructive testing

This paper concerns the predictive maintenance (PM), which is in increasing use in many fields. It proposes the use of a non-destructive testing technique, namely ultrasonic guided waves, in order to reinforce the PM of a given installation. To achieve PM, various sensors are required to collect data, which are needed as inputs of predictive models. These sensors are employed to monitor the health (said also in other words degradation) of the installation to be maintained in operation. In this study two models based on EDM-ARIMA and neural network are built and applied on data collected on a tube via the ultrasonic guided waves. Residual useful life of the structure is provided in different cases.

INTRODUCTION

The rise of digitalization and industry 4.0 promises an optimization of industrial production based on data. Numerous players in industrial analytics guarantee to effectively tackle issues just by gathering and analyzing data. Lately, more and more data is acquired and more and more analytical tools become available [START_REF] Domingo Llorente Rivera | Towards a Predictive Maintenance Maintenance System of a Hydraulic Pump[END_REF].

Nowadays, mastering the operation of industrial systems represents a big challenge for maintenance operators, since these systems are more and more complex with technological development. Therefore, many companies are now implementing predictive maintenance strategies which is a proactive approach that aims to anticipate faults and breakdowns of an industrial system, based on monitoring modules, diagnosis, prognosis and decision support in the form of an intervention calendar.

In terms of predictive maintenance, the last stage of the Structural Health Monitoring (SHM) process, which is the prognostics, is used to anticipate the evolution of the defect or degradation in the future. A SHM system is defined as an integrated part of the structure to be monitored and should ideally be able to detect, locate and evaluate damage in the structure, estimate its severity and monitor its evolution with time. The utilization of SHM will result in (i) diminished basic load by changing design principles; (ii) decreased maintenance expenses; and (iii) expanded accessibility, however it is recognized at the same time that SHM technologies will be presented only if they have direct financial advantages. Major design aspects that govern any SHM system are (i) changes of the monitored physical properties due to damage; (ii) transducer-structure interaction to reliably measure the expected changes with the required resolution; and (iii) signal processing and analysis tools to extract the required information on structural integrity [START_REF] Veidt | Non-destructive evaluation (NDE) of aerospace composites: structural health monitoring of aerospace structures using guided wave ultrasonics[END_REF].

Prognostics essentially attempt to anticipate how much time remains before a fault or failure will happen, giving the current state of the asset. In other words, prognostics are based on predicting the residual useful life (RUL) of a system based on condition monitoring of the asset and subsequent analysis of the information and data collected for a given task occurrence. In many situations, the complexity of the systems observed does not allow to derive robust and precise prognosis models. However, historical data that capture the life signal behavior of measured signals or extracted features from the incipient fault stage to equipment failure are often available. In such cases, datadriven methods that model how such signals and features evolve can be utilized to generate predictions of RUL.

Data collection of physical factors (temperature, vibration, etc.) can be achieved using vibration analysis, thermography etc. [START_REF] Galar | Essential Electronic Tools Efficiency[END_REF].

In the present work, data were collected using Ultrasonic Guided Wave (UGW). UGW are mechanical waves, at frequencies over the scope of human perceptible frequencies that propagate along a stretched structure while guided by its limits [START_REF] Yaacoubi | Ultrasonic guided waves for reinforced plastics safety[END_REF], its use to the field of preventive maintenance is in persistent increment. Today, UGW is assuming a significant economic position, particularly in the fields where the wellbeing of people is of higher significance than financial cost (i.e.: aeronautics, trains ...). They have been generally used to assess the trustworthiness of engineering structures, because they are capable to spread over significant separations and fantastic affectability to the nearness of deformities in the proliferation way [START_REF] Zhang | An integrated numerical model for investigating guided waves in impact-damaged composite laminates[END_REF]. Furthermore, they contain rich information MOSIM'20 -November 12-14, 2020 -Agadir -Morocco regarding the mechanical properties along their path of propagation and possess high sensitivity to changes in material and geometric in homogeneities due to damage. UGW is applied as an active as well as passive monitoring techniques [START_REF] Croxford | Strategies for Guided-wave Structural Health Monitoring[END_REF], where the latter is comparable to an autonomous acoustic emission measurement system. But the majority of GW SHM systems in development use active transducers to excite specific stress wave modes to interrogate the component. The prognostics rely on the approach of projection methods, which projects the current level of degradation in the future. It's essentially about a time series prediction problem. Many approaches are proposed for time series forecasting like: Artificial Neural Network (ANN), exponential smoothing model (ETS), Auto Regressive Integrated Moving Average (ARIMA), machine learning, Support vector Regression (SVR), etc. [START_REF] Hongzhan | Hybrid of ARIMA and SVMs for Short-Term Load Forecasting[END_REF][START_REF] Suhermi | Roll motion prediction using a hybrid deep learning and ARIMA model[END_REF][START_REF] Bhardwaj | A comparative study of wave-let-based ANN and classical techniques for geophysical time-series forecasting[END_REF][START_REF] Sibaramapanigrahi | A hybrid ETS-ANN model for time series forecasting[END_REF].

In this study, we focus on two approaches. The first one is a hybrid approach called Empirical Mode Decomposition-ARIMA (EMD-ARIMA), since EMD is suitable for nonlinear non-stationary data [START_REF] Huang | The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-Stationary Time Series Analysis[END_REF][START_REF] Feng | Country risk forecasting based on EMD and ELM: evidence from BRICS countries[END_REF]. EMD is used to decompose the original data into components in order to reduce the difficulty to realize the high-precision predictions, and ARIMA predicts components respectively [START_REF] Liu | An EMD-recursive ARIMA method to predict wind speed for railway strong wind warning system[END_REF]. The second one is a statistical approach, mainly used in machine learning, called Auto-Regressive Neural Networks (ARNNET). Its concept is inspired from the biological neural network of the central nervous system. Neural Networks are better used to estimate functions based on a huge volume of training data. These are capable of modelling complex and non-linear functions [START_REF] Sena | TO FORECAST PER CAPITA DISPOSABLE INCOME[END_REF].

The remainder of the paper is organized as follows. The next section provides some basic theoretical aspects of EMD-ARIMA and ARNNET. Section 3 describes the datasets used in the study and discusses the methodology used in the time series modeling and forecasting. The results obtained in this study are presented in section 4.

MATHEMATICAL MODEL

EMD-ARIMA Models

In this section we will discuss a hybrid EMD-ARIMA method of prediction. The computational framework of this method is given in Figure 1. As shown in Figure 1, the method consists in: (a) using the EMD to decompose the original time data into components with different features (sub-layers); (b) applying ARIMA models to each sub-layer and adopt the built models to forecast and (c) summarize the predictions of the sub-layers to get the final prediction for the original data [START_REF] Liu | An EMD-recursive ARIMA method to predict wind speed for railway strong wind warning system[END_REF].

EMD

EMD is an adaptive time series decomposition technique based on Hilbert-Huang transformation (HHT), which is suitable for nonlinear non-stationary data [START_REF] Huang | The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-Stationary Time Series Analysis[END_REF]. EMD can decompose the original data into a series of components with different frequencies, namely, intrinsic mode function (IMF) [START_REF] Feng | Country risk forecasting based on EMD and ELM: evidence from BRICS countries[END_REF] using the equation as follows:

𝑋(𝑡) = ∑ 𝐶 𝑖 (𝑡) + 𝑅(𝑡) 𝑛 𝑖=1 (1)
Where{𝐶 𝑖 (𝑡)}, (i=1...n) is the IMFs in different decompositions, {R (t)} is the residue and n is the number of IMF's. The computational steps of the EMD are given as follows:

Step 1: Identify all the local extrema (maxima and minima) of series {X (t)}.

Step 2: Connect all the local maxima by a cubic spline to generate its upper envelop {Xup(t)}, and all the local minima to obtain the lower envelop {Xlow(t)}.

Step 3: Calculate the mean envelop {M(t)} as follows:

𝑀(𝑡) = [𝑋 𝑢𝑝 (𝑡) + 𝑋 𝑙𝑜𝑤 (𝑡)] 2 (2)
Step 4: Extract the details as follows:

𝑍(𝑡) = 𝑋(𝑡) -𝑀(𝑡) (3)
Step 5:

Check whether {Z(t)} is an IMF: (a) if it is then set C(t)=Z(t) and set X(t)= R(t)= X(t) -C(t) ; (b) if not,
replace {X(t)} with {Z(t)} then repeat the steps 2-4 until the termination criterion is satisfied. The equation of the termination condition can be represented as follows:

∑ [𝑍 𝑗-1 (𝑡) -𝑍 𝑗 (𝑡)]² [𝑍 𝑗-1 (𝑡)]² 𝑚 𝑡=1 ≤ 𝜹 (4) (j=1, 2..., t=1,2,…,m)
Where "m" is the number of the data points, "𝜹" is the terminated parameter, and "j" is the time of iterative calculation.

Step 6: Repeat steps 1-5 until all IMFs are found [START_REF] Liu | An EMD-recursive ARIMA method to predict wind speed for railway strong wind warning system[END_REF]. MOSIM'20 -November 12-14, 2020 -Agadir -Morocco

ARIMA

Once all IMFs are found, we apply ARIMA model to each IMF to fit the time series. An ARIMA model uses autoregressive moving average (ARMA) model to fit stationary time series. When the time series are non-stationary, they must be converted into stationary series by difference transformation before using ARMA for modeling. This is what "I" in ARIMA refers to; integrated. The difference transformation can be divided into the order difference transformation and the periodic difference transformation. The first-order difference transformation is:

∇𝑌 𝑡 = 𝑌 𝑡 -𝑌 𝑡-1 = (1 -𝐵)𝑌 𝑡 (5)
where, ∇ is the difference operator, 𝑌 𝑖 value, B is the backshift operator.

Occasionally, the differenced data won't appear to be stationary and it is necessary to difference the data a second time to obtain a stationary series, and the second-order difference transformation is:

∇ 2 𝑌 𝑡 = ∇Y 𝑡 -∇𝑌 𝑡-1 = 𝑌 𝑡 -2𝑌 𝑡-1 + 𝑌 𝑡-2 (6)
The periodic difference transformation is as follow:

∇ 𝑠 𝑌 𝑡 = 𝑌 𝑡 -𝑌 𝑡-𝑠 (7)
where, ∇ 𝑠 is the periodic difference operator and s is the number of periods.

When the differentiated series is stationary, we can use ARIMA model and it can be formulated as follows [START_REF] Ohyver | ARIMA Model for forecasting the Price of Medium Quality Rice to Anticipate Price Fluctuations[END_REF]:

𝜑(𝐵)(1 -𝐵) 𝑑 𝑌 𝑡 = 𝜃(𝐵)𝑒 𝑡 ( 8
)
Where,𝜑(𝐵) = 1 -𝜑 1 𝐵 -⋯ -𝜑 𝑝 𝐵 𝑝 , 𝜃(𝐵) = 1 -𝜃 1 𝐵 -⋯ -𝜃 𝑞 𝐵 𝑞 , 𝜑 𝑖 (i=1..p) and 𝜃 𝑖 (i=1..q) are the model parameters that are estimated using the autocorrelation function (ACF) and the partial autocorrelation function (PACF); at is the error; p and q are the autoregressive and the moving average order, respectively and d is the order of differentiation [START_REF] Hongzhan | Hybrid of ARIMA and SVMs for Short-Term Load Forecasting[END_REF].

Auto-Regressive Neural Networks

In this section we will introduce the Auto-Regressive Neural Network model. It is a network of multiple nodes working in parallel as shown in Figure 2. The inputs are combined through a linear function defined by :

𝑧 𝑗 = 𝑔 (∑ 𝑦 𝑖 𝑤 𝑖𝑗 -𝛽 𝑗 𝑖 ) ( 9
)
Where 𝑧 𝑗 is the output of the j th node, 𝑦 𝑖 is the i th input, 𝛽 𝑗 is the bias of the j th node, 𝑤 𝑖𝑗 is the i th weight in the j th node.

The result of the combination is then passed through the non-linear sigmoid activation function g() which transforms the input to the node to an output which in turn acts as an input for the nodes in the next layer. The activation function is given by [START_REF] Akhter | A prediction-based approach for stock returns using autoregressive neural networks[END_REF][START_REF] Oskar | A SIMPLE AUTO-REGRESSIVE NEURAL NETWORK FOR TIME-SERIES[END_REF] :

𝑔(𝑥) = 1 1 + 𝑒 -𝑥 (

EMD-ARIMA

Mainly, ARIMA's forecasting task is selecting the appropriate model order, which is the value of p, d and q. In this study, the following steps are used to identify the model used for the prediction:

Step 1: As explained in section 2, we adopt the Empirical Mode Decomposition (EMD) to decompose the degradation indicator samples into Intrinsic Mode Functions (IMFs). Figure 4 is an example of data decomposition into IMFs.

Step 2: In the step of model identification, all of the decomposed degradation indicator has their stationarity checked to decide the "d" value in ARIMA (p, d, q). If an IMF has been found non-stationary, a difference will be executed for this IMF. In this study, in order to check the stationarity of time series we compare the p-values of the series by applying the Augmented Dickey Fuller test (ADF Test). If the p-value, using the function is lower than 0.05, means the time series is stationary and vice versa [START_REF] Fanoodia | Reducing demand uncertainty in the platelet supply chain through artificial neural networks and ARIMA models[END_REF][START_REF] Kim | Unit Roots in Economic and Financial Time Series: A Re-Evaluation at the Decision-Based Significance Levels[END_REF].

Step 3: We compute and examine the Autocorrelation function (ACF) and partial autocorrelation function (PACF) of each IMF to decide whether to include an Auto-Regressive (AR) component or Mobile Average (MA) one or both [START_REF] He | Epidemiology and ARIMA model of positive-rate of influenza viruses among children in Wuhan, China: A nine-year retrospective study[END_REF].

Step 4: In order to obtain the model and predict the data for each IMF, we split the data into two parts, which are training data and validation data. For the torsional mode, we consider the first 910 data as training data. The remaining data is used as validation data. On the other side, for the flexural mode, we consider the first 950 data as training data and the other 109 data as validation data. For prediction accuracy, the Root Mean Squared Error (RMSE) is introduced to appraise and compare the different simulation results. The RMSE is calculated using the following equation.

𝑅𝑀𝑆𝐸 = √ ∑ (𝑥 𝑛+𝑖 -𝑥 ̂𝑛(𝑖))² 𝑁 𝑖=1 𝑁 ( 11
)
Where N denotes the length of the validation data, 𝑥 𝑛+𝑖 denotes the i-th actual value of validation data, 𝑥 ̂𝑛(𝑖) denotes the i-th forecast. It is worth noting that, p, q and d may be also selected via an information criterion such as the Akaike information criterion (AIC) calculated as follow [START_REF] Ramos | Performance of state space and ARIMA models for consumer retail sales forecasting[END_REF]:

𝐴𝐼𝐶 = -2 log(𝐿) + 2𝑘 ( 12
)
Where L is the likelihood of the model and k is the total number of parameters and initial states that have been estimated. This criterion is used to measure the goodness of fit and determine the choice of optimal lags in the model. A negative value of AIC means that the model is well built [START_REF] Young | Recursive Estimation and Time Series Analysis: an Introduction[END_REF].

Step 5: Once the model is chosen for each IMF, we sum all the forecasted data using ARIMA and compare it to the actual values of the validation data.

Step 6: We set a threshold of the degradation indicator and determine at what time the forecasted values reach it.

Figure3. EMD results of the original data of the torsional mode. Each IMF is normalized by its maxima. MOSIM'20 -November 12-14, 2020 -Agadir -Morocco

ARNNET

The effective execution of a Neural Network is reliant on proper training of the network with large training sets. In the process of training, the algorithm learns to adjust the weights and bias associated with each connection. In this paper we consider the notation NNAR (p, k) to indicate that in the model there are p lagged inputs and k nodes in the hidden layer [START_REF] Sena | TO FORECAST PER CAPITA DISPOSABLE INCOME[END_REF][24]. The steps followed to build the neural network model are the following:

Step 1: Determine the order of auto-regression (p) for the dataset. In Figure 5, PACF plot is drawn for the dataset.

The abscissa represents the lags and the ordinate represents the partial autocorrelation values for the corresponding lags. We search for the last signifiant spike. In our case, the value of p is 19 and 25 consecutively for Torsional and Flexural data.

Step 2: Data has been partitioned into training and validation data. We used the same partition as the one used for EMD-ARIMA method.

Step 3: Build the ARNNET of the data using the function "nnetar" in R.

Step 4: We set a threshold of the degradation indicator and determine at what time the forecasted values reach it.

RESULTS AND DISCUSSION

The results of the EMD of the original degradation indicator for the torsional mode is given in Figure 4. From Figure 4, it can be seen that the time series have been converted successfully seven IMFs and a residue for the torsional mode. Concerning the flexural mode, the time series was decomposed into eight IMFs and a residue.

The accuracy of some of the estimated results of the tested ARIMA models applied for each IMF prediction as well as their AIC is given in Table 1 and 2. The model with minimum AIC is often the best model for forecasting.

Table1. Accuracies and AIC of some tested models for the flexural mode.

Table 2. Accuracies and AIC of some tested models for the torsional mode. MOSIM '20 -November 12-14, 2020 -Agadir -Morocco The accuracy of the estimated results of the chosen ARIMA model applied for each IMF for both torsional and flexural mode is highlighted also in Table 1 and2. For the torsional mode, if the original data is not decomposed and we apply ARIMA models, the suitable model turns out to be ARIMA (4, 1, 2) with AIC= 545.62 and RMS= 3.17. For our data, the RMSE value using EMD-ARIMA is equal to 1.26, which is 60% lower than the RMSE given by applying ARIMA directly. Same as the flexural mode, if we apply ARIMA directly without decomposing the data, the suitable ARIMA model is ARIMA (4, 2, 2) with AIC=-2255.66 and RMS= 1.21. For our data, the RMSE value using EMD-ARIMA is equal to 0.514, which is 57.5% lower than the RMSE given by applying ARIMA directly. This result shows that EMD-ARIMA is better than ARIMA which is not surprising since it decreases the forecast difficulty by treating each IMF at once. Producing estimates of confidence interval is also an important aspect. We also evaluated the performance of the different ARIMA models in producing forecast intervals. Table 3 and 4 shows the confidence intervals at 95% of some of the tested ARIMA models comparing to the actual values. From table 3, we can note that the suitable ARIMA model for the torsional mode comes generally with the tighter confidence interval. However, from Table 4, we can note that the suitable ARIMA model for the flexural mode comes generally with the larger confidence interval.

As mentioned in section 3.2, the last step is to fix a threshold of the degradation indicator and see when it is reached. In our case, the threshold is fixed to the highest value reached by the original data, which is 1.1387e-13 for the torsional mode, and 3.037e-14 for the flexural mode. According to the forecasted results using the chosen model, these values are reached at day 157 for the torsional mode and the day 147 for the flexural mode, from the start of the test. Based on the original data, the actual value should be 141 days for both modes. We note that there is a 16and 6-days gap consecutively for the torsional and flexural mode. In Figure 6 and 7, the original curve is represented in black, the curve in light blue represents the forecasted data with the chosen model. The dark blue and the green curve are simultaneously the prediction of one of the tested models and the prediction using directly ARIMA on the original data.

ARIMA (p, d, q)

For the same threshold, we adopt the ARNNET to predict the RUL. According to the forecasted results, this value is reached at day 137 for the Flexural mode which is 4 days before the actual value. However, the forecasted values of the Torsional data didn't reach the threshold.

The RMSE of the forecasted data is 0.11 and 0.148 consecutively for Torsional and Flexural mode.

In Figure 8 and 9, the original curve is represented in black and the curve in green represents the forecasted data using ARNNET method.

The tables 5 and 6 summarize the values of the residual life as well as the RUL found in this research for both torsional and flexural mode. We can conclude that EMD-ARIMA provides better results than ARIMA for both modes.

CONCLUSION

Accurate residual useful life forecasting can have a great impact on risk management and predictive maintenance. The purpose of this work is to adopt and compare two methods of time series prediction. The EMD-ARIMA and the ARNNET model were used to predict the residual useful life of tubular products and test its accuracy. The data provided in this study are the evolution of the degradation indicator of a torsional and flexural mode over time. How to effectively model these series and how to improve the quality of forecasts are the outstanding questions.

As one of the most popular forecasting methods, ARIMA turns out to be less effective than the hybrid method EMD-ARIMA. EMD-ARIMA decomposes the time series into IMFs and treats each IMF independently in order to reduce the forecasting complexity.

The hybrid method forecast is able to follow the pattern of the real data better than the ARIMA model but with a certain gap. It shows that the hybrid methodology could significantly increase the forecast accuracy. MOSIM '20 -November 12-14, 2020 -Agadir -Morocco The ARNNET is computationally more intensive than the EMD-ARIMA as it has to estimate the weight and bias corresponding to each node, to come up with better results only for the Flexural mode. However, the ARNNET did fit a good model to the torsional data.

Figure 1 .

 1 Figure 1. Framework of EMD-ARIMA models.

  Our study aims to predict the RUL of tubular structures in order to schedule predictive maintenance, so we are brought to predict the evolution of the degradation indicator (DI) over time. The datasets used in this study are two time series of the degradation indicator for two mode of propagation of UGW in pipes which are torsional and flexural. The dataset contains 1059 observations for each the torsional and flexural modes. This degradation indicator was established over the collected UGW signals basing on singular value decomposition of the matrix of reference signals acquired from the pipe. The time series' plot of the datasets is shown in Figure 3.

Figure 3 .Figure 2 .

 32 Figure 3. Degradation indicator for torsional and flexural mode. Each signal is normalized by its maxima.

Figure 4 .

 4 Figure 4. IMF's of the Torsional mode.

Figure 5 .

 5 Figure 5. PACF of the Torsional and Flexural data.

Figure 7 .

 7 Figure 7. Examples of forecasted data for flexural mode.Normalized by the maxima of the original data.

Figure 9 .

 9 Figure 9. Forecasted data of the Torsional mode using ARNNET.

Figure 8 .

 8 Figure 8. Forecasted data of the Flexural mode using ARNNET.

Table 5 .

 5 Residual life and RMSE for the torsional mode.

	Flexural mode	Actual Value	Chosen EMD-ARIMA	Tested EMD-ARIMA	ARIMA (4,2,2)	NNETAR
	RMSE		0.126	0.16	0.317	0.11
	RUL	141	147	152	ϕ	ϕ
	Torsional mode	Actual value	Chosen EMD-ARIMA	Tested EMD-ARIMA	ARIMA (4,1,2)	NNETAR
	RMSE		0.514	0.73	1.21	0.148
	RUL	141	157	157	ϕ	137

Table 6 .

 6 Residual life and RMSE for the flexural mode.

Confidence interval

Lower Upper IMF1

(5,0,1) -0.04154 0.11242 (5,0,4) -0.04677 0.107095 (5,