Audrey Barranger
email: audrey.barranger@ifremer.fr

Christophe Klopp

Barbara Le Bot

Gaëlle Saramito

Lise Dupont

Stéphanie Llopis

Claudia Wiegand

Françoise Binet

Insights into the molecular mechanisms of pesticide tolerance in the Aporrectodea caliginosa earthworm

Keywords: Aporrectodea, transcriptomics, pesticides residues, mitochondrial respiratory chain, agrosystem, soil biodiversity

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

Pesticides are specifically designed to regulate the population dynamics of unwanted organisms such as weeds, pathogens or pests. Post-war they became a major foundation of worldwide agricultural performance [START_REF] Tang | Risk of pesticide pollution at the global scale[END_REF]. As they are applied directly to soil as sprays, drenches, and granules or, increasingly, as seed coatings [START_REF] Gunstone | Pesticides and Soil Invertebrates: A Hazard Assessment[END_REF], most topsoils in agrosystems (up to 83% in Europe (V. [START_REF] Silva | Pesticide residues in European agricultural soils -A hidden reality unfolded[END_REF]) contain one or more pesticide residues, causing adverse effects to ecosystem and human health [START_REF] Alpizar | Global Chemicals Outlook II-From Legacies to Innovative Solutions: Implementing the 2030 Agenda for Sustainable Development[END_REF][START_REF] Bernhardt | Synthetic chemicals as agents of global change[END_REF]. Soils are arguably the most complex and biodiverse ecosystems on Earth, containing nearly a quarter of the planet's biodiversity, which play crucial roles in many ecosystem functions and services [START_REF] Bardgett | Belowground biodiversity and ecosystem functioning[END_REF].

In agricultural landscapes, soil sustainability and fertility depend on its biological component [START_REF] Swift | Biodiversity and ecosystem services in agricultural landscapes-are we asking the right questions? Agriculture[END_REF]. Earthworms, the first animal biomass on Earth, support numerous ecological functions: by foraging, they fragment and decompose organic litter [START_REF] Phillips | Global data on earthworm abundance, biomass, diversity and corresponding environmental properties[END_REF]; by creating biopores, they affect soil structure, water permeability and aeration [START_REF] Bastardie | Assessment of earthworm contribution to soil hydrology: a laboratory method to measure water diffusion through burrow walls[END_REF][START_REF] Pérès | Relationships between earthworm communities and burrow numbers under different land use systems[END_REF], by ingesting soil, they stimulate microbial communities [START_REF] Kersanté | Interactions of earthworms with Atrazine-degrading bacteria in an agricultural soil[END_REF] and regulate pesticide biodegradation [START_REF] Monard | Combined effect of bioaugmentation and bioturbation on atrazine degradation in soil[END_REF][START_REF] Monard | Relationship between bacterial diversity and function under biotic control: The soil pesticide degraders as a case study[END_REF]. The survival of earthworm populations and the maintenance of a high biomass are crucial in view of these ecosystem services, which are part of their engineering function. Agricultural intensification (including habitat loss, soil tillage and agrichemical use) has, however, been identified as a major driving factor behind the observed decline of terrestrial invertebrates [START_REF] Forister | Declines in insect abundance and diversity: We know enough to act now[END_REF][START_REF] Hallmann | More than 75 percent decline over 27 years in total flying insect biomass in protected areas[END_REF][START_REF] Sánchez-Bayo | Worldwide decline of the entomofauna: A review of its drivers[END_REF][START_REF] Seibold | Arthropod decline in grasslands and forests is associated with landscape-level drivers[END_REF][START_REF] Tsiafouli | Intensive agriculture reduces soil biodiversity across Europe[END_REF]. Earthworms have suffered a drastic drop in both abundance and diversity: sustainably cropped soils (maize/wheat rotation) support ten times less biomass than grassland soils, and are less functionally diverse due to the scarcity of epigeic species and a higher proportion of immature individuals [START_REF] Curry | Intensive cultivation can drastically reduce earthworm populations in arable land[END_REF][START_REF] Riley | Soil structure, organic matter and earthworm activity in a comparison of cropping systems with contrasting tillage, rotations, fertilizer levels and manure use[END_REF]. These changes in soil macrofauna have long been attributed to tillage and soil compaction by agricultural machinery and to larger plots, but recent ecotoxicological studies strongly suggest harmful sublethal effects of pesticide mixtures circulating in soils at the (infra)individual, the populational and the community level [START_REF] Pelosi | Pesticides and earthworms. A review[END_REF][START_REF] Pelosi | Residues of currently used ACCEPTED MANUSCRIPT / CLEAN COPY pesticides in soils and earthworms: A silent threat? Agriculture[END_REF]. Although earthworm biodiversity in agroecosystems is declining [START_REF] Smith | Diversity and abundance of earthworms across an agricultural land-use intensity gradient[END_REF], some species persist despite the harsh conditions caused by soil disturbance and pesticide applications. This is particularly true of endogeic species such as genera Aporrectodea and Allolobophora sp (Givaudan, Binet, et al., 2014;Givaudan, Wiegand, et al., 2014;[START_REF] Pelosi | Reduction of pesticide use can increase earthworm populations in wheat crops in a European temperate region[END_REF]. Certain species may thus have developed tolerance mechanisms to pesticide exposure. Excessive concentrations of pollutants can be tolerated via two mechanisms: either by plastic physiological acclimation or by genetic adaptation over several generations [START_REF] Morgan | Microevolution and ecotoxicology of metals in invertebrates[END_REF]. The dependence of agriculture on the extensive use of pesticides may now cover a sufficiently long period (~70 years) that adaptations may have occurred in different species, the mechanisms of which are not yet fully understood. This is particularly so for earthworms, in contrast to insects [START_REF] Bass | The global status of insect resistance to neonicotinoid insecticides[END_REF][START_REF] Nuutinen | Comparing the organophosphorus and carbamate insecticide resistance mutations in cholin-and carboxyl-esterases[END_REF]. The evolutionary impacts ACCEPTED MANUSCRIPT / CLEAN COPY of contamination by pesticides on non-target species are rarely documented and therefore poorly understood, yet they may have important consequences on ecosystem functions. Agricultural soils can be managed in different ways, the major differences including tillage and the application of pesticides. Some crops are grown according to the principles of organic farming and others according to conventional agriculture, pesticide application and concentration varying greatly between and within each production method. Depending on their location, earthworm populations therefore have a history of contrasting exposures to chemicals. Because earthworm populations from plots farmed conventionally evolve in a contaminated environment, they could respond to residual contamination by pesticides through acclimation and/or adaptation mechanisms. Earthworm acclimation/adaptation processes have been evidenced so far for soil contamination by metals [START_REF] Fisker | Genetic adaptation of earthworms to copper pollution: Is adaptation associated with fitness costs in Dendrobaena octaedra?[END_REF][START_REF] Posthuma | Heavy-metal adaptation in terrestrial invertebrates: A review of occurrence, genetics, physiology and ecological consequences[END_REF]. The main question is whether the soil fauna (earthworms) can cope and develop tolerance mechanisms to chronic multi-residual contamination by pesticides. To date and to our knowledge, the process of acclimation in natura to long-term pesticide exposure of earthworms has only been observed in the species A. caliginosa (Givaudan, Binet, et al., 2014;Givaudan, Wiegand, et al., 2014). It has been shown that earthworms from soil managed in a conventional farming framework have improved their biotransformation and antioxidant response capabilities. Therefore, there is an urgent need to understand and predict in natura the sublethal effects and evolutionary consequences of agricultural pesticides on soil biodiversity. We continued to study the endogeic species Aporrectodea caliginosa as a model earthworm because it is dominant in most European temperate agrosystems [START_REF] Boag | Distribution, prevalence and intensity of earthworm populations in arable land and grassland in Scotland[END_REF][START_REF] Curry | Relationships between earthworm populations and management intensity in cattle-grazed pastures in Ireland[END_REF][START_REF] Jordan | Effect of management practices and cropping systems on earthworm abundance and microbial activity in Sanborn Field: a 115-year-old agricultural field[END_REF][START_REF] Lamandé | Changes of pore morphology, infiltration and earthworm community in a loamy soil under different agricultural managements[END_REF][START_REF] Nuutinen | Comparing the organophosphorus and carbamate insecticide resistance mutations in cholin-and carboxyl-esterases[END_REF][START_REF] Pérez-Losada | Phylogenetic assessment of the earthworm Aporrectodea caliginosa species complex (Oligochaeta: Lumbricidae) based on mitochondrial and nuclear DNA sequences[END_REF]. A. caliginosa is thus representative of cultivated fields in temperate regions and has recently been identified as a relevant model species for research in soil ecotoxicology [START_REF] Bart | Aporrectodea caliginosa, a relevant earthworm species for a posteriori pesticide risk assessment: current knowledge and recommendations for culture and experimental design[END_REF][START_REF] Pelosi | Searching for a more sensitive earthworm species to be used in pesticide homologation tests-a meta-analysis[END_REF]. To document any evolutionary effects of pesticides, the current study examined two working hypotheses: H1) Populations of A. caliginosa earthworms from conventionally managed fields have developed a tolerance to pesticides compared with those from organically managed fields; H2) This tolerance to pesticides is expressed by a faster and/or more intense response of the molecular pathways involved in biotransformation/detoxification, oxidation and genral response to stress. In order to investigate these hypotheses, we sampled A. caliginosa earthworms from sites with either low-level contamination (organic) or high-level contamination (conventional) to detect any constitutive expression of differentially expressed molecular pathways among these populations. These populations (F0) were then experimentally exposed in the laboratory to pesticide stress to identify differences in molecular responses when again exposed to pesticide (in this case epoxiconazole, a persistent fungicide frequently applied in European agriculture to protect cereals and sugar beet against Septoria sp. and rust). State-of-the art omics technology (RNA sequencing) and bioinformatics were used to characterize the molecular mechanisms of tolerance in a non-targeted way. Overall, our study sheds new light on the molecular mechanisms of tolerance to chemical pollution in soil of terrestrial invertebrates and helps to clarify how species are able to cope with pesticide contamination in agrosystems.

