Stephen Mulero
email: stephen.mulero@outlook.com

Eve Toulza
email: eve.toulza@univ-perp.fr

Anaïs Loisier
email: anais.loisier@gmail.com

Meryl Zimmerman
email: zimmermann.meryl@gmail.com

Jean-François Allienne
email: allienne@univ-perp.fr

Joséphine Foata
email: foata_j@univ-corse.fr

Yann Quilichini
email: quilichini_y@univ-corse.frjpp

Jean-Pierre Pointier
email: pointier@univ-perp.fr

Olivier Rey
email: olivier.rey@univ-perp.fr

Jérôme Boissier
email: boissier@univ-perp.fr

Yann Quilichini

Malacological survey in a bottle of water: a comparative study between manual sampling and environmental DNA metabarcoding approaches

Keywords: biomonitoring, Corsica, environmental DNA, malacology, metabarcoding, snail communities

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

INTRODUCTION

Anthropogenic activities contribute to habitat fragmentation [START_REF] Strayer | Freshwater biodiversity conservation: recent progress and future challenges[END_REF] or pollution [START_REF] Blettler | Freshwater plastic pollution: Recognizing research biases and identifying knowledge gaps[END_REF] that have adverse effects on ecosystems [START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF], including an important decrease of biodiversity (Waldron et al., 2017). Freshwater ecosystems appear to be particularly threatened by human activities leading to local species extinction [START_REF] Blettler | Freshwater plastic pollution: Recognizing research biases and identifying knowledge gaps[END_REF]. According to the IUCN Red list, 46% of these ecosystems worldwide are endangered or vulnerable [START_REF] Janssen | European Red List of Habitats Part 2. Terrestrial and freshwater habitats[END_REF]. These ecosystems host 10% of all known species despite representing only 0.8% of the Earth's surface (Balian et al., 2008, Strayer and[START_REF] Strayer | Freshwater biodiversity conservation: recent progress and future challenges[END_REF]. They also provide important ecological services for the development of human populations [START_REF] Carpenter | State of the World's Freshwater Ecosystems: Physical, Chemical, and Biological Changes[END_REF]. Among freshwater species, gastropods have received particular attention for several reasons. First, they are suffering massive extinction notably due to habitat loss/degradation, the introduction of alien species and, in some cases, overexploitation [START_REF] Bouchet | Lists of protected land and freshwater molluscs in the Bern Convention and European Habitats Directive: are they relevant to conservation[END_REF][START_REF] Johnson | Conservation status of freshwater gastropods of Canada and the United States[END_REF][START_REF] Lydeard | The global decline of nonmarine mollusks[END_REF][START_REF] Strong | Global diversity of gastropods (Gastropoda; Mollusca) in freshwater[END_REF]. Second, they play a key role in the functioning of most freshwater habitats. They are involved in the cycle of several biochemicals especially nitrogen [START_REF] Hill | Nitrogen processing by grazers in a headwater stream: riparian connections[END_REF] and are the prey of several fishes (Dillon, 2000b).

Freshwater snails are also effective shredders and processors of vegetal matter in freshwater environment hence contributing to local matter and nutrient cycling (Dillon, 2000a). Finally, several freshwater gastropods are of medical or veterinary importance as they constitute intermediate hosts for several pathogens (e.g. Trematoda) [START_REF] Lu | Snail-borne parasitic diseases: an update on global epidemiological distribution, transmission interruption and control methods[END_REF]. Under the ongoing global changes, these diseases are currently (re-)emerging worldwide at an alarming rate [START_REF] Jones | Global trends in emerging infectious diseases[END_REF][START_REF] Kincaid-Smith | Emerging Schistosomiasis in Europe: A Need to Quantify the Risks[END_REF]. Hence, monitoring these snail species is crucial to address risk maps and thus prevent the transmission of these infectious diseases to nearby human and animal populations [START_REF] Mulero | Environmental DNA for detecting Bulinus truncatus: A new environmental surveillance tool for schistosomiasis emergence risk assessment[END_REF][START_REF] Jones | Confirmation of Galba truncatula as an intermediate host snail for Calicophoron daubneyi in Great Britain, with evidence of alternative snail species hosting Fasciola hepatica[END_REF].

Unfortunately, the monitoring of freshwater snails is challenging because: (i) there is no standardized monitoring protocol [START_REF] Tallarico | Freshwater Gastropods as a Tool for Ecotoxicology Assessments in Latin America*[END_REF]; (ii) their populations are spatially and temporally highly dynamic [START_REF] Lamy | Testing metapopulation dynamics using genetic, demographic and ecological data[END_REF]; and (iii) classical methods used to collect and to identify snail species are time-consuming, laborious and need malacological expertise (Vinarski and Kramarenko, 2015). In this context, the use of more sensitive molecular detection approaches such as environmental DNA (eDNA) are promising [START_REF] Bohmann | Environmental DNA for wildlife biology and biodiversity monitoring[END_REF].

According to Taberlet et al., environmental DNA is a "complex mixture of genomic DNA from many different organisms found in an environmental sample" [START_REF] Taberlet | Environmental DNA: For biodiversity research and monitoring[END_REF], Taberlet et al., 2012). This complex and heterogeneous community of DNA fragments originating from all surrounding organisms is either free in the environment or associated with organic or inorganic particles in water, soil or atmosphere [START_REF] Pedersen | Ancient and modern environmental DNA[END_REF], Taberlet et al., 2012). The use of eDNA approaches for monitoring targeted snail species already exist to detect invasive species such as Potamopyrgus antipodarum and Crepidula fornicata [START_REF] Goldberg | Environmental DNA as a new method for early detection of New Zealand mudsnails (Potamopyrgus antipodarum)[END_REF][START_REF] Miralles | Development and validation of eDNA markers for the detection of Crepidula fornicata in environmental samples[END_REF][START_REF] Clusa | eDNA for detection of five highly invasive molluscs. A case study in urban rivers from the Iberian Peninsula[END_REF]; or gastropod species that transmit parasitic diseases such as Galba truncatula, Austropeplea tomentosa, Bulinus truncatus and Oncomelania hupensis [START_REF] Jones | Detection of Galba truncatula, Fasciola hepatica and Calicophoron daubneyi environmental DNA within water sources on pasture land, a future tool for fluke control?[END_REF][START_REF] Mulero | Environmental DNA for detecting Bulinus truncatus: A new environmental surveillance tool for schistosomiasis emergence risk assessment[END_REF][START_REF] Fornillos | Detection of Schistosoma japonicum and Oncomelania hupensis quadrasi environmental DNA and its potential utility to schistosomiasis japonica surveillance in the Philippines[END_REF]. In such a biodiversity assessment context, only one eDNA tool have very recently been developed to characterize communities of freshwater bivalves [START_REF] Prié | Environmental DNA metabarcoding for freshwater bivalves biodiversity assessment: methods and results for the Western Palearctic (European sub-region)[END_REF]. However, still is lacking a general eDNA-based tool to characterize the communities of freshwater gastropods.

To fill this gap, we here present the development of a ready-to-use eDNA-based diagnostic protocol for monitoring freshwater snail communities and its application in Corsica island. This study is constructed around four objectives (i) Developing a step-by-step environmental metabarcoding protocol. (ii) Comparing this protocol with a classical visual malacological survey approach in 19 sites from 10 rivers across Corsica Island (France). (iii) Providing an updated vision of the local malacofauna at the Island scale with focus on both invasive species and species with medical importance. (iv) Discussing the pros and cons of the emerging eDNA approaches compared to classical prospecting methods.

MATERIALS AND METHODS

Malacological survey and water sampling

The field work was conducted during August and September 2018 in the Corsica Island. This summer end period fulfil optimal conditions for the presence of freshwater snails populations that reach their highest density [START_REF] Kincaid-Smith | Emerging Schistosomiasis in Europe: A Need to Quantify the Risks[END_REF]. Overall, the sampling campaign consisting of both a water sampling for eDNA and a malacological survey were performed at 23 sites distributed over 13 different rivers (Figure 1). These sites are known for human recreational activities and are annually monitored by the French agency "Agence Régionale de Santé" de Corse (ARS, 2019).

Water sampling

Water samples for subsequent eDNA analyses were collected before the malacological survey and following a previously optimized protocol [START_REF] Mulero | Environmental DNA for detecting Bulinus truncatus: A new environmental surveillance tool for schistosomiasis emergence risk assessment[END_REF]. At each sampling site, three liters of water were filtered in a single position from the river streambed just under the water surface. The water was collected downstream the prospected site to capture flowing eDNA from upstream transects, covering potentially the overall gastropod communities at a given site. Water filtrations were achieved using a sterile filtration unit with a 0.45 µm PES membrane of 90 mm diameter (VWR PES filter unit 1 L, model 514-0301) connected to a manual vacuum pump (Mityvac model MV8500). Afterwards, each filtration membrane was immediately stored in a 50 mL tube with a sterile Longmire buffer [START_REF] Longmire | Use of "lysis buffer" in DNA isolation and its implication for museum collections[END_REF]. Membranes were then stored in a dark cold room at 4°C until subsequent molecular analyses. During all the sampling process, precautions were taken to avoid any contamination [START_REF] Taberlet | Environmental DNA[END_REF]: operators used disposable sterile gloves and all reusable materials (i.e. scissors and forceps) were decontaminated before and after each site by successively placing instruments in a 10% bleach bath, in a 90% ethanol bath and then in a DNA AWAY ™ bath (Thermo Scientific) finalized by a flame sterilization [START_REF] Taberlet | Environmental DNA[END_REF].

Malacological survey

At each site, the malacological survey was achieved along a total of 300 meters by the same operators. Three transects of 50 meters (downstream, at the referenced site and upstream) were prospected along the two shores of the river (3 x 2 transects). Within each transect, all freshwater snails were sampled manually or by scooping the grass on the water bench using a colander. Snails were then grouped according to their morphology and taxonomically identified at least at the genus, at best at the species level based on their morphological traits. The number of individuals of each species was semi-quantitatively estimated. If more than 100 snails were collected, we have attributed a ">100" score to the given species. At each site, GPS coordinates and several variables were collected including: the sampling day, the half day of sampling (i.e. A.M. or P.M.), water temperature & pH, water flow and altitude.

Molecular approaches

All pre-PCR molecular steps were performed under a sterile PCR hood decontaminated before and after each use as follow: the working surface was successively washed using 90% ethanol, 10% bleach and a DNA AWAY ™ solution (Thermo Scientific ™ ) and then exposed to UV light for 30 minutes. The reusable materials were decontaminated following the same protocol as for the field sampling and treated under UV exposure for 30 minutes.

DNA extractions

Total eDNA from each membrane was extracted using the DNeasy® Blood & Tissue kit (QIAGEN) following an adapted protocol [START_REF] Mulero | Environmental DNA for detecting Bulinus truncatus: A new environmental surveillance tool for schistosomiasis emergence risk assessment[END_REF]. Briefly, membranes were subdivided into four equal parts to allow the use of 2 ml microcentrifuge tubes, these four parts were individually lysed at 65°C for one hour in a solution containing 567 µl of ATL buffers and 63 µl of proteinase K solution. Then a solution containing 630 µl of AL buffers and 630 µl of 100% ethanol was added in each tube. The remaining steps of the protocol were performed following the manufacturer instructions, except for that DNA was ultimately eluted in 100 µl (2 * 50 µl) of AE buffer pre-heated to 65°C. The resulting DNA extracts were then preserved at -20°C until subsequent steps.

In parallel, DNA was extracted from one individual of each of the 11 snail species collected during the malacological survey to produce a synthetic "mock community". To this aim, we used the E.Z.N.A.® Tissue DNA kit (OMEGA bio-tek, Inc) following the "tissue protocol". These snail DNA extractions were used as positive PCR and sequencing controls in the following steps.

To refine the molecular species assignations obtained from subsequent NGS analyses, we sequenced a 16S rRNA region for six snail species collected in Corsica. The 16Sbr-H primers were used following the PCR conditions described in [START_REF] Saito | Phylogeography of freshwater planorbid snails reveals diversification patterns in Eurasian continental islands[END_REF]. The resulting amplicons contain the barcode region used for the metabarcoding. The PCR products were then sequenced on an ABI 3730xl sequencer at the GenoScreen platform (Lille, France). Sequences were obtained for Ancylus fluviatilis (MT361136), Gyraulus sp. (MT361134), Physa acuta (MT361133), Pisidium casertanum (MT361132), Theodoxus fluviatilis (MT361131) and Galba truncatula (MT361135).

PCR primers

The primer pair used in this study (Gast01) was previously developed and tested in silico [START_REF] Taberlet | Environmental DNA: For biodiversity research and monitoring[END_REF] but never tested and/or applied for empirical purposes. These primers target a 60-70 bp fragment of 16S mitochondrial rRNA using the Gast01F (5' CCGGTCTGAACTCAGATCA 3') and the Gast01R (5' TTTGTGACCTCGATGTTGGA 3') primer pair. Based on their in silico analyses, these primers theoretically amplify at least 1,280 species distributed among 456 genera, mainly related to Gastropoda. In our study, we first validated that these primers were able to produce a size-expected amplicon using qPCR for the 11 snail species collected during the malacological survey and eight additional snail species from previous field works (Table 1). These 19 DNA extracts were amplified by qPCR on a LightCycler®480 qPCR device following the same protocol as for the first PCR step of the library preparation, see 2.2.3 below.

Library preparation and MiSeq amplicon sequencing

Five negative controls and six positive controls were processed along with the environmental samples. Negative extraction controls consisted in ultrapure water samples processed following the same DNA extraction protocol as for the filtration membranes in each extraction runs (N = 2).

PCR negative controls consisted in PCR reactions performed using water as template in each PCR run (N = 3). Positive controls consisted in two categories of mock communities. The first category consisted in an equimolar pool of 12 DNA extracts obtained from individuals of each snail species identified in Corsica and during previous field missions (Table 1). These mock communities are useful to detect PCR competition among different snail species. The second category of mock community consisted in pooling an equimolar quantity of PCR products individually obtained from 16 freshwater snail species (i.e. including the 11 identified snail species collected in Corsica and five additional non-endemic species bred in our laboratory). These controls are useful to detect possible biases during the sequencing process.

Individual libraries were generated for each membrane section of the 23 filtration membranes and for each control, hence, representing a total of 103 libraries (23 x 4 = 92 eDNA samples + 11 controls). NGS libraries were prepared following the Illumina two-step PCR protocol, using mitochondrial 16S primers with Illumina adapters for the first locus-specific PCR.

Both PCR were performed using the NEBNext ® Ultra ™ II Q5 ® Master Mix M0544L (NEW ENGLAND BioLabs, USA). The first PCR reactions were performed in a final volume of 35 µl containing 3.5 µl of extracted DNA (concentrations ranging from 0.05 ng/µl to 8.52 ng/µl), 17.5 µl of 2X Q5 Master Mix and 14 µl of indexed primer mix (final concentration 0.4 µM) and ran on a Techne TC-PLUS thermal cycler PCR device (Techne, UK) using the following program: 98°C for 30 seconds followed by 40 cycles at 98°C for 10 seconds, 55°C for 20 seconds, 65°C for 8 seconds and a final elongation step of 2 minutes at 65°C. At this point, all PCR products were checked on a 2.5% agarose gel stained with ethidium bromide after electrophoresis.

The PCR products were individually indexed in the second PCR step consisting in eight cycles using the Nextera ™ XT Index (Illumina, San Diego, USA) following the manufacturer's instructions. Finally, the libraries were normalized using SequalPrep ™ plates (Thermo Fischer Scientific, USA) before pooling. The pooled libraries were then purified following the JetSeq™ Clean protocol (Bioline, UK), checked on a Bioanalyzer High Sensitivity DNA kit (Agilent, USA) and quantified using a Qubit fluorometric quantification (Thermo Fisher Scientific, USA). Pairedend sequencing (2 x 250 cycles) was performed with a MiSeq Reagent Kit v2 on an Illumina MiSeq ™ instrument at the Bio-Environnement platform (University of Perpignan Via Domitia, France).

Data analyses

Bioinformatic pipeline

The resulting amplicon sequence dataset was processed using the Find Rapidly OTUs with Galaxy Solution (FROGS) pipeline implemented in Galaxy [START_REF] Escudie | FROGS: Find, Rapidly, OTUs with Galaxy Solution[END_REF] at the Genotoul platform (Toulouse, France). (i) The amplicon dataset was first pre-processed, according to Gast01 primer specificities (i.e. amplicon size of 60 -70 bp), we filtered out the sequences so as to keep amplicon sizes from 95 to 120 nucleotides (i.e. corresponding to the size of the targeted amplicon with Illumina adapters). (ii) The sequences kept were next clustered in operational taxonomic units (OTUs) using the swarm algorithm and using denoising and an aggregation distance of three nucleotides. (iii) The dataset was filtered out for chimeras using VSEARCH [START_REF] Rognes | VSEARCH: a versatile open source tool for metagenomics[END_REF].