ACCEPTED MANUSCRIPT / CLEAN COPY

Materials and Methods

Study area, agricultural context and earthworm populations

This study took place in agricultural landscapes that are representative of the Atlantic biogeographical region of Europe with a temperate climate (mean rainfall of 696 mm year -1 , mild annual mean temperature of 12.1°C). Agricultural land management is focused on intensive cropping dedicated to dairy and livestock productions and is mainly conventional, based on systematic use of chemicals (mineral fertilizers, pesticides) but with an increasing trend toward organic management.

Earthworm A. caliginosa, an endogeic species distributed worldwide in temperate areas, was chosen as an ecologically-relevant biological model [START_REF] Klobučar | Aporrectodea caliginosa, a suitable earthworm species for field based genotoxicity assessment?[END_REF]. We studied two populations of A. caliginosa, one from a conventionally cropped site (henceforth referred to as "population from the highly contaminated sampling site" or "the HC population") and one from an organically cropped site (henceforth referred to as "population from the lowly contaminated sampling site" or "the LC population") (GPS coordinates: HC sampling site: 48.111180, -1.775199; LC sampling site: 48.127222, -1.725714; approximate coordinates to maintain the farmer's privacy). The cropped sampling sites studied all lie within the same geological basin (silty Basin of Rennes, Brittany, France). Soils are loamy (conventional and organic sites, respectively: Clay 14.0% and 17.9%; Silt 71.2% and 63.2%; Sand 14.8% and 18.9%; organic matter 1.86% and 2.89%; pHw 6.5 and 7.2). The conventionally managed site had been cropped under rotations of corn/cereals (wheat, triticale)/protein crops (peas and fava beans) for at least 22 years, and treated annually with pesticides. The organic site had been under a rotation of grassland (3 years), corn, wheat, rapeseed, and barley without any pesticides for 29 years. Similar tillage and farming practices were applied to both the LC and HC sites. In order to study earthworm tolerance to pesticide pressure in their environment, adults from both the HC and LC populations were reared in another common soil (henceforth termed "rearing soil") whose properties were similar to the LC and HC sites (15.4% clay, 72.2% silt, 12.4% sand, 2.63% organic matter, pHw 6.6) and was chemical-free [i.e, taken from a site that had been managed organically since 1991 under a rotation of grassland (3 years), corn, wheat, rapeseed, and barley (GPS coordinates: 48.132686, -1.718434)].

Characterization of soil and A. caliginosa populations in the LC and HC sites

Soil and earthworm sampling

In 2018, soil and A. caliginosa individuals were collected from the LC and HC sampling sites for general characterization and for the experimental epoxiconazole exposure. For pesticide residue analysis, soil and A. caliginosa individuals were concomitantly sampled at three points 50 m apart. A total of six soil cores per sampling point were taken with an auger at a depth of 20 cm and pooled as a composite soil sample. At least three adult specimens of A. caliginosa were manually collected per sampling point and pooled together for pesticide analysis (once gut-voided). For genetic characterization of A. caliginosa populations, 15 individuals were manually collected from each site.

ACCEPTED MANUSCRIPT / CLEAN COPY

In order to compare HC and LC population abundance, earthworms were sampled using chemical extraction followed by hand-sorting [START_REF] Bouché | Lombriciens de France[END_REF][START_REF] Duriez | An adapted method for sampling earthworms at night in wildlife studies[END_REF]. A diluted expellant solution of allyl isothiocyanate (AITC) was applied on the soil within a 100x100 cm metal frame i.e., 1 m 2 . AITC stock solution 5 g L -1 in isopropanol (propan-2-ol) was further water diluted to a concentration of 0.1 g L -1 [START_REF] Zaborski | Allyl isothiocyanate: an alternative chemical expellant for sampling earthworms[END_REF]. After collecting emerging individuals, a block of soil 31.7x31.7 cm² and 25 cm deep was excavated from where the expellant solution was applied, and the remaining earthworms were hand-sorted. Four sampling replicates, 50 m apart, were performed per site.

Analyses of pesticide residues in HC/LC soil and A. caliginosa

A total of 73 molecules covering herbicides, fungicides, and insecticides were targeted (Table S1) following a survey of farmers on their crop protection practices (i.e., frequency and application doses).

In order to quantify the residues of selected pesticides, a QuEChERS extraction method was applied to soil and earthworm samples, followed by liquid chromatography analysis coupled with tandem mass spectrometry (LC-MS/MS) as described in [START_REF] Daniele | Development of a method for the simultaneous determination of multi-class pesticides in earthworms by liquid chromatography coupled to tandem electrospray mass spectrometry[END_REF].

Genetic characterization of Aporrectodea caliginosa populations

The A. caliginosa species complex contains at least three cryptic lineages/species (L1, L2 and L3) that can be distinguished using a DNA barcoding approach [START_REF] Shekhovtsov | Different dispersal histories of lineages of the earthworm Aporrectodea caliginosa (Lumbricidae, Annelida) in the Palearctic[END_REF]. The cytochrome c oxidase subunit I mitochondrial gene (COI) fragment was used to determine to which lineage belonged the two studied populations. The COI is proposed as a standard DNA barcode for animals [START_REF] Hebert | Biological identifications through DNA barcodes[END_REF]. The fragment was amplified according to Folmer et al., 1994 and details are given in supporting information S1.

Onsite population densities of A. caliginosa

After being collected at the sampling sites, A. caliginosa individuals were weighed back at the laboratory without emptying their gut content (fresh weight basis), counted and taxonomically identified based on morphological criteria [START_REF] Bouché | Lombriciens de France[END_REF]. The abundance and biomass of A. caliginosa populations were calculated for both the HC and LC sites.

Experimental epoxiconazole exposure

After being collected from the HC and LC sites, A. caliginosa earthworms were acclimatized for 14 days in the common rearing soil in a climate chamber (air temperature: 15°C; light day/night cycle: 16/8 h; air moisture: 80±5%). Upon retrieval, the 25 cm top soil dedicated to earthworm rearing was airdried until it reached 14% of WHC then sieved (2 mm sieve) and kept in sealed containers (100 l) prior to use. The rearing soil was then subjected to chemical analyses for pesticides (300 molecules, GIRPA lab Beaucouzé, France) and for polycyclic aromatic hydrocarbons (PAH) and metal trace elements by ACCEPTED MANUSCRIPT / CLEAN COPY the INRAE soil analysis laboratory (Arras, France). These facilities are certified by the French Ministry of the Environment and have a COFRAC accreditation. The condition that the soil had to meet in order to be used as a rearing soil for our experiment was the absence or minimal contamination by these pollutants (see Table S2).

Prior to the start of the exposure, each A. caliginosa individual was rinsed with distilled water, gently dried on filter paper and weighed. Ten adults were then transferred into each mesocosm, insuring a similar mean earthworm weight at each treatment. The sampling design is detailed in Table S3. Soil mesocosms consisted of polypropylene boxes (175 mm x 159 mm x 117 mm) whose lid was pierced with tiny holes to ensure sufficient aeration. The mesocosms were filled with 1 kg of either epoxiconazole spiked or control soil, and 10 g of dry horse manure was added onto the soil surface for earthworm feeding. Epoxiconazole is a persistent fungicide that was chosen due to 1) its high frequency of application in European agriculture (V. [START_REF] Silva | Pesticide residues in European agricultural soils -A hidden reality unfolded[END_REF] and 2) its detection in the highly contaminated site under study (HC). A. caliginosa earthworms were exposed to epoxiconazole as BASF's "OPUS®" commercial formulation (125 g active ingredient 1 -1 ), a broad-spectrum contact and systemic fungicide to protect cereals and sugar beet against Septoria sp. and rust. OPUS® was diluted in distilled water at the recommended dose for farmers of 0.17 mg kg -1 , considering a field application rate of 125 g ha -1 and assuming a single application with a homogenous distribution and no crop interception in the soil's top 5 cm [START_REF] Dittbrenner | Physiological and behavioural effects of imidacloprid on two ecologically relevant earthworm species (Lumbricus terrestris and Aporrectodea caliginosa)[END_REF]. Soil was spiked by manually adding the diluted pesticide solution or distilled water (for the controls) on the soil with a 14% water content until reaching a final soil water content of 25%, and by renewing it after 14 days to sustain the chemical pressure. Exposure lasted 2 days for the transcriptomic responses and 28 days in total for all other endpoints. In the following text, the term "basal" refers to control individuals or conditions, and the term "stress" refers to epoxiconazole-exposed individuals or conditions.

Epoxiconazole analysis in experimental soil and A. caliginosa

Soil samples and gut-voided earthworms were separately freeze-dried, ground, sieved to 250 µm and then stored at -18°C until epoxiconazole analysis. Epoxiconazole was extracted using a pressurized liquid ASE (Accelerated Solvent Extractor) 350 (Dionex Corporation, Sunnyvale, USA) followed by GC/MS/MS analysis according to [START_REF] Mercier | A multi-residue method for the simultaneous analysis in indoor dust of several classes of semi-volatile organic compounds by pressurized liquid extraction and gas chromatography/tandem mass spectrometry[END_REF] with specific adjustments for both soil and earthworm tissues, details are given in Supporting information S2. To assess the earthworms' epoxiconazole accumulation, bioaccumulation factors (BAFs) defined as the ratio of total earthworm epoxiconazole concentration to total soil epoxiconazole concentration at T0 were calculated. Both concentrations were based on dry weight.

ACCEPTED MANUSCRIPT / CLEAN COPY

A. caliginosa transcriptomic responses

At T0 and after two days of exposure to epoxiconazole, A. caliginosa earthworms were individually frozen and ground (whole body, gut voided) with a CryoMill (Retsch) under liquid nitrogen conditions.

Frozen earthworm powders were stored at -80°C until use. Aliquots containing 30 mg of individual earthworm powder were used for total RNA extraction using a NucleoSpin RNA kit (Macherey-Nagel) according to the manufacturer's protocol. Genomic DNA was removed by DNase digestion using a column then total RNA was eluted in RNase-free H2O. RNA purity and quantity were assessed using a Nanodrop (ND-1000, ThermoFisher), and RNA integrity was checked using a Bioanalyzer 2100 (Agilent, CA, USA). Total RNA samples were stored at -80°C. Samples were sent to GeT-PlaGe core facility (INRA Toulouse) in dry ice. RNA sequencing, de novo assembly, annotation, quantification, differential expression analysis, and enrichment analyses are described in supporting information S3.

Analysis of A. caliginosa metabolism

The weight of 30 gut-voided earthworms from each of the control and epoxiconazole exposure groups was individually monitored during the experiment at T0 and after 7 and 28 days of exposure.

Respirometry measurements were taken on T0 and after 7 and 28 days of exposure. Ten earthworms from each treatment group were removed from the soil microcosm, rinsed, gently blotted dry on filter paper, and placed in a 250 ml glass jar hermetically closed for two hours. CO2 was measured by a microgas chromatograph (3000A, SRA Instruments) equipped with a single PoraPLOT U capillary column coupled with a thermal conductivity detector.

Statistical analyses

Statistical analyses were carried out using the statistical software R (RStudio, R 4.0.3). Details on transcriptomics analysis are given in supporting information S3. For epoxiconazole analysis in soil and earthworms, weight monitoring, and respirometry measurements, all data were checked for normality distribution (Shapiro-Wilk test) and homogeneity of variances (Bartlett's test), with visual examination of QQ-plots. When assumptions were met, a one-way ANOVA was run with Tukey's post hoc tests.

The non-parametric Kruskal-Wallis test was used if assumptions were not met; treatment groups were compared using a Dunn's pairwise comparison with Bonferroni correction. All tests were performed with a 0.05 significance level.

Results

Characterization of A. caliginosa populations from the sampling sites

Aporrectodea caliginosa COI genetic diversity

A total of nine haplotypes (H1-H9) were obtained for the 30 individuals sampled. All sequences were assigned to Aporrectodea caliginosa (match score > 99%) and several haplotypes (H2, H4, H8 and H9) ACCEPTED MANUSCRIPT / CLEAN COPY matched A. caliginosa L3 in The Barcode of Life Data System (BOLD). The maximum likelihood tree confirmed the BOLD assignment and further revealed that all study haplotypes belonged to the same well-supported clade formed by the A. caliginosa L3 lineage (Figure 1). A higher haplotype diversity was observed in the LC population whereas a higher nucleotide diversity was observed in the HC population (Figure 1). This higher nucleotide diversity is due to the presence in the HC population of two divergent haplotypes (H1 and H7, Figure 1) that are absent in the LC population.

ACCEPTED MANUSCRIPT / CLEAN COPY

Onsite abundance and biomass of A. caliginosa populations

The mean abundance of A. caliginosa was more than three times higher in the LC site than the HC one (56.5 ± 28.1 and 15.3 ± 7.0 individuals.m -2 in LC and HC, respectively). The corresponding mean A. caliginosa biomasses were 18.5 ± 9.8 and 3.0 ± 0.9 g.m -2 in the LC and HC sites, respectively.

Onsite soil and earthworm pesticide contamination

The multi-residual analysis showed that all samples contained at least one and up to 20 and 11 pesticides for soil and earthworms, respectively. Out of the 73 targeted pesticides, 28 and 15 had been detected at least once in the soil and A. caliginosa, respectively. Soil and A. caliginosa from the HC site contained a higher number and higher concentrations of pesticides than the LC site. The highest concentrations in soil from the HC site were for the herbicides aclonifen, pendimethalin and clomazone, with 417.0, 15.6, and 113.0 ng/g dw. The highest concentrations of pesticides in earthworms from the same site were for metolachlor ESA (herbicide) and imidacloprid (insecticide) with 747.0 and 99.2 ng/g ww, respectively).

In both soil and A. caliginosa earthworms from the LC site, only five (atrazine, carbendazim, flusilazole, metconazole, fluxapyroxad) and three molecules (imidacloprid, tebuconazole and metolachlor ESA), respectively, were detected and all were at concentrations below 5.0 ng/g ww (Table S4).

Establishing the de novo transcriptome assembly of earthworm species A. caliginosa

A total of 2,068,522,212 raw reads were generated, with an average of 23,729,267 (± 1,892,470) reads for each sample. Of these, a total of 139,950,716 reads were used for the transcriptome assemblies. The clean reads were assembled into 64,556 contigs with an average length of 2,010 pb ranging from 203 to 28,628 pb and N50 of 2,862 pb. The assembled transcriptome size was 129,787,458 bp. The assembly quality was then assessed by calculating the read on contig re-alignment rate for every sample, and by processing the contigs with BUSCO (version 2) using the metazoa_odb9 database (Simao et al. 2015). This analysis was to check the presence and completeness of a set of expected single-copy proteincoding genes for a given branch of the evolution tree. Among the 978 metazoa_odb9 proteins searched by BUSCO (3.0.2), 964 were found to be complete in single or multiple copies. This represents 98.6% of the proteins expected in the genome. The remaining ones were missing (1.1%) or fragmented (0.3%).