Singletons and underrepresented clusters (i.e. <10 sequences) were removed as they correspond most likely to PCR error than rare sequences (Valentini et al., 2016[START_REF] Prié | Environmental DNA metabarcoding for freshwater bivalves biodiversity assessment: methods and results for the Western Palearctic (European sub-region)[END_REF]. (iv) Each OTU was assigned to a species or taxon through a two-step MEGABLAST affiliation procedure.

The first MEGABLAST analysis was restricted to Mollusca species. The second MEGABLAST analysis was performed without restricting parameters to check the robustness of the previous affiliations. The 10 best hits were kept for subsequent analysis. (v) The resulting OTUs were filtered following two criteria: first, only OTUs presenting a minimal blast coverage of 84% of amplicon length were kept; second, only OTUs presenting a pairwise identity above 89% with the affiliated sequence were kept, these thresholds are based on previously published protocols [START_REF] Macher | Comparison of environmental DNA and bulk-sample metabarcoding using highly degenerate cytochrome c oxidase I primers[END_REF] and on intraspecific sequences variations among snail strains. The remaining OTUs were considered as "unassigned".

Lastly, to confirm that a given OTU was detected in a specific sample, we have adapted our validation criteria depending on OTUs abundances. An abundant OTU (i.e. ≥1000 copies considering all the samples) was considered present in a specific membrane when reaching at least 10 copies (overall the four membrane sections). Regarding poorly abundant OTU (i.e. < 1000 copies considering all the samples), we considered a sample positive if four or more copies of such OTU were distributed among two quarters of a same membrane.

RESULTS

Malacological survey

Overall, 11 distinct freshwater snail species were identified among the 23 sampled sites (Table 1-S1). These species were already registered in Corsica [START_REF] Inpn | National Inventory of Natural Heritage [Online[END_REF], IUCN, 2020) where 18 species of freshwater mollusk have been described so far, five of which inhabiting lentic ecosystems (Table S1) and [START_REF] Mouthon | Les mollusques dulcicoles -Données biologiques et écologiques -Clés de détermination des principaux genres de bivalves et de gastéropodes de France[END_REF][START_REF] Inpn | National Inventory of Natural Heritage [Online[END_REF]. The species richness ranged from one (site 4) to six (site 23). The most common species was P. antipodarum (22/23 sites) and the rarest were T. fluviatilis (1/23 sites) and Bithynia tentaculata (3 individuals at site 13). Among the species identified, two are invasive (i.e. P. antipodarum and P. acuta) and two are intermediate hosts for Human and livestock parasites (i.e. B. truncatus and G. truncatula). Our semi-quantitative approach showed that P. antipodarum and A. fluviatilis were the most abundant species with 48.03±46.38 and 26.68±40.53 individuals, respectively (Table 1).

Molecular results

Over the 19 snail DNA extracts tested in qPCR with metabarcoding primers (Table 1), all but P. casertanum were positive. Regarding the equimolar pools of DNA, whatever pre-(Figure S2a) or post-PCR (Figure S2b) the sequencing of equimolar pools gave a number of reads for each species fitting to the expected ratio. This was higher evenness for post-PCR pools. Despite all our precautions, four among the 23 eDNA water samples (from sites 10, 13, 18 and 19 -see Figure 1) were possibly contaminated and were thus discarded. In fact, amplicons corresponding of A. fluviatilis, G. truncatula and B. truncatus were found in the negative PCR controls associated with these samples. Hence, the subsequent metabarcoding data analysis was performed on the 19 remaining sites distributed over 10 rivers of Corsica.

Excluding the four contaminated sites, the whole sequencing (controls + field samples) generated more than 14 M of clusters with an average of 97,000 sequences per library. After analyses, 8.7 M of sequences were affiliated (Figure 2). These sequences were grouped among 488 OTUs with 43 OTUs corresponding to Mollusca species. After removing the sequences from control samples, on the 7,290,067 remaining sequences, 6% (445,638) were affiliated to Mollusca, including 424,580 sequences affiliated to Gastropoda; 10% (715,308) were unaffiliated and 84% (6,129,121) were affiliated to non-mollusk species (Figure 2). Whatever the taxonomic group, the OTU affiliations were generally limited to the genus and rarely reached the species level (e.g. only four out of 11 freshwater snail species found by metabarcoding were affiliated to a single species name). However, the 16S sequences obtained from snail DNA extracts collected in Corsica, allowed to recover the corresponding species name for major Mollusca OTUs. All the genus or species detected using eDNA corresponded to species already identified in Corsica and no new genus was detected compared to our malacological survey. Hence, considering all prospected sites, we detected 61.1% (11/18) of historically known freshwater mollusks in Corsica and 84.6% (11/13) of species occurring in lotic systems (Table S1) and [START_REF] Inpn | National Inventory of Natural Heritage [Online[END_REF][START_REF] Mouthon | Les mollusques dulcicoles -Données biologiques et écologiques -Clés de détermination des principaux genres de bivalves et de gastéropodes de France[END_REF]. The detected snail species belong to three subclasses of the Gastropoda (i.e. Neritimorpha, Caenogastropoda, and Heterobranchia) (Figure 3). The non-detected subclasses (e.g.

Vetigastropoda and Patellogastropoda) are related to marine snail species. Beyond gastropods species, the present protocol allowed the detection of the bivalve P. casertanum.

Comparing the malacological survey to the eDNA monitoring

The malacological survey and the eDNA monitoring provided similar results (Figures 3 &4.). At the Island scale, all species identified during the malacological survey were also detected by eDNA monitoring. For some species such as the two Gyraulus species, the eDNA monitoring accuracy was limited to the genus level, one of these Gyraulus species has also not been affiliated by classical malacological survey either.

At the site level, the eDNA monitoring confirmed the detections obtained by malacological survey for 97.1% (67/69) of the detections. Only two occurrences of A. fluviatilis (site 17) and Gyraulus sp. (site 20) detected visually during the malacological survey, were not detected by eDNA monitoring (Figure 4). Conversely, the malacological survey confirmed the detections obtained by eDNA monitoring for 77% (67/87) of the detections (Figure 4). Moreover, eDNA monitoring detected on average 1.8±2.84 more species per site than the malacological survey although this difference was not significant (Wilcoxon rank test; W = 24, P = 0.11). Based on the eDNA approach, B. truncatus was detected at 11 sites while only visually detected at 4 sites (Figures 4 & 5b). At the river scale, the results obtained from direct malacological survey and from the eDNA monitoring approach gave even more congruent results (Figure 3). Among all the freshwater species detected in Corsica, two species are hosts for the transmission of trematodes of medical and veterinary importance (i.e. B. truncatus, Figure 5b and G. truncatula, Figure 5e) and two other are invasive species (i.e. P. antipodarum, Figure 5d and P. acuta, Figure 5f). The two invasive species display a wide distribution range (17 and 9 sites among the 19 monitored sites irrespective of the method) while the two species of medical/veterinary importance were mainly distributed in the Southern part of the Island. We were able to compare the current distributions of B. truncatus and P. antipodarum obtained within this study with those obtained in 1966, during previous malacological survey realized in Corsica for these two species (Doby et al., 1966a). The distribution of B. truncatus has greatly diminished over the last 60 years and is currently limited to 4 rivers of the Porto-Vecchio region (i.e. Cavu, Solenzara, Osu and Tarcu; [START_REF] Berry | Schistosomiasis haematobium, Corsica, France[END_REF] for which we identified these species with the two approaches (Figure 5a

DISCUSSION

One of the key steps for monitoring ecosystems and to quantify the effect of anthropogenic activities on biodiversity is to develop new tools enabling rapid and exhaustive characterization of the overall intrinsic species diversity. This is particularly true for freshwater ecosystems that undergo rapid changes and that are particularly threatened nowadays [START_REF] Carpenter | State of the World's Freshwater Ecosystems: Physical, Chemical, and Biological Changes[END_REF]. So far, the use of environmental DNA (eDNA) has been shown to be an efficient tool for monitoring fish [START_REF] Hanfling | Environmental DNA metabarcoding of lake fish communities reflects long-term data from established survey methods[END_REF], arthropod [START_REF] Krol | How Does eDNA Compare to Traditional Trapping? Detecting Mosquito Communities in South-African Freshwater Ponds[END_REF] and bivalve [START_REF] Prié | Environmental DNA metabarcoding for freshwater bivalves biodiversity assessment: methods and results for the Western Palearctic (European sub-region)[END_REF] communities in freshwater environments. Here, we present the first application of an eDNA monitoring method enabling the non-invasive characterization of entire gastropods communities from freshwater samples.

Benchmarking the eDNA monitoring

Our comparative approach between a classical (i.e. visual) prospection and an eDNA-based malacological survey highlights the ability of eDNA to characterize faithfully freshwater snail communities from water samples. Considering the 19 sites successfully diagnosed, 97.1% (67/69) of total occurrences reported by the malacological survey were confirmed using eDNA monitoring.

A noticeable benefit of eDNA approaches compared to malacological survey is the less timeconsuming sampling effort needed per site; only one sample of three liters realized by one operator (≈ 30 minutes from sampling to sample preservation) is needed [START_REF] Mulero | Environmental DNA for detecting Bulinus truncatus: A new environmental surveillance tool for schistosomiasis emergence risk assessment[END_REF]. In comparison, the time needed to realize six measure units in classical malacological survey is much longer (>2 hours) depending on field practicability and experimenter skill level. However, whatever the visual taxonomical identification of specimen under binocular or the molecular lab experiment for metabarcoding are both time consuming activities. In the present study, we did not try to optimize the volume of water sampled. However, in the case of larger sampling campaigns it might be interesting to determine the lowest volume of water necessary to ensure the same diagnosis efficiency. We recently shown, using a targeted approach, that lower water volumes down to one liter allows detecting B. truncatus without losses in efficacy [START_REF] Mulero | Environmental DNA for detecting Bulinus truncatus: A new environmental surveillance tool for schistosomiasis emergence risk assessment[END_REF].

Using the eDNA monitoring approach, we tended to detect more species per site than with the malacological survey. Indeed, over the 87 detection events found using the eDNA monitoring, 77% of them (67) were corroborated by visual observations during the malacological survey. This apparent higher richness found by eDNA metabarcoding approach compared to classical survey was already reported in previous studies [START_REF] Deiner | Environmental DNA metabarcoding: Transforming how we survey animal and plant communities[END_REF][START_REF] Hanfling | Environmental DNA metabarcoding of lake fish communities reflects long-term data from established survey methods[END_REF], Valentini et al., 2016[START_REF] Prié | Environmental DNA metabarcoding for freshwater bivalves biodiversity assessment: methods and results for the Western Palearctic (European sub-region)[END_REF]. Three non-exclusive hypotheses could explain the higher efficacy of eDNA compare to classical methods. First, DNA transportation through the river flow allows collecting eDNA from upstream the targeted site. In fact, in lotic systems, eDNA is known to disperse in a given way with waterflow [START_REF] Hanfling | Environmental DNA metabarcoding of lake fish communities reflects long-term data from established survey methods[END_REF][START_REF] Deiner | Environmental DNA reveals that rivers are conveyer belts of biodiversity information[END_REF]. Second, eDNAbased methods are expected to outcompete visual surveys for detecting cryptic or rare species (e.g. [START_REF] Dejean | Improved detection of an alien invasive species through environmental DNA barcoding: the example of the American bullfrog Lithobates catesbeianus[END_REF]. Cryptic freshwater mollusk species are particularly difficult to detect and previous studies estimated that about 20% of the species present at a single site are undetected after a single prospecting session during malacological surveys [START_REF] Dubart | Modeling competition, niche, and coexistence between an invasive and a native species in a two-species metapopulation[END_REF]. Moreover, detecting biases are also common in malacology because populations of snails are highly dynamic in space and time [START_REF] Lamy | Testing metapopulation dynamics using genetic, demographic and ecological data[END_REF]. Interestingly, our eDNA approach has been particularly useful to detect P. casertanum compared to classical malacological survey, this species being known to be difficult to detect visually because of his preferred niche (i.e. under boulders) and its small size. Third, false positives constitute one major limitation of eDNA-based protocols [START_REF] Cristescu | Uses and Misuses of Environmental DNA in Biodiversity Science and Conservation[END_REF]. These biases are generally related to either primer amplification biases or inappropriate field and/or laboratory protocols [START_REF] Cristescu | Uses and Misuses of Environmental DNA in Biodiversity Science and Conservation[END_REF]. Regarding primer amplification biases, further applications of the current tool are needed to evaluate these errors considering that it is impossible to test the primers with the DNA of all known organisms in vitro. In our study, such biases are likely to be limited because our controlled mock communities were detected with no apparent amplification & sequencing biases between species. Concerning DNA contamination, the four sites that could have been contaminated were discarded. Moreover, the detections found using eDNA monitoring yet not confirmed by malacological survey, followed an ecologically realistic distribution. For G. truncatula, the two sites found positive only by eDNA monitoring were distributed on the Gravona River, downstream from a site where G. truncatula was detected by the two approaches (i.e. site 5). The two other sites were distributed on the Rizzanese river in which Lymnaeid snails hosts of F. hepatica were previously identified [START_REF] Gretillat | Epidemiology of Certain Trematode Diseases of Domestic Animals in Corsica (Bovine Bilharziosis and Bovine and Ovine Distomiasis). Observations Conducted during a Mission Accomplished during the Autumn of 1962[END_REF]. Noteworthy, this species prefers temporary water body, but could be found in river shore and associated trickle [START_REF] Dreyfuss | Les habitats de Lymnaea truncatula Müller (Mollusque) le long de deux rivières[END_REF], which might explain the observed relatively low densities.

Regarding B. truncatus, the detection differences are the most noticeable between the two approaches, the seven additional sites were mainly distributed in the Gravona River and nearby rivers. For each of these rivers, B. truncatus was previously identified (Doby et al., 1966a). These non-randomly distributed detections in rivers were this species is historically known to be present support rather the non-detection of certain species by malacological survey, than false positive eDNA detection. However, a thorough field survey of these rivers is needed to validate our results.

Environmental DNA for characterizing communities

A recent study focused on freshwater bivalves has proposed an eDNA metabarcoding protocol for the characterization of these communities [START_REF] Prié | Environmental DNA metabarcoding for freshwater bivalves biodiversity assessment: methods and results for the Western Palearctic (European sub-region)[END_REF]. The main differences between our study and the latter one are the volume of water sampled (30 L vs 3 L) and the use of several primer combinations in the study of Prié et al. It is interesting to notice that both approaches tend to the same conclusions: eDNA detects more species than a classical malacological survey and the sensitivity of both protocols is similar. Indeed, in a single sampling, at their spatial scale [START_REF] Prié | Environmental DNA metabarcoding for freshwater bivalves biodiversity assessment: methods and results for the Western Palearctic (European sub-region)[END_REF] have detected 60% of the known bivalve community and we have found 61% of known gastropods in Corsica. Despite different protocols and different targeted groups this show replicability but also the limitation of the eDNA approach. The present eDNA metabarcoding protocol allowed to qualitatively characterize the Corsican gastropod fauna at global scale.

However, as shown in previous eDNA-based metabarcoding studies, the accuracy of such characterization highly depends on the reliability of currently available databases together with sequence size that sometimes results in a limitation for the taxonomic affiliation to the genus taxonomic level especially for closely related species [START_REF] Krol | How Does eDNA Compare to Traditional Trapping? Detecting Mosquito Communities in South-African Freshwater Ponds[END_REF][START_REF] Deiner | Environmental DNA metabarcoding: Transforming how we survey animal and plant communities[END_REF]. Here, we couldn't distinguish the two distinct species belonging to the Gyraulus genus identified visually. These species could also be difficult to differentiate morphologically, that could explain confusions with Planorbis planorbis on public NCBI databases. At the eDNA era, such limitations clearly call for the improvement and reinforcement of public databases aiming at characterizing biodiversity at the taxonomic and molecular scale [START_REF] Henry | Integrating ongoing biodiversity monitoring: potential benefits and methods[END_REF][START_REF] Balint | Accuracy, limitations and cost efficiency of eDNA-based community survey in tropical frogs[END_REF]. To this end, we submitted 6 new 16S rRNA sequences of freshwater snail species identified in Corsica.

However, to overcome these limitations in the immediate future, [START_REF] Krol | How Does eDNA Compare to Traditional Trapping? Detecting Mosquito Communities in South-African Freshwater Ponds[END_REF] recently suggested to use a combination of barcode to optimize taxonomic assignation. Recent studies also called for the need of reinforcing the collection of metadata in eDNA studies that could influence the detection of eDNA [START_REF] Nicholson | An analysis of metadata reporting in freshwater environmental DNA research calls for the development of best practice guidelines[END_REF]. In this way, we collected environmental variables at all sampling sites including some that are generally neglected (e.g. pH, waterflow, water temperature; [START_REF] Nicholson | An analysis of metadata reporting in freshwater environmental DNA research calls for the development of best practice guidelines[END_REF]. However, the sampling of a wider diversity of freshwater environment is still needed to accurately assess the effect of these variables on our ability to detect species.

Environmental DNA for the risk assessment of parasitic diseases

Beyond the technological development, our study provides an updated picture of the distribution of snails communities in Corsica. Among the 11 species detected, two snail species are of medical/veterinary importance. An updated distribution map for these two species is of prime interest for determining the risk of the associated waterborne diseases. In particular, in 2013 an outbreak of Schistosomiasis emerged in Corsica, this disease being caused by a trematode parasite using B. truncatus as host [START_REF] Boissier | Outbreak of urogenital schistosomiasis in Corsica (France): an epidemiological case study[END_REF][START_REF] Ramalli | Persistence of schistosomal transmission linked to the Cavu river in southern Corsica since 2013[END_REF]. Bulinus truncatus was widely distributed in Corsica in the middle of the 20 th century and at this time some authors warned authorities about the potential for such a local emergence of Schistosomiasis (Doby et al., 1966b, Doby et al., 1966a). Today, this species distribution is limited in four rivers in Southern Corsica (Cavu, Solenzara, Tarcu and Osu; [START_REF] Berry | Schistosomiasis haematobium, Corsica, France[END_REF]. Using our eDNA monitoring tool, we detected the presence of B. truncatus in the Gravona River, although possibly at low densities and/or upstream the sampled sites. This result call for an urging field validation of the local presence of this species to prevent potential transmission to Humans in this highly frequented river near Ajaccio, the main agglomeration of Corsica.

Fasciola hepatica is another trematode species already present in Corsica that uses G. truncatula as intermediate host [START_REF] Oviedo | The intermediate snail host of Fasciola hepatica on the mediterranean island of Corsica[END_REF]. This parasite is responsible for the zoonotic Fasciolosis disease generally infecting livestock and Humans (Valero et al., 2002). The present study provides an up-to-date geographic distribution of G. truncatula in Corsica, hence, delineating pastures presenting a potential risk of transmission to livestock and Humans.

Interestingly, the eDNA monitoring tool also detected the presence of G. truncatula in the Rizzanese River, a region in which lymnaeids snails hosts of F. hepatica were previously identified [START_REF] Gretillat | Epidemiology of Certain Trematode Diseases of Domestic Animals in Corsica (Bovine Bilharziosis and Bovine and Ovine Distomiasis). Observations Conducted during a Mission Accomplished during the Autumn of 1962[END_REF] although not locally detected by visual inspection in this study. In a sanitary context, these results also call for thorough monitoring of this species in this particular river.

Environmental DNA for studying bioinvasions

Invasive snail species needs attention as they can drastically disturb the receiving ecosystems and impact the ecology of indigenous species by altering biotic composition of invaded ecosystems [START_REF] Cowie | Invertebrate invasions on Pacific Islands and the replacement of unique native faunas: a synthesis of the land and freshwater snails[END_REF][START_REF] Strayer | Alien species in fresh waters: ecological effects, interactions with other stressors, and prospects for the future[END_REF]. In this context, we have detected two invasive snail species: Physa acuta and Potamopyrgus antipodarum. Physa acuta originating from North America (Vinarski, 2017) has invaded Europe likely through the activity of aquarium keepers and maritime trade (Vinarski, 2017). Physa acuta was first reported in Corsica in 1977 [START_REF] Holoyak | Données de malacologie continentale transmises par[END_REF]. Today, we observe that this species has a wider distribution in Corsica than reported in public databases [START_REF] Holoyak | Données de malacologie continentale transmises par[END_REF], GBIF, 2020).

Potamopyrgus antipodarum originates from New Zealand [START_REF] Goldberg | Environmental DNA as a new method for early detection of New Zealand mudsnails (Potamopyrgus antipodarum)[END_REF] and has been introduced in Europe in 1859 probably via human activities such as fish farming (Alexandre da [START_REF] Da Silva | Modelling species distributions to predict areas at risk of invasion by the exotic aquatic New Zealand mudsnail Potamopyrgus antipodarum (Gray 1843)[END_REF]. This parthenogenetic snail is widely distributed in Corsica and needs attention considering the important density that these populations can reach (e.g. >400,000/m²; [START_REF] Hall | Exotic snails dominate nitrogen and carbon cycling in a highly productive stream[END_REF]. This snail species is known to negatively impact nutrient cycling or can competitively exclude native species [START_REF] Goldberg | Environmental DNA as a new method for early detection of New Zealand mudsnails (Potamopyrgus antipodarum)[END_REF]. Potamopyrgus antipodarum was first identified in Corsica in 1961 and have rapidly invaded rivers of the south-west side of the island (Doby et al., 1966a). Today, we observe that this species is widely distributed and abundant all around the island (Figure 5d). Interestingly, the researchers that carried out the previous malacological survey hypothesized that this species could competitively exclude other freshwater snails' such as snail intermediate hosts for parasites (Doby et al., 1966a). In the same study, the authors observed that B. truncatus hosting either Schistosoma bovis (an animal parasite historically present in Corsica) or S. haematobium (a human parasite currently present in Corsica) was present all around the island. The considerable restriction of B. truncatus distribution in parallel with the increase distribution of P. antipodarum seems to agree with Doby's et al. predictions. However, environmental metabarcoding approaches are currently limited for quantifying these interactions in natura, and more generally when they are applied for conservation, and parasitological risk assessment purposes because these study areas require not only presence/absence information but also abundance data (Yates et al., 2020). The relationship between eDNA and species abundance is subjected to unexplained variations (Yates et al., 2020) mainly because eDNA is heterogeneously released in the environment [START_REF] Jo | Effect of water temperature and fish biomass on environmental DNA shedding, degradation, and size distribution[END_REF] and DNA persistence depend on several environmental biotic and/or abiotic factors [START_REF] Barnes | Environmental conditions influence eDNA persistence in aquatic systems[END_REF][START_REF] Taberlet | Environmental DNA: For biodiversity research and monitoring[END_REF].

Conclusion

In view of these results, the eDNA monitoring using metabarcoding approaches appears particularly useful for biomonitoring studies. Here we presented the first application of eDNA metabarcoding protocol for the characterization of entire gastropods communities from water samples. This eDNA-based tool provides very similar results when compared to the classical malacological survey and is useful for the detection of habitat hosting whether invasive snail species or species of medical and veterinary importance, while requiring fewer sampling efforts than the visual prospecting. Beyond these qualitative data, our results have given more support to previous hypotheses about a potential competitive interaction between B. truncatus and P.

antipodarum. Indeed, this interaction needs to be experimentally validated but could play an important role in the endemization process of the urogenital bilharziasis in Corsica. Hence, the environmental metabarcoding is particularly suited for a first analysis of environments for which data on malacofauna are scarce. In a regular use context, the temporality of records could enable an accurate characterization of local snail communities in space and time leading to the early detection of phenomena such as biological invasions, biodiversity losses, and hence, paving the way for complementary studies in other fields.

Gyraulus laevis

Lentic, lotic 20, 23 0.06 (± 0.34) from 0 to 3 Gyraulus sp. Lentic,lotic 6,7,12,14,15,17,20,23 2.99 (± 8.66) from 0 to 50

Radix balthica

Lentic, Lotic 1, 2 0.11 (± 1.

  & b). Contrarily, P. antipodarum has invaded most of the sites we have prospected and is now present on the east side of Corsica (Figure 5c & d).

Figure 1 .

 1 Figure 1. Map referencing the 23 sites prospected across Corsica Island for malacological survey and for water eDNA sampling. All these sites are lotic and distributed among 13 different rivers (darker lines). These sites are annually monitored by the French agency ARS, coordinates and environmental variables available upon request. The sites in black were sampled but removed from the analysis due to potential DNA contaminations (see text for details).

Figure 2 .

 2 Figure 2. Major OTUs ranking based on sequence abundances. These sequences were obtained from eDNA samples collected in the field and do not contain those from experimental controls. Each taxon is based on merged OTUs (i.e. all OTUs affiliated to a same species are merged in one category).

Figure 3 .

 3 Figure 3. Simplified phylogeny of molluscs from Cunha & Giribet, (2019) highlighting all species reported in Corsica Island during the present study and species used to test primers (mock communities). Horizontal barplots represent species occurrences recorded at the river scale (N = 10) based on the malacological survey (in orange) and on the eDNA metabarcoding results (in blue). Numbers into brackets show the number of rivers positives for both methods on the total of positive rivers for the method showing the highest number of positive for a given species.

Figure 4 .

 4 Figure 4. Schematic representation of all detected snail species (or genera) based on the malacological survey (left side) and the eDNA monitoring tool (right side) at the 19 analysed sites (sites 10, 13, 18 and 19 were discarded because of possible DNA contaminations). Columns represent species and are arranged symmetrically between the two detection methods. The number in each cell is the semi-quantitative abundance of each species at each site.

  recovered from the work of Mouthon, 1982 and confirmed with the INPN 2020 database. Lentic = still water, Lotic = flowing water. The displayed median headcount for each species reported by the malacological survey consider the six measure units conducted on each sites as separate values.

Figure S2 .

 S2 Figure S2. Number of sequences generated per species in equimolar controls, a. results for pooled DNA (N = 12) and b. results for pooled PCR products (N = 16). The first three columns are associated to the three technical replicate realized for each pool category, the fourth

Table 1 .

 1 Snail species used for assessing the detection range of the Gast01 PCR primers in vitro.

	Snail species	Origin	Date	Provided by	Average headcount per site
	Bulinus forskalii	South-Africa	2017	Huyse T., KU Leuven, Belgium	_
	Bulinus globosus	Senegal	2012	Huyse T., KU Leuven, Belgium	_
	Bulinus africanus	South-Africa	2017	Huyse T., KU Leuven, Belgium	_
	Radix natalensis a,b	South-Africa	2017	Huyse T., KU Leuven, Belgium	_
	Galba schirazensis a,b	Cuba	2013	Hurtrez-Boussès S., IRD Montpellier, France	_
	Galba cubensis a,b	Cuba	2013	Hurtrez-Boussès S., IRD Montpellier, France	_
	Pseudosuccinea columella a,b	France	2011	Hurtrez-Boussès S., IRD Montpellier, France	_
	Biomphalaria glabrata a,b Brazil	2018		_
	Bulinus truncatus a,b	Corsica, France	2018		2.28 (± 13.95) from 0 to >100
	Bithynia tentaculata b	Corsica, France	2018		0.02 (± 0.28) from 0 to 3
	Radix balthica b	Corsica, France	2018		0.11 (± 1.21) from 0 to 13
	Ancylus fluviatilis a,b	Corsica, France	2018		26.68 (± 40.53) from 0 to >100
	Galba truncatula b	Corsica, France	2018		0.06 (± 0.38) from 0 to 3
	Gyraulus laevis a,b	Corsica, France	2018		0.06 (± 0.34) from 0 to 3
	Gyraulus sp. a,b	Corsica, France	2018		2.99 (± 8.66) from 0 to 50
	Physa acuta a,b	Corsica, France	2018		17.38 (± 34.29) from 0 to >100
	Potamopyrgus antipodarum a,b	Corsica, France	2018		48.03 (± 46.38) from 0 to >100
	Pisidium casertanum a,b	Corsica, France	2018		1.45 (± 7.10) from 0 to 50
	Theodoxus fluviatilis b	Corsica, France	2018		0.87 (± 5.26) from 0 to 50
	a DNA extracts used for the preparation of equimolar pools of DNA (N = 12). b DNA extracts used
	for the preparation of equimolar pools of PCR products (N = 16). The displayed median headcount
	for each species reported by the malacological survey consider the six measure units conducted on
	each sites as separate values.

Table S1 .

 S1 Freshwater mollusks that were historically and presently detected in Corsica.

					Averag
	Species name		Habitat	Identified in sites	e headcount per
					site
	Galba truncatula		Lotic, Lentic	2, 5, 11, 16	0.06 (± 0.38) from 0 to 3
					2.28 (±
	Bulinus truncatus		Lotic, Lentic	1, 2, 3, 11	13.95) from 0 to
					>100
	Bithynia tentaculata		Lotic, Lentic	13	0.02 (± 0.28) from 0 to 3
					26.68
	Ancylus fluviatilis		Lotic	All except 16	(± 40.53) from 0
					to >100
	Potamopyrgus antipodarum		Lotic, Lentic	All except 4 and 12	48.03 (± 46.38) from 0 to >100
	Physa acuta		Lentic, Lotic	3, 6, 7, 8, 9, 15, 16, 21, 23	17.38 (± 34.29) from 0 to >100
	Pisidium casertanum	ponds	Lotic, ditches,	8, 11, 15, 21, 22, 23	1.45 (± 7.10) from 0 to 50

  expected composition of the sample in number of distinct DNA ordered per genus.

	10481 25502 10752 4018 9 column give the 70% 80% 90% 100%	8765 22568 9680 3409 1 1	7038 21032 8385 3000 1 2	1 1 1 1
	60%	34137	32774	28749	1
	50%				2
	40%	41045	36241	31574
	30%	5854	5830	4834	2
	20%	27390	26584	23368	1
	0% 10%	20749	19773	17611	1 1
		Rep 1	Rep 2	Rep 3	Expected
	Ancylus sp. fluviatilis	Biomphalaria sp.		Bulinus sp.
	Galba sp.		Gyraulus sp.		Physa acuta
	Pisidium sp.	Potamopyrgus antipodarum	Pseudosuccinea columella
	Theodoxus fluviatilis
	100%	7395	6294	6530	1
	50% 60% 70% 80% 90%	12208 15499 15850 12393 11515 18335	9086 11901 13083 9553 9043 14222	9591 12233 13944 10101 9523 15033	2 1 1 1 1 2
	0% 10% 20% 30% 40%	12609 22685 7495 10591 26244	9476 17109 5551 7816 19263	10078 17895 6123 8254 20256	1 1 1 1 3
		Rep 1	Rep 2	Rep 3	Expected
	Ancylus sp. fluviatilis	Biomphalaria sp.		Bithynia tentaculata
	Bulinus sp.		Galba sp.		Gyraulus sp.
	Physa acuta		Pisidium sp.		Potamopyrgus antipodarum
	Pseudosuccinea columella	Radix sp.		Theodoxus fluviatilis
					44

Acknowledgements

We thank the French agency "Agence Régionale de Santé de Corse" (ARS) for the data relative to the prospected sites. This study is set within the framework of the "Laboratoire d'Excellence (LABEX)" TULIP (ANR-10-LABX-41).

Funding

This study was funded by the Région Occitanie, the European "Fonds Européen de Développement Régional" (FEDER); the French Agency for Food, Environmental and Occupational Health & Safety (ANSES) [PNREST 2019/1/059 Molrisk]; and the French "Agence Nationale de la Recherche" (ANR) [Program HySWARM ANR-18-CE35-0001].

Availability of data and materials

The dataset supporting the conclusions of this article is included within the article.

Authors Contributions

JB., OR. and SM. conceived the study and the sampling design. AL., MZ., JF., YQ and SM. carried out the malacological survey and the field sampling. SM., JFA. and ET. processed the samples and performed molecular experiments and amplicon sequencing. SM., JB., OR. and ET., analyzed and interpreted the data. JPP. provided the malacological expertise. SM., JB. and OR. drafted the manuscript. JB., OR., SM., ET., AL., MZ., JFA., JF., YQ. and JPP. revised the manuscript. All authors read and approved the final manuscript. and f. Physa acuta, an invasive snail species from North America. Black dots are sites were the given species is absent for the two monitoring methods; Purple dots are sites were the given species was detected with the two monitoring methods; Blue dots are sites were the given species was detected using the eDNA monitoring only; Orange dots are sites were the given species was detected using the malacological survey only; Yellow dots are sites were the given species was identified in our study in 2018 and in the study of Doby et al. in 1966.

Disclosures

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


R Ehrenreich

S Allart

H Sousa

Cegla

J V Seidel

M Lendl

V Bourrier

D Ehrenreich

R Allart

S G Sousa

H M Cegla

X Bonfils

U Conod

A Grandjean

A Wyttenbach

N Astudillo-Defru

D Bayliss

Kevin Heng

B Lavie

C Lovis

C Melo

F Pepe

D Ségransan

S Udry

Hot Exoplanet Atmospheres Resolved

Keywords: Planetary Systems -Planets and satellites, atmospheres, individual, WASP-127b -Techniques, spectroscopic -Instrumentation, spectrographs -Methods, observational

Introduction

Our knowledge of the exoplanet population has increased significantly over the last few years, with thousands of detected exoplanets for a wide range of masses, sizes, and system architectures, some of which are not found in our Solar System. Despite these advances and the consensus that exoplanets are ubiquitous in our Galaxy (Cassan et al. 2012;Dressing & Charbonneau 2013;Wright et al. 2012;Howard 2013;Batalha et al. 2013), the study of atmospheric composition and structure of exoplanets is still in its early days. One of the more important features of the currently known exoplanet population is the so-called Neptune desert, a dearth of Based on observations made at ESO 3.6 m telescope (La Silla, Chile) under ESO programmes 098.C-0306 and 100.C-0750 (PI Ehrenreich).

highly irradiated planets in the size range of Neptune (Lecavelier Des Etangs 2007; Mazeh et al. 2016). This feature is not due to observational bias; these planets should be easily accessible from both space and the ground. A proposed explanation of the Neptune desert is that these planets do not survive the photoevaporation that is driven by their host stars (Owen & Lai 2018;Owen 2019;Ehrenreich et al. 2015;Bourrier et al. 2018). The radius of one of the planets in the Neptune desert, WASP-127b (discovered by the WASP-South survey, [START_REF] Kurucz | SYNTHE spectrum synthesis programs and line data[END_REF], is comparable to that of Jupiter, but it has a super-Neptune mass, making it one of the puffiest planets to date. WASP-127b is shown in the mass-insolation space in Figure 1. Additionally, its comparatively large radius for a planet of a similar mass cannot be explained by the standard coreless model (see Fortney et al. 2007).

Date

Retained Spectra a Exp. Time Notes. (a) In parentheses: spectra in and out of transit, respectively. (b) Airmass at the beginning, centre, and end of transit. (c) Observed: 37 exposures. Partial transit. (d) Observed: 45 exposures. (e) Observed: 24 exposures. Partial transit, rejected from the analysis for insufficient out-oftransit baseline. (f) Observed: 42 exposures. WASP-127b orbits its bright G5-type host star (V=10.2) with a period of 4.18 days, but it receives a relatively low XUV flux despite its close-in position in the system (Chen et al. 2018). This makes photo-evaporation an unlikely explanation for its inflation. Potential alternatives include re-inflation by the host star, enhanced atmospheric opacity, and Ohmic or tidal heating (Perna et al. 2010;Leconte et al. 2010;Batygin & Stevenson 2010;Batygin et al. 2011;Rauscher & Menou 2013;Lithwick & Wu 2014;Lammer et al. 2016). At the time of writing, none of these concepts can be established as the clear explanation for the puffiness of WASP-127b, and more precise observations of its atmospheric composition (e.g. dos Santos et al. 2020) and system parameters are needed. Fortunately, WASP-127b is a prime candidate for spectroscopic observations, with a bright host star and an enormous atmospheric scale height (∼ 2350(∼ km, Lam et al. 2017). Previous observations with the 2.5m Nordic Optical Telescope (NOT) and the Gran Telescopio Canarias (GTC), see Palle et al. (2017) and Chen et al. (2018), have shown evidence for sodium, lithium, potassium, hazes in the form of a Rayleigh slope, and tentative indications for TiO and VO. Space-based observations using the Hubble Space Telescope (HST) and the Spitzer telescope conducted by Spake et al. (2019) were able to detect H 2 O and CO 2 bands. The observations restricted the metallicity of WASP-127b to super-solar, which has recently been confirmed, together with the cloudy nature of the atmosphere by an independent study using HST data (Skaf et al. 2020). Additionally, the results re-ported by Spake et al. (2019) are in line with the ground-based sodium observations. Žák et al. (2019), from here Z19, reported a sodium signal with an absorption of more than 1% in two transit observations of HARPS (High Accuracy Radial velocity Planet Searcher) data, which were observed as part of the HEARTS survey. The two transit observations comprised a total of four transits of WASP-127b. The analysis of the transits corresponds to a sodium detection of 0.1060 ± 0.0008 R p /R on the wavelength grid from Spake et al. (2019), which differs from the Spake et al. (2019) results at the 1 σ level (∼ 0.1035 ± 0.0010 R p /R , see Figure 21 in Spake et al. (2019)). In the following, we discuss the full available HARPS dataset for WASP-127b of four observation nights and show the effect of telluric sodium and low signal-to-noise ratio (S/N) traces stemming from the stellar sodium lines when planets with low transit depths are observed. Additionally, we present updated system parameters obtained from photometric observations with EulerCam at the Swiss telescope in la Silla, Chile, and the Transiting Exoplanet Survey Satellite (TESS).

Observations and data reduction

As part of the HEARTS survey (ESO programme: 098.C-0306, 100.C-0750; PI: Ehrenreich), we observed four transits of the inflated super-Neptune WASP-127b in front of its host star. The transits took place on 2017 February 28, 2017 March 20, 2018 February 13, and 2018 March 31 and were observed with the HARPS spectrograph at the ESO 3.6 m telescope in La Silla Observatory, Chile (Mayor et al. 2003). Because of visibility constraints of the target, two of the transits were only partially recorded and only have out-of-transit spectra before or after the transit itself. The out-of-transit spectra are defined as spectra that were recorded when no part of the planet was in front of its host star. All spectra in which the planet is partially or completely in front of the star are in-transit spectra. In four nights, a total of 141 spectra were taken, 86 in-transit and 55 out-of-transit. An observation log can be found in Table 1. For each night, we discarded the exposures in which the low S/N core of the stellar lines overlap the planetary orbital track by too much. This concerns 10, 6, 5, and 7 exposures in the four observation nights, respectively (see Section 5.1.2 and Appendix A, as well as Seidel et al. (2020), for further details). In night 1 a partial transit was recorded, of which one spectrum at the end (in-transit) was only partially taken and not discarded, but added to the transit. In the same night, three in-transit spectra showed significantly lower S/N, most likely due to passing clouds. These spectra were discarded for this analysis. In night two, seven consecutive in-transit spectra were rejected for low S/N, most likely again due to a passing large cloud. Of night three, in which a partial transit was recorded, two out-of-transit exposures were rejected because they were taken to test the best exposure time and not as part of the transit. Additionally, two in-transit spectra were rejected because their S/N was significantly lower than that of the remaining spectra (unknown origin). This left an insufficient number of out-of-transit spectra to properly correct the stellar spectrum, and night three was subsequently rejected from the analysis. In night four, four in-transit spectra were rejected because their S/N was low. This work focusses on the sodium doublet, which in the case of HARPS is located in echelle order 56, spanning the wavelength region from 5850.24 Å to 5916.17 Å. Fibre A was set on the target during the transit, and Fibre B observed the sky. The recorded HARPS spectra were extracted from the e2ds raw data by the HARPS Data Reduction Pipeline (DRS v3.5) 1 . The wavelength solution is given in the Solar System barycentric rest frame. Additionally, the data were corrected for the blaze, cosmic particles hitting the detector, and telluric contamination.

Telluric lines play a major role in transmission spectroscopy. Depending on the velocities of the observed system, the telluric lines can overlap with the stellar sodium doublet, which distorts the line shape and depth. To correct for this effect, we used molecfit version 1.5.1. (Smette et al. 2015;[START_REF] Kausch | [END_REF], an ESO tool for correcting telluric features in ground-based observations. For further details on how to apply molecfit, see Allart et al. (2017) The telluric lines were corrected down to the noise level for all airmasses. However, molecfit only corrects for local atmospheric profiles based on models and the measured weather data at the local site, but not for sky emission. It therefore cannot correct for telluric sodium emission lines. This potential contamination can be estimated from the data taken with Fibre B, which was set to the sky. We discuss the effect of telluric sodium in Section 5.1.

Simultaneous photometry with EulerCam

We obtained photometric observations with EulerCam, the charged coupled device (CCD) imager installed at the 1.2 m Euler telescope at La Silla Observatory, Chile. A detailed description of the instrument, as well as the processes needed to obtain relative aperture photometry, can be found in Lendl et al. (2012). For each transit event, stable reference stars were chosen iteratively. The light curves of full or partial transits thus obtained are shown in Figure 2. An overview of the main parameters for each observations can be found in Table 2. The photometric uncertainties of the EulerCam data do not reflect additional correlated noise due to stellar, atmospheric, or instrumental effects. To compensate for this, we added the values given in Table 2 1 http://www.eso.org/sci/facilities/lasilla/ instruments/harps/tools/archive.html quadratically to the uncertainties of each data set, enforcing a reduced χ2 of unity for each light curve. The combined 2 min EulerCam RMS is 0.000261. Additionally, one TESS light curve was observed in Sector 9, covering four individual transits, with a combined 2 min RMS of 0.000209. The PDC-SAP TESS light curve (Smith et al. 2012;Stumpe et al. 2014;[START_REF] Jenkins | Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series[END_REF] was included in our analysis, and a comparison of combined Eu-lerCam light curves and the TESS light curve is shown in Figure 3.

Revised system parameters

We analysed all available photometric transits with the 'code for transiting exoplanet analysis' (CONAN) (Lendl et al. 2020), which is based on the transit analysis code presented in Lendl et al. (2017) and uses the algorithms from Kreidberg et al. (2015) to calculate transit models. The fitted parameters (MCMC jump parameters) are listed in the top section of Table 3 and the model fit is shown as red lines in Figure 3.

We used a combined HARPS spectrum to derive the stellar spectroscopic parameters (T eff , log g, microturbulence, [Fe/H]), and respective uncertainties. For this purpose, we followed the ARES+MOOG method as described in Sousa (2014) and Santos et al. (2013). The ARES code 2 (Sousa et al. 2007(Sousa et al. , 2015) ) was used to consistently measure equivalent widths (EW) of iron lines that are included in the line list presented in Sousa et al. (2008). Briefly, ARES+MOOG performs a minimisation process searching for the ionisation and excitation equilibrium to find convergence on the best set of spectroscopic parameters. To compute the iron abundances, we used a grid of Kurucz model atmospheres [START_REF] Kurucz | SYNTHE spectrum synthesis programs and line data[END_REF]) and the radiative transfer code MOOG (Sneden 1973). The stellar mass, radius, and age were obtained through the use of the Padova stellar model isochrones (see. da Silva et al. 2006;Bressan et al. 2012) -http://stev.oapd.inaf.it/cgi-bin/param_1.3) by providing our spectroscopic parameters together with the Gaia DR2 parallax and the V magnitude. The limb-darkening values were derived using the routines by Espinoza & Jordán (2015). All updated stellar and planetary parameters are given in Table 3 and are used for the spectroscopic analysis in this work.

Transmission spectroscopy of WASP-127b

The goal of transmission spectroscopy is to extract the planetary signal in the in-transit spectra. The base of this procedure is to sum all in-transit spectra to produce the master-in and divide it by the master-out, thus taking out the stellar spectral lines from the transmission spectrum R (Brown 2001). However, groundbased observations are not suited for the direct application of this method because the flux changes with time as a result of the varying airmass. Each in-transit spectrum was therefore fitted by a third-order polynomial and then divided by the master-out in the stellar rest frame to guarantee comparable flux levels in all spectra. These corrected in-transit spectra were then shifted to the planet rest frame to align the planetary sodium signal, and they were summed together for higher S/N (Wyttenbach et al. 2015). The full description of this method can be found in Seidel et al. (2019).

Two of the four transits (nights one and two) from our programme have been analysed in Z19. Z19 analysed the two transits taken in 2017, which were at the time publicly available in the archive, corresponding to 82 spectra (52 in-transit). Based on this subset of our data, Z19 claimed a sodium detection at 5 in Z19)3 . In the following we reproduce this analysis and highlight how telluric sodium emissions and low S/N remnants can produce falsepositives.

Contamination sources in the sodium wavelength range

When we look at the data of each night individually, the sodium detection is curiously strong in night two, while a tentative sodium detection is seen in night one, but not in night four. Because at the time of analysis, Z19 only had access to the first two nights and because night one is a partial transit, their sodium detection is dominated by data from night two. An in-depth analysis of night two reveals an emission feature in the master-out sodium line core when compared with the master-out from the other three nights (see Figure 4). This emission causes the stellar line to appear more shallow. When the too shallow stellar sodium line is then used to extract the planetary sodium by division of the master-out, it creates an artificially deep residual because the emission varies in time and is not necessarily the same in and out of transit. This explains the unusually deep sodium feature in night two. This emission feature stems from telluric sodium, which was detected on fiber B in nights two and four.

To correct for the effect of the telluric sodium emission, we masked out the wavelength range of the emission feature in the stellar rest frame in each spectrum of nights two and four. The masking window was set at ±3.5 km s -1 of each line centre, based on the bin size of HARPS and the number of bins affected by the telluric sodium emission, see Figure 4. For WASP-127b, the planet signal in the stellar rest frame is shifted between -40 km s -1 and 20 km s -1 , with changes larger than one km s -1 from one exposure to the next. This implies that the masking will crop out part of the planetary sodium signal when it overlaps the masked region, in this case, in three to four exposures around mid-transit. However, we discuss Section 5.1.2 why spectra around mid-transit have to be rejected: they are affected by low S/N remnants. As a result, the masking for telluric sodium emission does not reduce the signal strength.

In Figure 5 we recreate the sodium transmission spectrum of Z19 (created from nights one and 2) at the left, but with our own telluric correction, and show the effect of masking the region that is affected by telluric sodium emission on their analysis at the right, with a significant decrease in the sodium feature. In total, the telluric sodium contamination makes up 43% of the signal strength in Z19, which reduces the detection levels from 8.33 to 4.8σ for the D1 and from 4.2 to 2.4σ for the D2 line. Out of all four nights, nights two and four were affected by telluric sodium emission and subsequently corrected in our analysis.

A stellar effect that might affect the transmission spectrum is the Rossiter-McLaughlin (RM) effect (Rossiter 1924;McLaughlin 1924;Queloz et al. 2000;Louden & Wheatley 2015;Cegla et al. 2016). This effect describes the impact of the Doppler-shift from the radial velocity of the local region of the stellar surface occulted by the planet at a given time of its transit. This veloc-Fig. 5. Percentage of extra absorption as shown in Z19, Figure 5 for the wavelength area of the sodium doublet. Grey points are at HARPS resolution, black points binned by 20. In the left panel is the transmission spectrum of WASP-127b created from the first two transits, that were available to the authors of Z19. On the right, the same dataset is taken and the core of the stellar sodium doublet masked out (we masked ±3.5 km s -1 of each line centre, amounting to roughly 8-9 bins per spectrum). In blue is a Gaussian fit on the grey data for visualisation. In line with Z19, we did not apply any spectra selection. The lower noise in comparison with Z19 can be attributed to our automatic correction for cosmics and different telluric corrections.

ity is below ∼ 2km s -1 for WASP-127. Wyttenbach et al. (2017) showed that the RM effect can be neglected for very slowly rotating stars. We thus concur with the reasoning of Z19 to neglect the RM effect in the analysis of the dataset presented here.

Stellar variability tracers

In addition to the telluric and stellar effects described before, stellar activity can affect the sodium feature. To detect any stellar effects, the Ca I and Mg I lines were monitored in Z19 (see Figure 3 in Z19), and no signal was found. Based on these results, Z19 assumed no effect of stellar activity in their analysis. However, we remark that these measurements have a low S/N, likely due to the inclusion of spectra that were taken while clouds passed or that were discarded for testing purposes (see Table 2 in Z19). Additionally, the plot of the Mg I line in Figure 3 of Z19 shows a spread of the data points at line centre greater than the noise level, with an uncertainty about twice as large at line centre than in the continuum. The high noise level in the Mg I line and the overall noise could effectively mask a signal and hide stellar activity. When we assume that no signal was hidden for Ca I and Mg I during the combined analysis of night one and 2, temporal stellar activity might still affect the transmission spectrum, but not be visible in the time-integrated spectrum. In this work, we present EulerCam light curves taken simultaneously with the spectroscopic data, see Figure 2, which show no photometric variability that might be attributable to stellar activity (e.g., spot modulations).

Low S/N trace

In each stellar spectrum, the most dominant feature is stellar sodium, where the flux drops significantly in the line core (see Figure 4 for an example). Because of the low flux, the photon Fig. 6. transmission spectrum of WASP-127b for the sodium doublet in the planetary rest frame. Upper panel: The transmission spectrum for all three (1,2,4) nights combined in full HARPS resolution is shown in grey, in black the same data is shown binned by x20. The data has been corrected for tellurics, cosmics, telluric sodium emission contamination and low S/N remnants. The rest frame transition wavelengths are marked with blue dashed lines. We attempted a Gaussian fit to both potential sodium lines jointly, however, the fit was only successful for the D2 line. The fit is shown as a red line. We, therefore, fitted the continuum to the D1 line instead. Lower panel: Residuals of the Gaussian fit.

noise is increased compared to the signal, resulting in a low S/N ratio in the wavelength region of the stellar sodium doublet. The planetary signal is extracted from the total flux by dividing each spectrum by the master-out in the stellar rest frame (SRF). In consequence, low S/N residuals remain at the former position of the stellar sodium doublet (Barnes et al. 2016;Borsa & Zannoni 2018;Seidel et al. 2020). When the spectra are then shifted to the planetary rest frame (PRF), the low S/N residuals are also shifted to different positions for each exposure (see Appendix A for 2D maps of the transits). However, for part of the in-transit spectra, the low S/N residuals lie in the wavelength range of the potential planetary sodium. In these regions, the residuals can mask a sodium signal due to the increase in noise or even produce a false sodium signature. We mitigate this effect by identifying the overlap between the low-S/N core of the stellar line (the FWHM of the stellar sodium feature at the wavelength sodium is known to absorb light) and the spectral range that potentially contains the planet sodium signal in each exposure. For WASP-127b this corresponds to a phase range from -0.008 to 0.008 in which this overlap would dominate any potential sodium signal. We subsequently discarded the affected exposures, see Table 1. A closer look at the 2D maps in Appendix A shows that the trace in night one and especially in night two is mainly in absorption and could thus introduce a false signal. This would account for the additional difference between the analysis in Z19 (barring the telluric sodium contamination) and our work. Night four is largely unaffected by the low S/N trace because the S/N is higher overall.

Sodium in WASP-127b

The resulting transmission spectrum for the three nights combined can be found in Figure 6. Using the Gaussian fit of the D2 line to estimate the detection level (see Hoeijmakers et al. 2020, for more information) results in a signal for the D2 line of 0.729 ± 0.454%, corresponding to 1.6σ. Although an absorption can be seen at the line centre of the D2 line, it was not detected in the D1 feature, and thus fitting Gaussians to both lines simultaneously was not possible. To calculate the significance of the overall sodium detection despite this, we compared the absorption in a 12Å window around the position of the sodium doublet (dashed blue lines in Figure 6), with the absorption in reference bands outside the sodium feature, as introduced in Wyttenbach et al. (2015). In the dataset presented here, the absorption in the sodium band is 0.456 ± 0.198%, corresponding to a significance of 2.3σ.

However, splitting the transmission spectrum into each night separately (see Figure 7) highlights that this feature stems largely from night one, with some absorption at the D2 line position in night two, but no indication of a sodium feature in night four. To quantify this assessment, we performed a Monte Carlo bootstrapping analysis based on Redfield et al. (2008). If the sodium detection is indeed produced by the planetary atmosphere and is not a spurious occurrence, no random selection of the spectra for a random in-transit and out-of-transit group of spectra should be able to reproduce the results. We created three different scenarios to test this hypothesis. An 'in-in' and 'out-out' scenario, where the random in-transit and out-of-transit spectra were all taken only from the real in-transit spectra or were all taken only from the real out-of-transit spectra, respectively. In the control case, 'in-out', the spectra were still drawn randomly, but for the random in-transit spectra only from the real in-transit spectra and vice versa.

The relative absorption depth distributions for the three analysed transits can be found in Figure 8. Both 'in-in' and 'out-out' are centred at 0.0%, as expected, meaning that any sodium detection cannot be reproduced solely from either the in-transit spectra or the out-of-transit spectra. However, the 'out-out' distribution has a large spread for night one (Figure 8), which results in a false-positive likelihood for the sodium detection in night one of 2.5% (Redfield et al. 2008). The large spread of the 'in-out' sample in the same night highlights the high noise in this specific night, which allows for a wide range of possibilities. No detection could be confirmed for night four.

Our results therefore show that a detection of sodium from this HARPS dataset can neither be confirmed nor confidently ruled out. It also remains unclear whether our inability to fit the D1 line stems only from the low S/N of the spectrum or if a physical mechanism reduces the strength of the D1 line disproportionally. Assuming that there is a sodium feature for WASP-127b, we propose a similar line of thought as for the ultra-hot Jupiter WASP-121b, a planet with a confirmed sodium detection and an interesting line ratio between the D2 and D1 line (Hoeijmakers et al. 2020): If sodium exists in the atmosphere of the planet and is extended enough, yet contained in the Roche radius, and additionally, if the number density of sodium is low, the outflow of sodium is optically thin in the line cores of the sodium doublet. This scenario is likely for WASP-127b, given that it is the puffiest exoplanet found to date. As shown in [START_REF] Draine | Physics of the Interstellar and Intergalactic Medium Dressing[END_REF], the line depth ratio of the D1 and D2 line on the ratio of the oscillator strengths in the thin limit, which is approximately 2 for the sodium doublet. Gebek & Oza (2020) linked the theoretical findings from [START_REF] Draine | Physics of the Interstellar and Intergalactic Medium Dressing[END_REF] and the puffy nature of some exoplanets and showed that when optically thin sodium gas is assumed across both line cores, the line ratio of the sodium doublet approaches 2 (see also Appendix A in Hoeijmakers et al. 2020). This hypothetical, optically thin, toroidal envelope of sodium around the planet could produce the line ratio needed to detect the D2 line, but might obscure the D1 line for WASP-127b.

Comparison with previous work

Placing our result in context with previous studies of WASP-127b highlights the need for further investigation of a potential sodium signature. WASP-127b has been observed multiple times at low resolution with ground-and space-based facilities. The first of these studies detected a strong Rayleigh slope and features compatible with TiO and VO in a most likely cloud-free atmosphere (Palle et al. 2017). However, sodium was detected tentatively in Palle et al. (2017) and was then claimed as a detection at the 5 σ level in the GTC follow-up Chen et al. (2018). Chen et al. (2018) detected potassium and lithium in addition to sodium, and the sodium detection was used to constrain the sodium abundance to super-solar. More recently, Spake et al. (2019) studied WASP-127b using the Hubble and Spitzer Space Telescopes and detected H 2 O, CO 2 , while additionally confirming the sodium detection from the first two studies with a supersolar abundance, despite an overall weaker signature in comparison, which might be due to the different resolutions of the instruments.

Lastly, as discussed in Section 5.1, Z19 presented a sodium detection based on the first two nights of this dataset and claimed a confirmation of the results presented with GTC instruments (Palle et al. 2017;Chen et al. 2018). We re-sampled our analysis on the HST/STIS wavelength grid from Figure 21 in Spake et al. (2019) showing all previously taken low-resolution datasets in R p /R . In the bin containing the sodium doublet, our data correspond to 0.1033 ± 0.0009R p /R , which is lower than the sodium detection in Spake et al. (2019) or Palle et al. (2017) and Chen et al. (2018) (see Figure 21 in Spake et al. (2019)), but they are compatible at the 1 σ level.

Conclusion

We studied the puffy super-Neptune WASP-127b in photometry and spectroscopy. We analysed seven transits from EulerCam and four transits from TESS and redetermined the stellar parameters from spectroscopy. The obtained updated system parameters were subsequently used to create the transmission spectrum. The spectroscopic data of four transits of WASP-127b were obtained with the HARPS spectrograph during the HEARTS survey, of which one transit was rejected because the out-of-transit baseline was insufficient. We showed that the previous analysis published in Z19 of two out of the four nights from this dataset resulted in a false-positive detection of neutral sodium based on contamination by stellar sodium and low S/N residuals. Our analysis of three nights of the same dataset showed an overall absorption in the 12Å band around the position of the sodium doublet at the 0.456 ± 0.198% level, corresponding to 2.3σ. This result is compatible with previous work conducted with space-and ground-based instruments, but it remains a nondetection. Curiously, we were able to fit a Gaussian to the D2 line, but not to the D1 line. This could be due to low S/N of the data or alternatively, if sodium exists in the atmosphere, due to a toroidal envelope of sodium in the optically thin limit. This curious finding shows the power of high-resolution spectroscopy in studying line shapes and ratios that would remain unnoticed at low resolution. We set the 3σ limit for the amplitude of a sodium detection with high-resolution spectrographs at 0.594% for future studies of the existence and potential ratio of the sodium lines. WASP-127b is a challenging target and a cautionary tale regarding the effect of low S/N remnants and telluric contamination. Any future analysis of this target has to carefully account for these complications.

Fig. 1 .

 1 Fig. 1. Exoplanets in the mass vs insolation space. The puffy super-Neptune WASP-127b shown as a big green dot. The data was retrieved from the NASA Exoplanet Archive: https://exoplanetarchive. ipac.caltech.edu/, retrieved 27 June 2020.

  , and for examples in the wavelength range of the sodium doublet, see Casasayas-Barris et al. (2019); Seidel et al. (2019); Chen et al. (2020); Casasayas-Barris et al. (2020); Hoeijmakers et al. (2020) or Seidel et al. (2020).

Fig. 2 .

 2 Fig. 2. Light curves obtained by EulerCam (bottom to top: 2014 April 28, 2016 March 09, 2017 February 28, 2017 March 20, 2018 February 13, 2018 March 31, 2019 February 24 offset by 0.02 for visibility. The computed models are shown in red.

Fig. 3 .

 3 Fig. 3. Light curves obtained by EulerCam (red model fit) and TESS (green model from detrended data and best-fit). The residuals from fitting are shown in the bottom, respectively.

Fig. 4 .

 4 Fig. 4. Normalised flux of the master out in the wavelength range of the sodium doublet. The stellar sodium doublet of night one is shown in blue, of night two in orange. The inset shows a zoom into the line centre of the D2 line of the same data set. Night one shows the expected line shape of the stellar sodium doublet, however, night two shows a sodium emission at line centre (width approx. 6-8 bins), which has to be masked in each spectrum to avoid a false-positive detection of the stellar sodium feature as planetary sodium.

Fig. 7 .

 7 Fig. 7. HARPS transmission spectrum of WASP-127b for the sodium doublet in the planetary rest frame for all three nights separately. See Figure 6 for details.

Fig. 8 .

 8 Fig.8. Distributions of the bootstrapping analysis for the 12 Å passband for the three used transits (night one, two, four) for 20000 random selections. The 'in-in' (magenta) and 'out-out' (light blue) distributions are centred at zero (no planetary sodium detection) and the 'in-out' distribution is shown in black.

Fig. A. 2 .

 2 Fig. A.2. Analysis of the overlap of the stellar and planetary sodium features for night two. See Figure A.1 for further details. The low S/N residual trace is visible for both the D2 and D1 line.

Fig. A. 3 .

 3 Fig. A.3. Analysis of the overlap of the stellar and planetary sodium features for night four. See Figure A.1 for further details. The low S/N residual trace is weak for this night.

Table 1 .

 1 Observation log.

	Article number, page 1 of 11

Table 2 .

 2 Overview of the EulerCam photometric observations.

	date	filter	RMS 2min [ppm] noise added
	2014 April 28	r'-Gunn 0.000959	0.000625
	2016 March 09	r'-Gunn 0.000910	0.000812
	2017 February 28 r'-Gunn 0.000704	0.000733
	2017 March 20	r'-Gunn 0.000721	0.000590
	2018 February 13 r'-Gunn 0.001069	0.000953
	2018 March 31	r'-Gunn 0.000946	0.000769
	2019 February 24 r'-Gunn 0.000946	0.000709

Table 3 .

 3 Planetary and stellar parameters obtained from EulerCam and TESS photometry.

	Jump parameters:
	Mid-transit time,[BJD] -2450000 6776.62124 +0.00023 -0.00028
	R p /R *	0.10103 0.00026 -0.00047
	a/R *	7.81 +0.11 -0.09
	Period [d]	4.17806203 0.00000088 -0.00000053
	RV amplitude, K [km s -1 ]	0.022 0.003 -0.002
	Primary parameters from spectral analysis:
	Stellar radius, R * [R ]	1.335 +0.025 -0.029
	Stellar mass, M * [M ]	0.949 +0.022 -0.019
	Derived parameters:
	Planetary radius, R p [R J ]	1.311 +0.025 -0.029
	Planetary mass, M p [M J ]	0.1647 +0.0214 -0.0172
	Planetary mean density, ρ p [ρ J ]	0.073 +0.010 -0.009
	Transit duration [d]	0.18137 +0.00035 -0.00058
	Transit depth, ∆F	0.01021 +0.00005 -0.00010
	Impact parameter	0.29 ± 0.04
	a/R *	7.81 +0.11 -0.09
	Orbital semi-major axis, a [au]	0.04840 +0.00136 -0.00095
	Stellar mean density, ρ * [ρ ]	0.401 +0.025 -0.026
	Stellar radius, R * [R ]	1.333 ± 0.027
	Stellar mass, M * [M ]	0.950 ± 0.020
	Inclination [deg]	87.84 +0.36 -0.33
	Eccentricity, e (fixed)	0.0
	Fixed quadratic limb-darkening parameters:
	u 1,r	0.422
	u 2,r	0.214
	u 1,t	0.350
	u 2,t	0.225
	c 1,r = 2u 1,r + u 2,r	1.057
	c 2,r = u 1,r -2u 2,r	0.208
	c 1,t = 2u 1,t + u 2,t	0.925
	c 2,t = u 1,t -2u 2,t	0.125
	the 8.33σ (D2) and 4.2σ level (D1) (see Figure

The last version of ARES code (ARES v2) can be downloaded at http://www.astro.up.pt/∼sousasag/ares

In Z19, the D1 and D2 lines have switched labels, which we assume is a typo.

Acknowledgements. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (project Four Aces; grant agreement No. 724427). This work has been carried out within the frame of the National Centre for Competence in Research 'PlanetS' supported by the Swiss National Science Foundation (SNSF). FPE and CLO would like to acknowledge the SNSF for supporting research with HARPS and EULER through the SNSF grants nr. 140649, 152721, 166227 and 184618. This work was supported by FCT -Fundação para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020 -Programa Operacional Competitividade e Internacionalização by the following grants: UID/FIS/04434/2019; UIDB/04434/2020; UIDP/04434/2020; PTDC/FIS-AST/32113/2017 & POCI-01-0145-FEDER-032113; PTDC/FIS-AST/28953/2017 & POCI-01-0145-FEDER-028953. A.W. acknowledges the financial support of the SNSF by grant number P400P2_186765. HMC acknowledges financial support from a UK Research and Innovation Future Leaders Fellowship. N. A.-D. acknowledges the support of FONDECYT project 3180063.

A. G. acknowledges support from the French CNRS and from the Agence Nationale de la Recherche (ANR grant GIPSE ANR-14-CE33-0018). A. G. has been supported by a grant from Labex OSUG@2020 (Investissements d'avenir -ANR10 LABX56). We thank L. A. dos Santos for his comments and insights during the preparation of this work.

Appendix A: Low S/N residuals