With respect to the assembly quality, the mean per sample read mapping rate was 97.02 ± 0.60% across libraries. RNA-seq raw sequences and the de novo assembled transcriptome assemblies have been deposited in NCBI under BioProject PRJNA883218. This Transcriptome Shotgun Assembly project has been deposited at DDBJ/EMBL/GenBank under the accession GKBY00000000. The version described in this paper is the first version, GKBY01000000.

ACCEPTED MANUSCRIPT / CLEAN COPY

Functional annotations of the de novo transcriptome assembly

A total of 25,725 contigs could be annotated, with reference to the sequences recorded in the NR protein database. The marine polychaete Capitella teleta (Phylum: Annelida) was the species found most often in the NR protein annotations, corresponding to 27.16% of the contigs (6,988 contigs). The other species found as best hit annotations were the leech Helobdella robusta (4,371 contigs, 16.99%) and the brachiopod Lingula anatina (3,146 contigs, 12.23%) (Figure 2A).

In total, 43,569 GO terms were identified for 21,593 annotated sequences under the three main ontologies. From these sequences, 18,160 (41.7%) were assigned to molecular functions (GO:0003674), 13,782 (31.6%) to biological processes (GO:0008150), and 11,626 (26.7%) to cellular components (GO:0005575). Within the molecular function category, binding (GO:0005488) and catalytic activity (GO:0003824) were the most represented GO terms with 12,956 sequences and 8,763 sequences, respectively. The top three GO terms for biological processes were cellular process (GO:0009987, 9802 contigs) including metabolic process (GO:008152, 9120 contigs), biological regulation (GO:0065007, 4040 contigs) and response to stimulus (GO:0050896, 2810 contigs). Furthermore, in the cellular component category, the predominant GO terms were grouped into cell (GO:0005623, 6785 contigs), ACCEPTED MANUSCRIPT / CLEAN COPY membrane (GO:0016020, 6245 contigs), organelle (GO:0043226, 4242 contigs) and protein-containing complex (GO:0032991, 2186 contigs) (Figure 2B).

In order to identify active functional pathways in A. caliginosa, contigs were mapped to the reference pathways in the KEGG database. In total, 16,933 contigs were mapped to 404 KEGG pathways and 5,847 KEGG orthologs. The KEGG annotations were helpful for identifying contigs related to xenobiotics biodegradation and metabolism-related pathways (Table S5). The pathways with the highest number of sequence hits to enzyme genes were Metabolism of xenobiotics by cytochrome P450 (map00980), Drug metabolismcytochrome P450 (map00982) and Drug metabolismother enzymes (map00983).

Molecular responses of LC and HC populations under basal and stress conditions

Analyses were performed to check the overall reproducibility and variation between earthworm samples belonging to the same population (HC or LC) and treatment (exposed (stress) or not (basal)). A heatmap of Pearson correlations based on log10 expression counts and an MDS plot were generated to evaluate the sample relatedness and identify outliers (Figure 3). The heatmap and MDS plot showed that individuals belonging to the same A. caliginosa population were highly correlated. However, one sample from the LC site was identified as an outlier (T0-LC-4) and was therefore removed from the downstream analysis. Under basal condition (LC.basal versus HC.basal), analysis showed that 939 contigs were differentially expressed with 569 contigs upregulated and 370 downregulated when comparing LC and HC populations that had a contrasting site history with respect to agricultural practices and pesticide use ACCEPTED MANUSCRIPT / CLEAN COPY (Figure 4A). Annotations enabled us to identify 226 proteins (24% of all the differentially expressed contigs) (Figure 4B, Table S6). Results of the gene ontology enrichment analysis for biological process and molecular function are presented in Figures 4D and4E. KEGG enrichment analysis revealed that one pathway was upregulated under the basal condition: protein digestion and absorption (map04974) and three pathways were downregulated: ECM-receptor interaction (map04512), pancreatic secretion (map04972) and protein digestion and absorption (map04974) (Table S8).

Under stress condition (LC.stress versus HC.stress), 761 contigs were differentially expressed by the two populations (321 downregulated and 440 upregulated) (Figure 4A) and 182 contigs were annotated (24%) (Figure 4B, Table S7). Results of the gene ontology enrichment analysis for biological process and molecular function are presented in Figures 4D and4E. KEGG analysis showed that three pathways were specifically downregulated under the stress condition: autophagy (map04136, map04138) and mitophagy (map04139) (Table S8).

ACCEPTED MANUSCRIPT / CLEAN COPY

Testing for differential molecular responses to an induced chemical stress (fungicide/epoxiconazole) in the LC and HC populations

After 2 days of the earthworms' exposure to epoxiconazole, 52 contigs were upregulated and 50 contigs downregulated in the LC population compared with 42 upregulated contigs and 52 downregulated contigs in the HC population.

In the LC population (LC.basal vs. LC.stress), 21% of differentially expressed contigs were annotated.

Epoxiconazole exposure induced both upregulation and downregulation of different immune-related ACCEPTED MANUSCRIPT / CLEAN COPY proteins (neuromacin, Serpin B6, C-type lectin domain family 4 member A, C1q domain-containing protein), and upregulation of the RNA splicing protein (splicing factor 3B subunit 5), stress response protein (Heat shock protein 90), and hydrolytic enzymes such as chitinase. Epoxiconazole exposure also induced downregulation of stress oxidative responses such as superoxide dismutase, and the sperm axoneme assembly protein (axoneme-associated protein mst101(2)-like) in the exposed group (Figure 5C). In the HC population (HC.basal vs. HC.stress), 13% of differentially expressed contigs were annotated. Expoxiconazole exposure induced upregulation of the chitin binding-related protein, cytochrome P450 family protein (CYP2 family), immune-related proteins (C1q domain-containing protein, lymphocyte cytosolic protein 2-like isoform X2), and the blood glycoprotein (VWFA domaincontaining protein). Exposure also induced downregulation of a protein involved in reduction-oxidation mechanisms (tyrosinase-like protein 1), chemoattractive glycoprotein ES20, an RNA biosynthetic process-related protein (c-Myc-binding protein like), the vesicle transport protein, and immune response protein (C1q domain-containing protein, putative deleted in malignant brain tumors 1 protein) (Figure 5C). Overall, both populations shared only six contigs (Figure 5A) in response to epoxiconazole exposure, just two of these being annotated (Figure 5B) (C1q domain-containing protein), and one BPs (chitin metabolic process) was enriched in both populations.

In order to understand the function of the differentially expressed contigs in the control and exposed groups in both the LC and HC populations, GO enrichment analysis was performed. In the LC population, epoxiconazole caused upregulation of 16 BPs linked to four main functions: cellular response to heat (GO:0034605), protein stabilization (GO:0050821), RNA splicing (GO:0000398), and the chitin metabolic process (GO:0006030). Moreover, four MF terms were upregulated in the exposed group linked to protein binding and chitin binding (Table S9).

In the HC population, exposure to epoxiconazole downregulated 10 BP GO terms linked to positive regulation of RNA biosynthetic process (GO:1902680) and endocytosis (GO:0006897) while 11 BP terms were upregulated, most being linked to the chitin metabolic process (GO:0006030), and response to xenobiotic stimulus (GO:0009410). Four MF GO terms downregulated in the exposed group were linked to scavenger receptor activity (GO:0005044) and transcription coactivator activity (GO:0003713), while 3 MF GO terms were upregulated; these were linked to extracellular matrix structural constituent (GO:0005201), chitin binding (GO:0008061) and oxidoreductase activity (GO:0016712) (Table S9).

ACCEPTED MANUSCRIPT / CLEAN COPY

Bioaccumulation and physiological responses after chemical stress

Epoxiconazole analysis in soil and A. caliginosa earthworms

Our analysis revealed that epoxiconazole concentrations in soil matched nominal concentrations on T0, except for a slightly lower concentration in soil containing earthworms from the HC site, probably due to heterogeneity in pesticide application. The fate of epoxiconazole during the experiment did not differ significantly between the two populations, with a similar decrease in concentration of 2% vs. 8% (day 7) and 11% vs. 15% (day 28) in the HC and LC mesocosms, respectively. Analyses of epoxiconazole in earthworm tissues showed a similar bioaccumulation factor (BAF around 2) on day 7 in the two populations. Lower bioaccumulation was observed on day 28 in the HC population, although not significant (BAF of 2 and 1.4 for LC and HC respectively) (Figure 6A).

A. caliginosa weight monitoring throughout the experiment

The initial mean weight of adult A. caliginosa earthworms (n=30) collected from the HC and LC sites was 0.562 ± 0.116 g and 0.556 ± 0.099 g, respectively. Survival was 100% in all the mesocosms.

ACCEPTED MANUSCRIPT / CLEAN COPY

Throughout the experiment, all earthworms gained weight similarly except those of the HC control that slightly decrease their initial weight. On day 28, a significant weight increase was observed in the HC exposed group (p=0.0355) compared to HC control group's weight (Figure 6B).

A. caliginosa metabolism monitoring throughout the experiment

At the beginning of the experiment (T0), earthworm respiration measured as CO2 production was higher in the LC population than in the HC one, although not significantly due to the high variability between individual earthworms. The metabolic rate of the HC populations increased significantly over time, whether exposed or not to expoxiconazole. Conversely, an unsignificant variation in energy dissipation was observed in the LC earthworms throughout the experiment (p>0.05) (Figure 6C).

ACCEPTED MANUSCRIPT / CLEAN COPY

Discussion

A. caliginosa facing pesticide soil contamination

Conventional farming leads to a residual stock of a cocktail of pesticides in soil; up to 20 different chemicals were detected in the soil of the HC site. Consequently, earthworms bioaccumulated a high level of pesticides in their tissue (up to 987.9 ng g -1 wet weight). Most were insecticides (neonicotinoids: imidacloprid, thiacloprid) or herbicides (mesotrione, metolachlor ESA). Very few studies have investigated the bioaccumulation of pesticides in the tissue of earthworms from arable land, though many have focused on metal bioaccumulation [START_REF] Beaumelle | Relationships between metal compartmentalization and biomarkers in earthworms exposed to fieldcontaminated soils[END_REF]. Recently, Pelosi et al., 2021 also found high concentrations of currently used pesticides in earthworms, within a range similar to our own observations. The population-level consequences of such chemical pressures can include different responses, such as a decline leading to future extinction with the loss of the associated ecological services, or survival mediated by either physiological or behavioral plasticity, or by adaptation through natural selection [START_REF] Chevin | When do adaptive plasticity and genetic evolution prevent extinction of a density-regulated population?[END_REF].

Genetic diversity among the studied populations

The differential tolerance observed in some population samples from reference and contaminated sites may represent a loss of species diversity through local contaminant-caused extinction of one or more members of a cryptic species complex rather than a within-species loss of putatively less-tolerant genotypes or haplotypes [START_REF] Rocha-Olivares | Differential tolerance among cryptic species: A potential cause of pollutant-related reductions in genetic diversity[END_REF]. Recent genotyping studies have begun to challenge current understanding of diversity through the identification of genetically distinct cryptic lineages within previously established earthworm morphospecies. There is some evidence that different earthworm lineages have different ecological properties, as observed in numerous other invertebrate species [START_REF] Spurgeon | Ecological drivers influence the distributions of two cryptic lineages in an earthworm morphospecies[END_REF]. [START_REF] Kille | DNA sequence variation and methylation in an arsenic tolerant earthworm population[END_REF] demonstrated that two cryptic lineages of earthworm L. rubellus react differently to soil contamination by arsenic, the specimens from lineage A employing genetic adaptation while the lineage B specimens use an epigenetic strategy. Although rarely performed in ecotoxicological study, checking for cryptic species in earthworms should thus be mandatory. For instance, A. caliginosa-considered a relevant model species for ecotoxicology research-is a complex of three highly divergent genetic lineages [START_REF] Fernández | Adding complexity to the complex: New insights into the phylogeny, diversification and origin of parthenogenesis in the Aporrectodea caliginosa species complex (Oligochaeta, Lumbricidae)[END_REF][START_REF] Pérez-Losada | Phylogenetic assessment of the earthworm Aporrectodea caliginosa species complex (Oligochaeta: Lumbricidae) based on mitochondrial and nuclear DNA sequences[END_REF][START_REF] Porco | Biological invasions in soil: DNA barcoding as a monitoring tool in a multiple taxa survey targeting European earthworms and springtails in North America[END_REF]. The nature of these lineages is not completely understood: The number of nucleotide substitutions between them in the mitochondrial COI gene is much higher than within the lineages [START_REF] Porco | Biological invasions in soil: DNA barcoding as a monitoring tool in a multiple taxa survey targeting European earthworms and springtails in North America[END_REF], and based on these criteria, the lineages may represent cryptic species. In this study, we showed that all the specimens of both populations belonged to the same L3 lineage, which is known to be present both in Western Europe (native populations) and North America (introduced populations). Thus, we confirmed using COI DNA barcoding that our ecotoxicological experiment was performed on a single species. As tolerance to environmental pollution in a population with a complex genetic structure can be achieved by various means, the COI data was also used to analyze the genetic variation ACCEPTED MANUSCRIPT / CLEAN COPY of the populations under study, keeping in mind that information obtained using COI mtDNA haplotypes is limited because only one locus is screened (Lazrek et al., 2011). We showed that the HC population had a higher nucleotide diversity than the LC population, a result mainly explained by the presence in this population of two divergent haplotypes that were absent in the LC population. Genetic diversity in earthworms can be affected by many factors, such as ecological preferences or life history traits, but also pollution, geographic distance, and geographic barriers for example [START_REF] Dupont | Genetic variation of the epigeic earthworm Lumbricus castaneus populations in urban soils of the Paris region (France) revealed using eight ACCEPTED MANUSCRIPT / CLEAN COPY newly developed microsatellite markers[END_REF]. To date, it is still difficult to draw general conclusions about spatial genetic variation in earthworms, especially at the fine scale [START_REF] Dupont | Dispersal constraints and fine-scale spatial genetic structure in two earthworm species[END_REF] as in our study. One of the possible environmental pressures that might modify genetic diversity between our populations could be the soil's contamination by pesticides, all other things being equal (i.e., pedo-climate context, soil properties and nutrients, metallic and HAP pollution). In contaminated environments, significant changes in genotype and allele frequencies may occur as a result of different processes: (i) an increase in mutation rates, (ii) selection of resistant genotypes, (iii) bottlenecks where population sizes are drastically reduced, and (iv) reduction of migration rates, preventing gene exchange among populations [START_REF] Dupont | Genetic variation of the epigeic earthworm Lumbricus castaneus populations in urban soils of the Paris region (France) revealed using eight ACCEPTED MANUSCRIPT / CLEAN COPY newly developed microsatellite markers[END_REF][START_REF] Van Straalen | Genetic variation in toxicantstressed populations: an evaluation of the "genetic erosion" hypothesis[END_REF].

Interpopulation differences: differentially expressed contigs under basal and stress conditions

Gene ontology analysis revealed that terms with a high number of differentially expressed contigs and a low p-value were linked to proteolysis (BP) and peptidase activity (MF) (Figure 4). These terms were mainly upregulated in the HC population rather than the LC population under both basal and stress conditions. Proteases or peptidases play a key role in protein turnover, along with transcription and expression of the active genome [START_REF] Wilkins | Insecticide resistance and intracellular proteases[END_REF]. Elevated proteolytic activities are one of the responses to pesticide exposure or in pesticide-resistant strains of insect pests (see [START_REF] Wilkins | Insecticide resistance and intracellular proteases[END_REF] for a detailed review on the subject). It has even been proposed as a basis for resistance detection [START_REF] Gong | Cloning and overexpression of CYP6F1, a cytochrome P450 gene, from deltamethrin-resistant Culex pipiens pallens[END_REF].

In our study, we can note the upregulation of different proteins linked to the ubiquitin-proteasome system (UPS) (A Chain A, Structure Of Usp15 Dusp-Ubl Deletion Mutant; ubiquitin carboxyl-terminal hydrolase 15 isoform X1). The UPS is a selective proteolytic system in which the conjugation of ubiquitin to damaged or unneeded proteins induces their degradation by the proteasome. Through overexpression of the proteasome Subunit Beta Type 6 (PSMB6), mosquitoes had acquired stable resistance to the insecticide deltamethrin (L. [START_REF] Sun | Identification of proteasome subunit beta type 6 (PSMB6) associated with deltamethrin resistance in mosquitoes by proteomic and bioassay analyses[END_REF].

The HC population also showed high lumbrokinase P2 enzyme upregulation. Earthworm fibrinolytic enzyme (EFE) has similar characteristics to alkaline trypsin-, chymotrypsin-and elastinase-like serine proteases with a strong fibrinolytic activity in the digestive tract of earthworms [START_REF] Zhao | Earthworm fibrinolytic enzyme[END_REF].

Trypsin and chymotrypsin serine proteases were overexpressed in a deltamethrin-resistant strain of mosquito Culex pipiens pallens ensuring a greater chance of survival under pesticide pressure [START_REF] Wu | Culex pipiens pallens: identification of genes differentially expressed in deltamethrin-resistant and-susceptible strains[END_REF]. We can note upregulation of cathepsin L1, mitogen-activated protein kinase (MAPK), and ACCEPTED MANUSCRIPT / CLEAN COPY matrix metalloproteinase-18-like in the HC population, which were also implicated in responses or resistance to pesticides [START_REF] Wilkins | Insecticide resistance and intracellular proteases[END_REF]. Different explanations have been proposed for increased proteolytic activity linked to pesticide resistance or exposure (A. [START_REF] Silva | Insecticide Resistance Mechanisms in the Green Peach Aphid Myzus persicae (Hemiptera: Aphididae) I: A Transcriptomic Survey[END_REF]. Peptidases may be involved in protein degradation to interact with damaged protein or to fulfill higher energy demands, which is usually a response to stress [START_REF] Pedra | Genome-wide transcription profile of field-and laboratory-selected dichlorodiphenyltrichloroethane (DDT)-resistant Drosophila[END_REF]. Peptidases may also play a role during protein biosynthesis or in modification of enzyme conformation related, for example, to the metabolic machinery required to detoxify insecticides [START_REF] Ahmed | Comparison of proteolytic enzyme activities in adults of insecticide resistant and susceptible strains of the housefly M. domestica L[END_REF]. In our study this elevated peptidase activity in the HC population could be due to: 1) an actual response to onsite pesticide exposure of HC earthworms (they were exposed to pesticides at the sampling site shortly before being collected, although there was a 14-day acclimation period in the lab in soil without pesticides, see materials and methods section), 2) or conversely an additional stress response due to the absence of pesticides in the rearing soil and/or 3) a response due to acclimation or adaptation mechanisms. Indeed, we can observe that under the stress condition there are fewer differentially expressed contigs linked to proteolysis in the HC population. This could be explained by the fact that earthworms from the HC population were already adapted to the induced stress due to their long history of exposure to chemicals. The metabolic differences in terms of respirometry (CO2 production) and growth performance that we observed between the HC and LC populations support this explanation.

Although the results are noisy, a general trend in differentiated metabolic rate was observed between the LC and HC populations irrespective of experimental treatments. It is also interesting to note that earthworms from the HC population were significantly smaller (mean weight of 0.196 g) than individuals from the LC population (mean weight of 0.327g) in concordance with the high metabolic activity observed in gene expression results, as already observed in other studies suggesting energetic costs in organisms coping with soil pollution [START_REF] Fisker | Genetic adaptation of earthworms to copper pollution: Is adaptation associated with fitness costs in Dendrobaena octaedra?[END_REF][START_REF] Holmstrup | Body metal concentrations and glycogen reserves in earthworms (Dendrobaena octaedra) from contaminated and uncontaminated forest soil[END_REF][START_REF] Wiegand | Bioaccumulation of paraquat by Lumbriculus variegatus in the presence of dissolved natural organic matter and impact on energy costs, biotransformation and antioxidative enzymes[END_REF]. These results are in accordance with (Givaudan, Wiegand, et al., 2014), who showed that the metabolic rate increased after fungicide exposure in earthworms naturally pre-exposed to pesticides.

Other proteins linked to proteolysis were upregulated in the HC population; serine/threonine kinase Unc-51-like kinase-1 (Ulk1) is thought to be essential for inducing autophagy, a protective intracellular bulk degradation mechanism [START_REF] Li | Serine/threonine kinase Unc-51-like kinase-1 (Ulk1) phosphorylates the co-chaperone cell division cycle protein 37 (Cdc37) and thereby disrupts the stability of Cdc37 client proteins[END_REF] that is activated by various high stresses in cellular response to reactive oxygen and nitrogen species as well as toxic proteins [START_REF] Li | Serine/threonine kinase Unc-51-like kinase-1 (Ulk1) phosphorylates the co-chaperone cell division cycle protein 37 (Cdc37) and thereby disrupts the stability of Cdc37 client proteins[END_REF]. The earliest identified degradation route for proteins is within lysosomes; foreign proteins and cell materials can be transferred to the lysosome by a number of mechanisms which include autophagy [START_REF] Wilkins | Insecticide resistance and intracellular proteases[END_REF].

These results are in accordance with KEGG analysis showing upregulation of pathways linked to these mechanisms in the HC population under both basal and stress conditions (lysine degradation, protein digestion and absorption, pancreatic secretion, autophagy, mitophagy).

Another interesting result observed regarding response linked to metabolism is the higher upregulation of three mitochondrial proteins in the HC population (NADH dehydrogenase subunit 5 (or NADHubiquinone oxidoreductase chain 5), complex I-B15 (NADH-ubiquinone oxidoreductase B15 subunit), ACCEPTED MANUSCRIPT / CLEAN COPY and cytochrome c oxidase subunit I (COI). At the cell level, mitochondria perform a variety of biochemical processes, but their main function is to produce most of the cellular ATP via oxidative phosphorylation [START_REF] Stoker | Impact of pharmacological agents on mitochondrial function: A growing opportunity[END_REF]. Mitochondria are known to be impacted by pesticide exposure, even by non-conventional mitochondria-targeting pesticides [START_REF] Cowie | Molecular networks related to the immune system and mitochondria are targets for the ACCEPTED MANUSCRIPT / CLEAN COPY pesticide dieldrin in the zebrafish (Danio rerio) central nervous system[END_REF][START_REF] Leung | Mitochondria as a target of organophosphate and carbamate pesticides: Revisiting common mechanisms of action with new approach methodologies[END_REF][START_REF] Nicodemo | Fipronil and imidacloprid reduce honeybee mitochondrial activity[END_REF][START_REF] Van Pottelberge | Resistance mechanisms to mitochondrial electron transport inhibitors in a field-collected strain of Tetranychus urticae Koch (Acari: Tetranychidae)[END_REF]. Inhibition of electron transport at the mitochondrial respiratory chain has been a successful mode of action of pesticides [START_REF] Lümmen | Mitochondrial electron transport complexes as biochemical target sites for insecticides and Acaricids[END_REF][START_REF] Van Leeuwen | The molecular evolution of xenobiotic metabolism and resistance in chelicerate mites[END_REF]. A number of mitochondrial respiratory chain genes are associated with pesticide resistance and specifically insecticide resistance, including those mentioned above [START_REF] Van Leeuwen | Mitochondrial heteroplasmy and the evolution of insecticide resistance: Non-Mendelian inheritance in action[END_REF]. NADH dehydrogenase subunit 5 (or NADH-ubiquinone oxidoreductase chain 5) and complex I-B15 (NADH-ubiquinone oxidoreductase B15 subunit) are core and accessory subunits of the respiratory chain Complex I, respectively. Recently, [START_REF] Bajda | A mutation in the PSST homologue of complex I (NADH: ubiquinone oxidoreductase) from Tetranychus urticae is associated with resistance to METI acaricides[END_REF] reported the discovery of a mutation (H92R) in the PSST homolog of complex I in METI-I (mitochondrial electron transport inhibitor) of resistant phytophagous Tetranychus urticae mites. COI is one of three mitochondrial DNA encoded subunits of respiratory complex IV. An elevated expression of mitochondria cytochrome c oxidase as observed in the HC population was involved in the development of resistance to pesticides in several organisms: Chinese hamster ovary cells [START_REF] Alemany | Identification by RNA-based arbitrarily primed PCR of the involvement of cytochrome c oxidase in the development of resistance to methotrexate[END_REF], Blattella germanica [START_REF] Pridgeon | Overexpression of the cytochrome c oxidase subunit I gene associated with a pyrethroid resistant strain of German cockroaches, Blattella germanica (L.)[END_REF], A. aegypti [START_REF] Pridgeon | Permethrin induces overexpression of cytochrome c oxidase subunit 3 in Aedes aegypti[END_REF], Schistosoma mansoni, [START_REF] Pereira | Alterations in cytochrome-c oxidase expression between praziquantel-resistant and susceptible strains of Schistosoma mansoni[END_REF], and in whiteflies [START_REF] Yang | Transcriptomic and Proteomic Responses of Sweetpotato Whitefly, Bemisia tabaci, to Thiamethoxam[END_REF]. It is interesting to note that this COI gene was used to analyze the genetic diversity of the LC and HC populations and showed that two divergent haplotypes were present in the HC population.

In our study, the significant upregulation of some mitochondrial genes may be a compensatory mechanism for electron transport in the respiratory chain through a plastic response or a genetic adaptation. Higher mutation rates in animal mtDNA together with limited DNA repair mechanisms render these genes susceptible to rapid evolution through random drift or natural selection. Traits encoded by mtDNA have thus the potential to evolve, and reach fixation, very rapidly [START_REF] Van Leeuwen | Mitochondrial heteroplasmy and the evolution of insecticide resistance: Non-Mendelian inheritance in action[END_REF]. It is also interesting to note that multiple lines of evidence have recently indicated that the cytosolic (Ubiquitin-proteasome system) UPS plays a crucial role in the quality control of mitochondrial proteins. The effects of the UPS go beyond the removal of damaged proteins and include the adjustment of mitochondrial proteome composition, the regulation of organelle dynamics, and the protection of cellular homeostasis against mitochondrial failure [START_REF] Bragoszewski | Control of mitochondrial biogenesis and function by the ubiquitin-proteasome system[END_REF]. However, the correlation between mitochondrial dysfunction and UPS dysfunction after pesticide exposure needs further investigation.

Interestingly, a shift regarding pathways linked to ion binding (metal ion binding, calcium ion binding) was observed; these pathways were upregulated in the LC population (downregulated in the HC population) under the basal condition and became upregulated in the HC population (downregulated in the LC population) under the stress condition. Further investigations should be performed to understand ACCEPTED MANUSCRIPT / CLEAN COPY this response, but it is well known that ion channels remain the primary target of most of the small molecule insecticides, and are involved in pesticide resistance [START_REF] French-Constant | The molecular and population genetics of cyclodiene insecticide resistance[END_REF]).

In the HC population, upregulation of CytochromeP450 CYP2J2 was observed under the stress condition. In humans, this cytochrome is involved in phase I xenobiotics metabolism [START_REF] Xu | Cytochrome P450 2J2: distribution, function, regulation, genetic polymorphisms and clinical significance[END_REF]. This response was only observed in the HC population and could be a specific mechanism in response to pesticide exposure developed in this population.

Both populations expressed different responses to additional stress

It was also observed that the LC and HC populations displayed almost the same number of differentially expressed contigs in response to epoxiconazole exposure, but they only shared six contigs, indicating a specific response for each population. Unfortunately, the lack of annotation did not allow us to draw a particular pathway in response to exposure. This global result could be due to the fact that the response was analyzed only 2 days after induced epoxiconazole exposure, which is perhaps too short to observe a differentiated response from the two populations. A higher number of differential gene expressions in fungicide-exposed arthropods was observed after 4 days of exposure, compared with 2 and 7 days [START_REF] Simões | Using time-lapse omics correlations to integrate toxicological pathways of a formulated fungicide in a soil invertebrate[END_REF]. Another possible explanation could be that the choice of epoxiconazole was maybe not the best pesticide to trigger a clear response in earthworms as evidenced in other studies:

Epoxiconazole induced either no or moderate effects in earthworms at an environmentally relevant exposure [START_REF] Bart | Differences in sensitivity between earthworms and enchytraeids exposed to two commercial fungicides[END_REF]Givaudan, Binet, et al., 2014;Givaudan, Wiegand, et al., 2014;[START_REF] Nelieu | Evaluation of epoxiconazole bioavailability in soil to the earthworm Aporrectodea icterica[END_REF][START_REF] Pelosi | Sublethal effects of epoxiconazole on the earthworm Aporrectodea icterica[END_REF]).

In the LC population, processes involved in response to stress included upregulation of heat shock protein 90. HSPs are known to be rapidly synthesized in response to environmental stressors, including pesticide stress (Y. [START_REF] Sun | Characterizing heat shock protein 90 gene of Apolygus lucorum (Meyer-Dür) and its expression in response to different temperature and pesticide stresses[END_REF]. Different proteins were also up-or downregulated in response to epoxiconazole without a clear explanation of the exact nature of these changes, e.g., downregulation of the neuromacin propeptide, whose functions include bacterial killing, symbiostasis in the gut, immune defense, and regeneration of the damaged nerve cord [START_REF] Bruno | Worms' antimicrobial peptides[END_REF]. The expression of genes related to these immune processes has been found to change following pesticide exposure [START_REF] Costa | New perspectives on cytokine pathways modulation by pesticide exposure[END_REF]. The BTB domain-containing protein was also downregulated, which is a versatile protein-protein interaction motif that participates in a wide range of cellular functions, including transcriptional regulation, cytoskeleton dynamics, ion channel assembly and gating, and targeting proteins for ubiquitination [START_REF] Stogios | Sequence and structural analysis of BTB domain proteins[END_REF]. Contigs linked to the Ran GTPase binding function were downregulated. Ran (GTP-binding nuclear protein Ran) is a small GTPase involved in important cellular activities such as nucleocytoplasmic transport, mitotic spindle assembly, and nuclear envelope formation; it is known to be differentially expressed in response to stress [START_REF] Bo | Analysis of Ran related to pesticide resistance in Drosophila Kc cells[END_REF][START_REF] Nachury | The direction of transport through the nuclear pore can be inverted[END_REF]. Conversely, Serpin B6-like was upregulated in response to epoxiconazole exposure. Serpin B6 inhibits Cathepsin G, thereby inhibiting its functions to clear pathogens, regulate inflammation by ACCEPTED MANUSCRIPT / CLEAN COPY modifying the chemokines, cytokines, cell surface receptors, and C components, control blood pressure, and induce thrombogenesis [START_REF] Mangan | The role of serpins in vertebrate immunity[END_REF].

In the HC population, upregulation of the cytochrome P450 family protein belonging to the CYP2 family was evidenced again, as already observed in the interpopulation differences (upregulation of CYP2J2). This result confirmed the activation of a specific pathway involved in the phase I xenobiotic metabolism only in the HC population, requiring in-depth analysis of this mechanism.

As observed in the LC population, different proteins were modulated in response to epoxiconazole without an obvious explanation of their functions. We can note downregulation of tyrosinase-like protein 1; in humans, Tyrp1 is a melanocyte-specific gene product involved in melanin synthesis. This response could be a direct effect of epoxiconazole. In Mytilus galloprovincialis, evidence of tissue damage was given by a massive deposit of melanin (melanosis) after exposure to an insecticide, highlighting the inflammation processes [START_REF] Stara | Acute effects of neonicotinoid insecticides on Mytilus galloprovincialis: A case study with the active compound ACCEPTED MANUSCRIPT / CLEAN COPY thiacloprid and the commercial formulation calypso 480 SC[END_REF]. Downregulation of the c-Myc-binding protein-like was also observed; this is a proto-oncogene that functions as a transcription factor, thought to regulate expression of more than 15% of cellular genes [START_REF] Dang | The c-Myc target gene network[END_REF]. Its expression promotes cell proliferation and genomic instability by accelerating cells through G1 and S phases of the cell cycle, abrogating cell cycle checkpoints, and increasing cell metabolism [START_REF] Mladinic | The use of FISHcomet to detect c-Myc and TP 53 damage in extended-term lymphocyte cultures treated with terbuthylazine and carbofuran[END_REF].

Different proteins linked to immune responses (lymphocyte cytosolic protein 2-like isoform X2, C1q domain-containing protein) were also modulated in response to epoxiconazole exposure. Pesticides are known to be immunotoxic and interfere with specific immunological functions of each type of immune cell (see [START_REF] Lee | Adverse effects of pesticides on the functions of immune system[END_REF] for a review).

A chemoattractive glycoprotein (ES20) was found to be downregulated in response to epoxiconazole exposure. Alarm pheromones have been detected in earthworms [START_REF] Ressler | Alarm pheromone in the earthworm Lumbricus terrestris[END_REF], which can deter other members of the species but can act as a chemoattractant to other animals such as snakes [START_REF] Jiang | Purification and characterization of a chemoattractant from electric shock-induced earthworm secretion, its receptor binding, and signal transduction through the vomeronasal system of garter snakes[END_REF]. In both populations, the chitin metabolic pathway was enriched due to the upregulation of chitinase. This response could be explained by the fact that epoxiconazole is known to induce stimulation of chitinase in wheat as a defense against fungi [START_REF] Siefert | Epoxiconazole-induced stimulation of the antifungal hydrolases chitinase and β-1, 3glucanase in wheat[END_REF].

Implications of long-term pesticide use for soil health

The knowledge required to assess the implications of pesticide residue to soil biota and their repercussions on soil security is far from being complete, letting a wide gap in the way regulations for pesticide approval and policies on environmental health are formulated. Within regulatory risk assessment of pesticides, there is particularly a need for more research on long-term impacts as highlighted by "EFSA (European Food Security Agency) Scientific Opinion addressing the state of the science on risk assessment of plant protection products for in-soil organisms" (EFSA, 2017). Our findings suggest on-going evolutionary responses of soil fauna to long-term anthropogenic pressure.

Soil fauna adaptation to soil chemical pressure is a pivotal question. As soil engineer, earthworm tolerance to pesticides may hold important support for the soil sustainability and the agroecosystem ACCEPTED MANUSCRIPT / CLEAN COPY resilience. However, the costs to the population of such an acclimation/adaptation, especially in terms of reproduction, are still unknown. Our work advocates for further investigations that explicitly integrate the evolutionary effect of pesticides on soil biota and calls also for considering long-term effect of pesticides in deriving environmentally safe concentrations or registering new pesticide active ingredients.

Conclusion

Overall, we demonstrated for the soil engineer A. caliginosa, that two same lineage populations deriving from similar pedo-climatic environments, constitutively exhibit distinct gene expression pathways after they long-lived in differently managed agricultural soils. The main difference observed was regarding metabolism, with pathways linked to proteolytic activities and to the mitochondrial respiratory chain.

Other specific responses were observed and should be further investigated, such as the cytochrome P450 family CYP2 or responses linked to ion binding. It is important to keep in mind that these results need further in-depth investigations. Indeed, with the low annotation rates of the earthworm transcriptome (about 20%, which is standard for de novo assemblies), what we have discovered in the current study might be just the tip of the iceberg. The recent publication of the genome of A. caliginosa [START_REF] Perry | Molecular insights into high-altitude adaption and acclimatisation of Aporrectodea caliginosa[END_REF] and a future annotated A. caliginosa genome/transcriptome can empower the approach pursued in this study and will also aid future research. Finally, studies that integrate responses from the molecular scale to the individual and population scale should be pursued in order to understand the physiological processes involved and to predict their potential repercussions at higher levels of biological organization.

Credit author

Figure 1 .

 1 Figure 1. (A) Maximum likelihood tree estimated from the COI sequences of the study Aporrectodea caliginosa specimens (H1 -H9) and reference sequences of A. caliginosa spp. from Genbank. The TN93+G+I model was used as the model for nucleotide substitution. Branch lengths are proportional to the estimated number of nucleotide substitutions. The numbers above the branches are bootstrap values. (B) Genetic diversity in A. caliginosa populations. The sample size (N) and the following measures of COI genetic diversity are given: number of haplotypes (Nh), number of polymorphic sites (S), haplotype diversity (Hd) and nucleotide diversity (π).

Figure 2 .

 2 Figure 2. (A) Top hit species distribution on the basis of best sequence alignments and lowest E values and (B) Distribution of Gene Ontology (GO) assignments of assembled A. caliginosa contigs. GO categories are shown on the x-axis grouped into three main categories: biological processes, cellular components and molecular functions. The y-axis indicates the percentage and total number of genes in each category.

Figure 3 .

 3 Figure 3. Multidimensional scaling plots (MDS) generated using the limma plotMDS function where A) distances correspond to leading log-fold change between samples and B) distances between samples correspond to the leading biological coefficient of variation (BCV). In both plots, dimension 1 separates LC and HC populations, indicating the paired nature of the samples in each population. (T0: beginning of the experiment; J2: two days after exposure, Basal: control individuals, Stress: exposed individuals)

Figure 4 .

 4 Figure 4. Venn diagram depicting the total number of differentially expressed contigs (A) and Venn diagram depicting the total number of annotated differentially expressed contigs (B) when comparing LC and HC populations under stress and basal conditions, as well as the number of common contigs in both conditions. Heatmap of the 59 differentially expressed contigs shared (FDR<0.05) by both LC and HC populations (C). Comparison of enriched biological processes (D) and enriched molecular functions (E) between LC and HC populations under basal and stress conditions (numbers indicate the number of significant contigs corresponding to a given function, followed by the number of annotated contigs).

Figure 5 .

 5 Figure 5. Venn diagram depicting the total number of differentially expressed contigs (A) Venn diagram depicting the total number of annotated differentially expressed contigs (B) and Heatmap of differentially expressed contigs (C) in the LC and HC populations after epoxiconazole exposure (stress condition).

Figure 6 .

 6 Figure 6. Epoxiconazole concentrations in soil (n=3) and earthworms (n=5) as well as earthworm bioaccumulation factor (BAF) on days 0, 7, and 28. Results are expressed as a mean ± standard deviation (A). Weight change of earthworms from the HC and LC sites on days 7 and 28. Values are the mean out of 30 replicates (N=30) ± standard deviation (SD) (B). Metabolic rate (μg CO2 g -1 worm (fresh weight) h -1 ) of earthworms from the HC and LC sites exposed to epoxiconazole (STRESS) or not (BASAL) on days 0, 7, and 28. A dot on the boxplot indicates a mean value (N= 10 worms per group) (C). Different letters (a or b) denote statistical differences (p<0.05).

  Statement Audrey Barranger: Conceptualization, Methodology, Formal analysis, Investigation, Data Curation, Visualization, Writingoriginal draft, Project administration, Supervision. Christophe Klopp: Methodology, Formal analysis, Data Curation, Writingreview & editing. Barbara Le Bot: Methodology, Formal analysis, Resources, Writingreview & editing. Gaëlle Saramito: Methodology, Formal analysis, Data Curation. Lise Dupont: Methodology, Formal analysis, Data Curation, Writing review & editing. Stéphanie Llopis: Investigation. Claudia Wiegand: Funding acquisition, Writing review & editing. Françoise Binet: Conceptualization, Methodology, Investigation, Resources, Supervision, Writingreview & editing, Funding acquisition, Project administration.

Acknowledgments

The present study was funded by the Brittany region (France), the François Sommer Foundation (BUZHUG Project -18XZ316-01D) and the Regional Directorate for the Environment, Planning and Housing (DREAL) through the PHYTOSOL project (N°EJ 2201157402). The authors would like to thank Philippe Leroy (INRAE) and Christian Mogis (GAEC La Mandardière) for providing history of ACCEPTED MANUSCRIPT / CLEAN COPY the cultivation and access to the sampling sites. Thanks too to Yves-Laurent Dusabyinema (Genotoul, Toulouse) for RNA sequencing. The authors are also grateful to Gaëlle Daniele and Emmanuelle Vulliet (ISA, Lyon) for multi-residual pesticide analyses. Finally, we wish to thank Coup de Puce-Expansion for their professional editing services by a native English speaker.

Data availability statement

RNA-seq raw sequences and the de novo assembled transcriptome assemblies have been deposited in NCBI under BioProject PRJNA883218. This Transcriptome Shotgun Assembly project has been deposited at DDBJ/EMBL/GenBank under the accession GKBY00000000. The version described in this paper is the first version, GKBY01000000.